diff --git a/atlas/templates/atlas/series_viewer.html b/atlas/templates/atlas/series_viewer.html index 2b08f208..a295998b 100755 --- a/atlas/templates/atlas/series_viewer.html +++ b/atlas/templates/atlas/series_viewer.html @@ -250,16 +250,11 @@
Viewer / Compare Preview
-
-
Downsample preview comparison
-
-
+
Downsample compare preview loads into this viewer as stack options.
@@ -350,6 +345,15 @@ +
+ + + Thickness controls slab merge depth. Spacing controls distance between generated slices. +
@@ -385,6 +389,7 @@ const seriesPk = {{ series.pk }}; const totalImages = {{ image_url_array_and_count.1 }}; const viewerHeightKey = `series_viewer_height_${seriesPk}`; + const truncateViewerHeightKey = `truncate_viewer_height_${seriesPk}`; // Restore viewer height from localStorage function restoreViewerHeight() { @@ -438,10 +443,46 @@ endHidden.value = endInput.value; } + function getTruncateViewerElement() { + return document.getElementById(truncateApiKey); + } + + function persistTruncateViewerHeight() { + const modalViewer = getTruncateViewerElement(); + if (!modalViewer) { + return; + } + const height = window.getComputedStyle(modalViewer).height; + localStorage.setItem(truncateViewerHeightKey, height); + } + + function applyTruncateViewerHeight(el) { + if (!el) { + return; + } + const saved = localStorage.getItem(truncateViewerHeightKey); + if (saved) { + el.style.height = saved; + } + } + + function observeTruncateViewerResize() { + const modalViewer = getTruncateViewerElement(); + if (!modalViewer) { + return; + } + applyTruncateViewerHeight(modalViewer); + const observer = new ResizeObserver(() => { + persistTruncateViewerHeight(); + }); + observer.observe(modalViewer); + } + // Modal optimization handlers const truncateModal = document.getElementById('truncate-series-modal'); if (truncateModal) { truncateModal.addEventListener('show.bs.modal', function () { + observeTruncateViewerResize(); setTimeout(() => { try { if (window.mountDicomViewers) { @@ -530,24 +571,27 @@ return; } - const compareWrap = document.getElementById('downsample-compare-wrap'); - const compareViewer = document.getElementById('downsample-compare-viewer'); - if (!compareWrap || !compareViewer) { + const modalViewer = getTruncateViewerElement(); + if (!modalViewer) { return; } - compareWrap.classList.remove('d-none'); + const existingHeight = window.getComputedStyle(modalViewer).height; - // Recreate element to force a clean mount with new stacks. + // Recreate the same modal viewer element to force a clean remount with compare stacks. const fresh = document.createElement('div'); - fresh.id = 'downsample-compare-viewer'; + fresh.id = truncateApiKey; fresh.className = 'dicom-viewer-root w-100'; - fresh.style.height = '240px'; + fresh.style.height = existingHeight; + fresh.style.minHeight = '360px'; + fresh.style.resize = 'vertical'; + fresh.style.overflow = 'auto'; fresh.style.background = '#222'; fresh.setAttribute('data-auto-cache-stack', 'false'); fresh.setAttribute('data-named-stacks', namedStacks); - compareViewer.replaceWith(fresh); + modalViewer.replaceWith(fresh); + observeTruncateViewerResize(); try { if (window.mountDicomViewers) { diff --git a/atlas/views.py b/atlas/views.py index f116ceb6..0390c853 100755 --- a/atlas/views.py +++ b/atlas/views.py @@ -376,21 +376,130 @@ def _infer_native_spacings(dicom_items): return float(row_spacing), float(col_spacing), float(z_spacing) -def _resample_volume_z_nearest(volume, native_z_spacing, target_z_spacing): - """Nearest-neighbor resample along z-axis to requested spacing.""" +def _extract_recon_geometry(dicom_items): + """Build sorted geometry context for reconstruction from (image, ds, arr) items.""" import numpy as np - if target_z_spacing <= 0: - target_z_spacing = native_z_spacing + row_spacing, col_spacing, native_z_spacing = _infer_native_spacings(dicom_items) + first_ds = dicom_items[0][1] - ratio = float(native_z_spacing) / float(target_z_spacing) - target_slices = max(2, int(round(volume.shape[0] * ratio))) - src_idx = np.clip( - np.round(np.linspace(0, volume.shape[0] - 1, target_slices)).astype(int), - 0, - volume.shape[0] - 1, - ) - return volume[src_idx, :, :] + row_dir = np.array([1.0, 0.0, 0.0], dtype=float) + col_dir = np.array([0.0, 1.0, 0.0], dtype=float) + + iop = getattr(first_ds, "ImageOrientationPatient", None) + if iop and len(iop) >= 6: + try: + row_dir = np.array([float(iop[0]), float(iop[1]), float(iop[2])], dtype=float) + col_dir = np.array([float(iop[3]), float(iop[4]), float(iop[5])], dtype=float) + except Exception: + row_dir = np.array([1.0, 0.0, 0.0], dtype=float) + col_dir = np.array([0.0, 1.0, 0.0], dtype=float) + + if np.linalg.norm(row_dir) == 0 or np.linalg.norm(col_dir) == 0: + row_dir = np.array([1.0, 0.0, 0.0], dtype=float) + col_dir = np.array([0.0, 1.0, 0.0], dtype=float) + + row_dir = row_dir / np.linalg.norm(row_dir) + col_dir = col_dir / np.linalg.norm(col_dir) + + normal_dir = np.cross(row_dir, col_dir) + if np.linalg.norm(normal_dir) == 0: + normal_dir = np.array([0.0, 0.0, 1.0], dtype=float) + normal_dir = normal_dir / np.linalg.norm(normal_dir) + + decorated = [] + for idx, item in enumerate(dicom_items): + ds = item[1] + ipp = None + try: + ipp_attr = getattr(ds, "ImagePositionPatient", None) + if ipp_attr and len(ipp_attr) >= 3: + ipp = np.array([float(ipp_attr[0]), float(ipp_attr[1]), float(ipp_attr[2])], dtype=float) + except Exception: + ipp = None + + if ipp is None: + ipp = np.array([0.0, 0.0, idx * native_z_spacing], dtype=float) + + proj = float(np.dot(ipp, normal_dir)) + decorated.append((proj, ipp, item)) + + decorated.sort(key=lambda x: x[0]) + sorted_items = [d[2] for d in decorated] + sorted_ipps = [d[1] for d in decorated] + sorted_proj = np.array([d[0] for d in decorated], dtype=float) + + if len(sorted_proj) > 1: + diffs = np.diff(sorted_proj) + valid_diffs = np.abs(diffs[np.abs(diffs) > 1e-6]) + if valid_diffs.size > 0: + native_z_spacing = float(np.median(np.abs(valid_diffs))) + + if native_z_spacing <= 0: + native_z_spacing = 1.0 + + origin_ipp = sorted_ipps[0] if sorted_ipps else np.array([0.0, 0.0, 0.0], dtype=float) + source_positions_mm = sorted_proj - sorted_proj[0] if len(sorted_proj) > 0 else np.array([0.0], dtype=float) + + return { + "sorted_items": sorted_items, + "source_positions_mm": source_positions_mm, + "origin_ipp": origin_ipp, + "row_dir": row_dir, + "col_dir": col_dir, + "normal_dir": normal_dir, + "row_spacing": float(row_spacing), + "col_spacing": float(col_spacing), + "native_z_spacing": float(native_z_spacing), + } + + +def _aggregate_volume_along_z(volume, source_positions_mm, target_spacing, slab_thickness, mode): + """Aggregate a z-volume into slabbed slices using spacing/thickness in mm.""" + import numpy as np + + if target_spacing <= 0: + raise ValueError("Target spacing must be positive") + if slab_thickness <= 0: + raise ValueError("Slab thickness must be positive") + + source_positions_mm = np.asarray(source_positions_mm, dtype=float) + if source_positions_mm.size != volume.shape[0]: + raise ValueError("Source position count must match volume depth") + + start_pos = float(source_positions_mm[0]) + end_pos = float(source_positions_mm[-1]) + centers = np.arange(start_pos, end_pos + (target_spacing * 0.5), target_spacing, dtype=float) + if centers.size == 0: + centers = np.array([start_pos], dtype=float) + + half = slab_thickness / 2.0 + slabs = [] + for center in centers: + mask = np.abs(source_positions_mm - center) <= (half + 1e-6) + if not mask.any(): + nearest = int(np.argmin(np.abs(source_positions_mm - center))) + slab = volume[nearest : nearest + 1, :, :] + else: + slab = volume[mask, :, :] + + if mode == "max": + out = np.max(slab, axis=0) + elif mode == "min": + out = np.min(slab, axis=0) + else: + out = np.mean(slab, axis=0) + slabs.append(out) + + out_volume = np.stack(slabs, axis=0) + + if np.issubdtype(volume.dtype, np.integer): + info = np.iinfo(volume.dtype) + out_volume = np.clip(np.rint(out_volume), info.min, info.max).astype(volume.dtype) + else: + out_volume = out_volume.astype(volume.dtype, copy=False) + + return out_volume, centers @login_required @@ -548,6 +657,10 @@ def series_optimize_htmx(request, series_id): slice_thickness = request.POST.get("recon_slice_thickness", "") slice_spacing = request.POST.get("recon_slice_spacing", "") + recon_thickness_mode = (request.POST.get("recon_thickness_mode", "mean") or "mean").strip().lower() + + if recon_thickness_mode not in ("mean", "max", "min"): + return HttpResponse('
Invalid thickness merge mode.
') try: slice_thickness_val = float(slice_thickness) if slice_thickness else None @@ -555,6 +668,11 @@ def series_optimize_htmx(request, series_id): except ValueError: return HttpResponse('
Slice thickness and spacing must be numeric values.
') + if slice_thickness_val is not None and slice_thickness_val <= 0: + return HttpResponse('
Slice thickness must be greater than 0.
') + if slice_spacing_val is not None and slice_spacing_val <= 0: + return HttpResponse('
Slice spacing must be greater than 0.
') + dicom_items = [] for image in bounded_images: ds = _read_series_image_dataset(image) @@ -578,12 +696,51 @@ def series_optimize_htmx(request, series_id): if len(dicom_items) < 2: return HttpResponse('
Not enough consistently-sized slices for reconstruction.
') - volume = np.stack([item[2] for item in dicom_items], axis=0) - template_ds = dicom_items[0][1] - native_row_spacing, native_col_spacing, native_z_spacing = _infer_native_spacings(dicom_items) + geom = _extract_recon_geometry(dicom_items) + sorted_items = geom["sorted_items"] + volume = np.stack([item[2] for item in sorted_items], axis=0) + + template_ds = sorted_items[0][1] + source_positions_mm = geom["source_positions_mm"] + origin_ipp = geom["origin_ipp"] + row_dir = geom["row_dir"] + col_dir = geom["col_dir"] + normal_dir = geom["normal_dir"] + native_row_spacing = geom["row_spacing"] + native_col_spacing = geom["col_spacing"] + native_z_spacing = geom["native_z_spacing"] target_spacing = slice_spacing_val if slice_spacing_val is not None else native_z_spacing - volume_for_recon = _resample_volume_z_nearest(volume, native_z_spacing, target_spacing) + slab_thickness = slice_thickness_val if slice_thickness_val is not None else target_spacing + + try: + volume_for_recon, recon_centers_mm = _aggregate_volume_along_z( + volume, + source_positions_mm, + float(target_spacing), + float(slab_thickness), + recon_thickness_mode, + ) + except Exception as exc: + logger.warning("Reconstruction aggregation failed for series {}: {}", series.pk, exc) + return HttpResponse('
Failed to build reconstruction slabs for the requested spacing/thickness.
') + + if volume_for_recon.shape[0] < 1: + return HttpResponse('
No reconstruction slices were generated.
') + + base_z_offset = float(recon_centers_mm[0]) if len(recon_centers_mm) > 0 else 0.0 + + def build_position(origin, row_vector, col_vector, stack_vector, row_index=0, col_index=0, stack_offset=0.0): + import numpy as np + + pos = ( + np.asarray(origin, dtype=float) + + (np.asarray(stack_vector, dtype=float) * float(stack_offset)) + + (np.asarray(row_vector, dtype=float) * (float(row_index) * float(native_row_spacing))) + + (np.asarray(col_vector, dtype=float) * (float(col_index) * float(native_col_spacing))) + ) + return [float(pos[0]), float(pos[1]), float(pos[2])] + created_series = [] for plane in recon_planes: @@ -591,12 +748,72 @@ def series_optimize_htmx(request, series_id): if plane_norm == "axial": recon_slices = [volume_for_recon[i, :, :] for i in range(volume_for_recon.shape[0])] pixel_spacing_out = [native_row_spacing, native_col_spacing] + spacing_between_slices_out = float(target_spacing) + image_orientation_out = [ + float(row_dir[0]), + float(row_dir[1]), + float(row_dir[2]), + float(col_dir[0]), + float(col_dir[1]), + float(col_dir[2]), + ] + + def position_for_index(idx): + return build_position( + origin_ipp, + row_dir, + col_dir, + normal_dir, + row_index=0, + col_index=0, + stack_offset=float(recon_centers_mm[idx]), + ) elif plane_norm == "coronal": recon_slices = [volume_for_recon[:, i, :] for i in range(volume_for_recon.shape[1])] pixel_spacing_out = [target_spacing, native_col_spacing] + spacing_between_slices_out = float(native_row_spacing) + image_orientation_out = [ + float(normal_dir[0]), + float(normal_dir[1]), + float(normal_dir[2]), + float(col_dir[0]), + float(col_dir[1]), + float(col_dir[2]), + ] + + def position_for_index(idx): + return build_position( + origin_ipp, + row_dir, + col_dir, + normal_dir, + row_index=idx, + col_index=0, + stack_offset=base_z_offset, + ) elif plane_norm == "sagittal": recon_slices = [volume_for_recon[:, :, i] for i in range(volume_for_recon.shape[2])] pixel_spacing_out = [target_spacing, native_row_spacing] + spacing_between_slices_out = float(native_col_spacing) + image_orientation_out = [ + float(normal_dir[0]), + float(normal_dir[1]), + float(normal_dir[2]), + float(row_dir[0]), + float(row_dir[1]), + float(row_dir[2]), + ] + + def position_for_index(idx): + return build_position( + origin_ipp, + row_dir, + col_dir, + normal_dir, + row_index=0, + col_index=idx, + stack_offset=base_z_offset, + ) else: continue @@ -615,13 +832,11 @@ def series_optimize_htmx(request, series_id): ds_new.SeriesInstanceUID = recon_series.series_instance_uid ds_new.PixelData = arr2d.tobytes() ds_new.PixelSpacing = [float(pixel_spacing_out[0]), float(pixel_spacing_out[1])] + ds_new.ImageOrientationPatient = image_orientation_out + ds_new.ImagePositionPatient = position_for_index(idx) - ds_new.SliceThickness = ( - float(slice_thickness_val) - if slice_thickness_val is not None - else float(target_spacing) - ) - ds_new.SpacingBetweenSlices = float(target_spacing) + ds_new.SliceThickness = float(slab_thickness) + ds_new.SpacingBetweenSlices = float(spacing_between_slices_out) out_io = io.BytesIO() ds_new.save_as(out_io, write_like_original=False) @@ -647,7 +862,12 @@ def series_optimize_htmx(request, series_id): ] ) return HttpResponse( - f'
Created {len(created_series)} reconstruction series: {links}
' + ( + '
' + f'Created {len(created_series)} reconstruction series ' + f'(spacing={target_spacing:.2f}mm, thickness={slab_thickness:.2f}mm, mode={escape(recon_thickness_mode)}): {links}' + '
' + ) ) return HttpResponse('
Unknown optimize operation.
')