feat(reconstruction): Enhance volume reconstruction with geometry extraction and thickness merge options

atlas/templates/atlas/series_viewer.html

@@ -250,16 +250,11 @@
                                     <div class="card-header fw-semibold">Viewer / Compare Preview</div>
                                     <div class="card-body p-2 d-flex flex-column gap-2">
                                         <div id="truncate_viewer" class="dicom-viewer-root w-100"
-                                             style="height: 320px; background: #222;"
+                                             style="height: 460px; min-height: 360px; resize: vertical; overflow: auto; background: #222;"
                                              data-images="{{ image_url_array_and_count.0 }}" 
                                              data-auto-cache-stack="false">
                                         </div>
-                                        <div id="downsample-compare-wrap" class="d-none">
+                                        <div class="small text-muted">Downsample compare preview loads into this viewer as stack options.</div>
-                                            <div class="small text-muted mb-1">Downsample preview comparison</div>
-                                            <div id="downsample-compare-viewer" class="dicom-viewer-root w-100"
-                                                 style="height: 240px; background: #222;"
-                                                 data-auto-cache-stack="false"></div>
-                                        </div>
                                     </div>
                                 </div>
                             </div>
@@ -350,6 +345,15 @@
                                                         <label class="form-label small mb-1" for="recon-slice-spacing">Slice spacing (mm)</label>
                                                         <input type="number" step="0.1" min="0.1" id="recon-slice-spacing" name="recon_slice_spacing" class="form-control form-control-sm" placeholder="e.g. 1.0">
                                                     </div>
+                                                    <div class="col-12">
+                                                        <label class="form-label small mb-1" for="recon-thickness-mode">Thickness merge mode</label>
+                                                        <select id="recon-thickness-mode" name="recon_thickness_mode" class="form-select form-select-sm">
+                                                            <option value="mean" selected>Mean (default)</option>
+                                                            <option value="max">Max intensity</option>
+                                                            <option value="min">Min intensity</option>
+                                                        </select>
+                                                        <small class="text-muted">Thickness controls slab merge depth. Spacing controls distance between generated slices.</small>
+                                                    </div>
                                                 </div>
 
                                                 <div class="mt-2">
@@ -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 modalViewer = getTruncateViewerElement();
-                    const compareViewer = document.getElementById('downsample-compare-viewer');
+                    if (!modalViewer) {
-                    if (!compareWrap || !compareViewer) {
                         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) {

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):
+def _extract_recon_geometry(dicom_items):
-    """Nearest-neighbor resample along z-axis to requested spacing."""
+    """Build sorted geometry context for reconstruction from (image, ds, arr) items."""
     import numpy as np
 
-    if target_z_spacing <= 0:
+    row_spacing, col_spacing, native_z_spacing = _infer_native_spacings(dicom_items)
-        target_z_spacing = native_z_spacing
+    first_ds = dicom_items[0][1]
 
-    ratio = float(native_z_spacing) / float(target_z_spacing)
+    row_dir = np.array([1.0, 0.0, 0.0], dtype=float)
-    target_slices = max(2, int(round(volume.shape[0] * ratio)))
+    col_dir = np.array([0.0, 1.0, 0.0], dtype=float)
-    src_idx = np.clip(
+
-        np.round(np.linspace(0, volume.shape[0] - 1, target_slices)).astype(int),
+    iop = getattr(first_ds, "ImageOrientationPatient", None)
-        0,
+    if iop and len(iop) >= 6:
-        volume.shape[0] - 1,
+        try:
-    )
+            row_dir = np.array([float(iop[0]), float(iop[1]), float(iop[2])], dtype=float)
-    return volume[src_idx, :, :]
+            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('<div class="alert alert-danger mb-0">Invalid thickness merge mode.</div>')
 
         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('<div class="alert alert-danger mb-0">Slice thickness and spacing must be numeric values.</div>')
 
+        if slice_thickness_val is not None and slice_thickness_val <= 0:
+            return HttpResponse('<div class="alert alert-danger mb-0">Slice thickness must be greater than 0.</div>')
+        if slice_spacing_val is not None and slice_spacing_val <= 0:
+            return HttpResponse('<div class="alert alert-danger mb-0">Slice spacing must be greater than 0.</div>')
+
         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('<div class="alert alert-warning mb-0">Not enough consistently-sized slices for reconstruction.</div>')
 
-        volume = np.stack([item[2] for item in dicom_items], axis=0)
+        geom = _extract_recon_geometry(dicom_items)
-        template_ds = dicom_items[0][1]
+        sorted_items = geom["sorted_items"]
-        native_row_spacing, native_col_spacing, native_z_spacing = _infer_native_spacings(dicom_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('<div class="alert alert-danger mb-0">Failed to build reconstruction slabs for the requested spacing/thickness.</div>')
+
+        if volume_for_recon.shape[0] < 1:
+            return HttpResponse('<div class="alert alert-warning mb-0">No reconstruction slices were generated.</div>')
+
+        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 = (
+                ds_new.SliceThickness = float(slab_thickness)
-                    float(slice_thickness_val)
+                ds_new.SpacingBetweenSlices = float(spacing_between_slices_out)
-                    if slice_thickness_val is not None
-                    else float(target_spacing)
-                )
-                ds_new.SpacingBetweenSlices = float(target_spacing)
 
                 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'<div class="alert alert-success mb-0">Created <strong>{len(created_series)}</strong> reconstruction series: {links}</div>'
+            (
+                '<div class="alert alert-success mb-0">'
+                f'Created <strong>{len(created_series)}</strong> reconstruction series '
+                f'(spacing={target_spacing:.2f}mm, thickness={slab_thickness:.2f}mm, mode={escape(recon_thickness_mode)}): {links}'
+                '</div>'
+            )
         )
 
     return HttpResponse('<div class="alert alert-danger mb-0">Unknown optimize operation.</div>')