diff --git a/atlas/tasks.py b/atlas/tasks.py index afe62320..18f43b90 100644 --- a/atlas/tasks.py +++ b/atlas/tasks.py @@ -199,19 +199,8 @@ def series_reconstruct_task( native_col_spacing = geom["col_spacing"] native_z_spacing = geom["native_z_spacing"] - target_spacing = float(slice_spacing_val) if slice_spacing_val is not None else native_z_spacing - slab_thickness = float(slice_thickness_val) if slice_thickness_val is not None else target_spacing - - volume_for_recon, recon_centers_mm = atlas_views._aggregate_volume_along_z( - volume, - source_positions_mm, - float(target_spacing), - float(slab_thickness), - recon_thickness_mode, - ) - - if volume_for_recon.shape[0] < 1: - raise ValueError("No reconstruction slices were generated") + requested_slice_spacing = float(slice_spacing_val) if slice_spacing_val is not None else None + requested_slab_thickness = float(slice_thickness_val) if slice_thickness_val is not None else None # Determine source acquisition plane from normal vector. dominant_normal_axis = int(np.argmax(np.abs(normal_dir))) @@ -219,9 +208,14 @@ def series_reconstruct_task( # Reorient the reconstructed volume into patient axes order [z, y, x] so # generated plane names are anatomically correct regardless of source plane. - source_axis_dirs = [normal_dir, row_dir, col_dir] - source_axis_spacings = [float(target_spacing), float(native_row_spacing), float(native_col_spacing)] - source_axis_sizes = [int(volume_for_recon.shape[0]), int(volume_for_recon.shape[1]), int(volume_for_recon.shape[2])] + # NOTE: DICOM stores ImageOrientationPatient as: + # - first triplet: direction across columns (x-axis in pixel grid) + # - second triplet: direction across rows (y-axis in pixel grid) + # In _extract_recon_geometry names are historically row_dir/col_dir, so we + # map axis spacing/indexing explicitly to avoid directional mixups. + source_axis_dirs = [normal_dir, col_dir, row_dir] + source_axis_spacings = [float(native_z_spacing), float(native_row_spacing), float(native_col_spacing)] + source_axis_sizes = [int(volume.shape[0]), int(volume.shape[1]), int(volume.shape[2])] source_for_patient_axis = {} sign_for_patient_axis = {} @@ -243,7 +237,7 @@ def series_reconstruct_task( spacing_y = source_axis_spacings[src_y] spacing_z = source_axis_spacings[src_z] - vol_zyx = np.transpose(volume_for_recon, (src_z, src_y, src_x)) + vol_zyx = np.transpose(volume, (src_z, src_y, src_x)) if sign_for_patient_axis[2] < 0: vol_zyx = np.flip(vol_zyx, axis=0) @@ -255,9 +249,9 @@ def series_reconstruct_task( def source_indices_from_patient_zyx(iz, iy, ix): src_indices = [0, 0, 0] - z_index = int(iz) - y_index = int(iy) - x_index = int(ix) + z_index = float(iz) + y_index = float(iy) + x_index = float(ix) if sign_for_patient_axis[2] < 0: z_index = source_axis_sizes[src_z] - 1 - z_index @@ -275,41 +269,106 @@ def series_reconstruct_task( src_k, src_r, src_c = source_indices_from_patient_zyx(iz, iy, ix) pos = ( np.asarray(origin_ipp, dtype=float) - + (np.asarray(normal_dir, dtype=float) * (float(src_k) * float(target_spacing))) - + (np.asarray(row_dir, dtype=float) * (float(src_r) * float(native_row_spacing))) - + (np.asarray(col_dir, dtype=float) * (float(src_c) * float(native_col_spacing))) + + (np.asarray(normal_dir, dtype=float) * (float(src_k) * float(native_z_spacing))) + + (np.asarray(col_dir, dtype=float) * (float(src_r) * float(native_row_spacing))) + + (np.asarray(row_dir, dtype=float) * (float(src_c) * float(native_col_spacing))) ) return [float(pos[0]), float(pos[1]), float(pos[2])] + spacing_x = source_axis_spacings[src_x] + spacing_y = source_axis_spacings[src_y] + spacing_z = source_axis_spacings[src_z] + + def _reduce_slab(slab): + if recon_thickness_mode == "max": + out = np.max(slab, axis=0) + elif recon_thickness_mode == "min": + out = np.min(slab, axis=0) + else: + out = np.mean(slab, axis=0) + + if np.issubdtype(volume.dtype, np.integer): + info = np.iinfo(volume.dtype) + out = np.clip(np.rint(out), info.min, info.max).astype(volume.dtype) + else: + out = out.astype(volume.dtype, copy=False) + return out + + def _slice_indices_for_axis(axis_length, axis_spacing, output_spacing, slab_thickness): + if output_spacing <= 0: + raise ValueError("Slice spacing must be greater than 0") + if slab_thickness <= 0: + raise ValueError("Slice thickness must be greater than 0") + + axis_positions_mm = np.arange(axis_length, dtype=float) * float(axis_spacing) + end_pos = axis_positions_mm[-1] if axis_positions_mm.size else 0.0 + centers_mm = np.arange(0.0, end_pos + (output_spacing * 0.5), output_spacing, dtype=float) + if centers_mm.size == 0: + centers_mm = np.array([0.0], dtype=float) + + half = slab_thickness / 2.0 + slabs = [] + centers_idx = [] + for center_mm in centers_mm: + mask = np.abs(axis_positions_mm - center_mm) <= (half + 1e-6) + if not mask.any(): + nearest = int(np.argmin(np.abs(axis_positions_mm - center_mm))) + idxs = np.array([nearest], dtype=int) + else: + idxs = np.where(mask)[0].astype(int) + + slabs.append(idxs) + centers_idx.append(float(center_mm / axis_spacing) if axis_spacing > 0 else float(idxs[0])) + + return slabs, centers_idx + plane_builders = {} for plane in recon_planes: plane_norm = plane.lower() if plane_norm == "axial": + plane_spacing = float(requested_slice_spacing) if requested_slice_spacing is not None else float(spacing_z) + plane_thickness = float(requested_slab_thickness) if requested_slab_thickness is not None else float(plane_spacing) + slab_indices, slab_centers_idx = _slice_indices_for_axis( + int(vol_zyx.shape[0]), float(spacing_z), plane_spacing, plane_thickness + ) plane_builders[plane_norm] = { - "count": int(vol_zyx.shape[0]), - "slice_fn": lambda idx: vol_zyx[idx, :, :], + "count": int(len(slab_indices)), + "slice_fn": lambda idx, slabs=slab_indices: _reduce_slab(vol_zyx[slabs[idx], :, :]), "pixel_spacing_out": [float(spacing_y), float(spacing_x)], - "spacing_between_slices_out": float(spacing_z), - "image_orientation_out": [0.0, 1.0, 0.0, 1.0, 0.0, 0.0], - "position_fn": lambda idx: position_from_patient_zyx(idx, 0, 0), + "spacing_between_slices_out": float(plane_spacing), + "slice_thickness_out": float(plane_thickness), + "image_orientation_out": [1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + "position_fn": lambda idx, centers=slab_centers_idx: position_from_patient_zyx(centers[idx], 0.0, 0.0), } elif plane_norm == "coronal": + plane_spacing = float(requested_slice_spacing) if requested_slice_spacing is not None else float(spacing_y) + plane_thickness = float(requested_slab_thickness) if requested_slab_thickness is not None else float(plane_spacing) + slab_indices, slab_centers_idx = _slice_indices_for_axis( + int(vol_zyx.shape[1]), float(spacing_y), plane_spacing, plane_thickness + ) plane_builders[plane_norm] = { - "count": int(vol_zyx.shape[1]), - "slice_fn": lambda idx: vol_zyx[:, idx, :], + "count": int(len(slab_indices)), + "slice_fn": lambda idx, slabs=slab_indices: _reduce_slab(vol_zyx[:, slabs[idx], :]), "pixel_spacing_out": [float(spacing_z), float(spacing_x)], - "spacing_between_slices_out": float(spacing_y), - "image_orientation_out": [0.0, 0.0, 1.0, 1.0, 0.0, 0.0], - "position_fn": lambda idx: position_from_patient_zyx(0, idx, 0), + "spacing_between_slices_out": float(plane_spacing), + "slice_thickness_out": float(plane_thickness), + "image_orientation_out": [1.0, 0.0, 0.0, 0.0, 0.0, 1.0], + "position_fn": lambda idx, centers=slab_centers_idx: position_from_patient_zyx(0.0, centers[idx], 0.0), } elif plane_norm == "sagittal": + plane_spacing = float(requested_slice_spacing) if requested_slice_spacing is not None else float(spacing_x) + plane_thickness = float(requested_slab_thickness) if requested_slab_thickness is not None else float(plane_spacing) + slab_indices, slab_centers_idx = _slice_indices_for_axis( + int(vol_zyx.shape[2]), float(spacing_x), plane_spacing, plane_thickness + ) plane_builders[plane_norm] = { - "count": int(vol_zyx.shape[2]), - "slice_fn": lambda idx: vol_zyx[:, :, idx], + "count": int(len(slab_indices)), + "slice_fn": lambda idx, slabs=slab_indices: _reduce_slab(vol_zyx[:, :, slabs[idx]]), "pixel_spacing_out": [float(spacing_z), float(spacing_y)], - "spacing_between_slices_out": float(spacing_x), - "image_orientation_out": [0.0, 0.0, 1.0, 0.0, 1.0, 0.0], - "position_fn": lambda idx: position_from_patient_zyx(0, 0, idx), + "spacing_between_slices_out": float(plane_spacing), + "slice_thickness_out": float(plane_thickness), + "image_orientation_out": [0.0, 1.0, 0.0, 0.0, 0.0, 1.0], + "position_fn": lambda idx, centers=slab_centers_idx: position_from_patient_zyx(0.0, 0.0, centers[idx]), } if not plane_builders: @@ -344,7 +403,7 @@ def series_reconstruct_task( ds_new.PixelSpacing = [float(cfg["pixel_spacing_out"][0]), float(cfg["pixel_spacing_out"][1])] ds_new.ImageOrientationPatient = cfg["image_orientation_out"] ds_new.ImagePositionPatient = cfg["position_fn"](idx) - ds_new.SliceThickness = float(slab_thickness) + ds_new.SliceThickness = float(cfg["slice_thickness_out"]) ds_new.SpacingBetweenSlices = float(cfg["spacing_between_slices_out"]) out_io = io.BytesIO() @@ -401,8 +460,8 @@ def series_reconstruct_task( return { "series_id": series.pk, "created_series": created_series, - "target_spacing": float(target_spacing), - "slab_thickness": float(slab_thickness), + "target_spacing": float(requested_slice_spacing) if requested_slice_spacing is not None else None, + "slab_thickness": float(requested_slab_thickness) if requested_slab_thickness is not None else None, "mode": recon_thickness_mode, "source_plane": source_plane, } \ No newline at end of file diff --git a/atlas/templates/atlas/task_overview.html b/atlas/templates/atlas/task_overview.html index 5cb97242..48162281 100644 --- a/atlas/templates/atlas/task_overview.html +++ b/atlas/templates/atlas/task_overview.html @@ -51,10 +51,28 @@ Atlas Task Overview -