Merge branch 'master' of https://github.com/tfunck/brainbuilder

Author: tfunck

brainbuilder/interp/acqvolume.py

@@ -93,7 +93,6 @@ def thicken_sections_within_chunk(
     tissue_type: str = "",
     gaussian_sd: float = 0,
     width: int = None,
-    use_adjust_section_means: bool = True,
 ) -> None:
     """Thicken sections within a chunk. A thickened section is simply a section that is expanded along the y axis to the resolution of the reconstruction.
 
@@ -114,13 +113,6 @@ def thicken_sections_within_chunk(
         "Error: source volume for thickening sections is empty\n" + source_image_fn
     )
 
-    if use_adjust_section_means:
-        adj_factors_dict = calculate_section_adjustment_factors(
-            chunk_sect_info["nl_2d_rsl"].values,
-            chunk_sect_info["sample"].values,
-            os.path.dirname(thickened_fn),
-        )
-
     if width is None:
         width = get_thicken_width(resolution, section_thickness, scale=1)
 
@@ -137,9 +129,6 @@ def thicken_sections_within_chunk(
         nl_2d_rsl = row["nl_2d_rsl"]
         section = nib.load(nl_2d_rsl).get_fdata().copy()
 
-        if use_adjust_section_means:
-            section += adj_factors_dict[y]
-
         if np.sum(section) == 0:
             print(f"Warning: empty frame {row_i} {row}\n")
 
@@ -153,6 +142,7 @@ def thicken_sections_within_chunk(
         rep = np.repeat(section.reshape(dim), y1 - y0, axis=1)
 
         rec_vol[:, y0:y1, :] += rep
+
         n[:, y0:y1, :] += 1
 
     # normalize by number of sections within range
@@ -393,7 +383,6 @@ def create_thickened_volumes(
                     chunk_sect_info,
                     resolution,
                     gaussian_sd=gaussian_sd,
-                    use_adjust_section_means=True,
                     width=width,
                 )
 

brainbuilder/interp/volinterp.py

@@ -11,7 +11,7 @@
 
 import brainbuilder.utils.ants_nibabel as nib
 from brainbuilder.interp.acqvolume import create_thickened_volumes
-from brainbuilder.utils.nl_deformation_flow import nl_deformation_flow_3d
+from brainbuilder.utils.nl_deformation_flow import nlflow_isometric
 from brainbuilder.volalign import align_3d, verify_chunk_limits
 
 
@@ -80,81 +80,6 @@ def _interpolate_missing_sections(
     return vol_interp
 
 
-def nl_deformation_flow_nii(acq_fin, output_dir, interp_acq_fin, clobber=False):
-    if not os.path.exists(interp_acq_fin) or clobber:
-        """
-        Process the acq volume by applying non-linear deformation flow.
-
-        Parameters
-        ----------
-        acq_fin : str
-            Path to the acq volume file.
-        output_dir : str
-            Directory to save the output files.
-        clobber : bool, optional
-            If True, overwrite existing files. The default is False.
-
-        Returns
-        -------
-        str
-            Path to the processed acq volume file.
-        """
-        acq_img = nib.load(acq_fin)
-        acq_vol = acq_img.get_fdata()
-        acq_vol[acq_vol < 0] = 0
-
-        origin = list(acq_img.affine[[0, 2], 3])
-        spacing = list(acq_img.affine[[0, 2], [0, 2]])
-
-        interp_acq_vol = nl_deformation_flow_3d(
-            acq_vol,
-            output_dir + "/nl_flow/",
-            origin=origin,
-            spacing=spacing,
-            clobber=clobber,
-        )
-
-        interp_acq_img = nib.Nifti1Image(
-            interp_acq_vol, acq_img.affine, direction_order="lpi"
-        )
-        interp_acq_img.to_filename(interp_acq_fin)
-
-
-def resample_interp_vol_to_resolution(
-    interp_acq_orig_fin,
-    interp_acq_iso_fin: str,
-    resolution: float,
-    clobber: bool = False,
-) -> np.array:
-    """Resample the interpolated volume to the specified resolution.
-    Parameters
-    ----------
-    interp_acq_vol : np.array
-        Interpolated volume.
-    acq_img : nib.Nifti1Image
-
-    """
-    if not os.path.exists(interp_acq_iso_fin) or clobber:
-        interp_acq_img = nib.load(interp_acq_orig_fin)  # type:ignore[assignment]
-        interp_acq_vol = interp_acq_img.get_fdata()
-
-        slice_thickness = interp_acq_img.affine[1, 1]
-
-        y_new = int(np.round(interp_acq_vol.shape[1] / (resolution / slice_thickness)))
-
-        interp_acq_vol = resize(
-            interp_acq_vol,
-            (interp_acq_vol.shape[0], y_new, interp_acq_vol.shape[2]),
-            order=1,
-        )
-
-        aff_iso = interp_acq_img.affine
-
-        aff_iso[1, 1] = resolution
-
-        nib.Nifti1Image(interp_acq_vol, aff_iso, direction_order="lpi").to_filename(
-            interp_acq_iso_fin
-        )
 
 
 def volumetric_interpolation(
@@ -210,20 +135,15 @@ def volumetric_interpolation(
         chunk_info_thickened = pd.read_csv(chunk_info_thickened_csv)
 
         print(f"Interpolation Chunk: {chunk}, Acquisition: {acq}")
+        cls_fin = chunk_info_thickened["thickened_cls"].values[0]
         acq_fin = chunk_info_thickened["thickened"].values[0]
 
-        interp_acq_iso_fin = acq_fin.replace("thickened", "interp-vol_iso")
-        interp_acq_fin = acq_fin.replace("thickened", "interp-vol_orig")
+        interp_acq_iso_fin, nlflow_tfm_list = nlflow_isometric(acq_fin, curr_output_dir, resolution, clobber=clobber)
 
-        print("Resampling", interp_acq_fin)
-        nl_deformation_flow_nii(
-            acq_fin, curr_output_dir, interp_acq_fin, clobber=clobber
-        )
+        interp_cls_iso_fin = nlflow_isometric(
+            cls_fin, curr_output_dir, resolution, tfm_list = nlflow_tfm_list, clobber=clobber
+        ) 
 
-        print("Resampling", interp_acq_iso_fin)
-        resample_interp_vol_to_resolution(
-            interp_acq_fin, interp_acq_iso_fin, resolution, clobber=clobber
-        )
 
         row = pd.DataFrame(
             {
@@ -232,6 +152,7 @@ def volumetric_interpolation(
                 "chunk": [chunk],
                 "acquisition": [acq],
                 "interp_nat": [interp_acq_iso_fin],
+                "interp_cls_nat": [interp_cls_iso_fin],
             }
         )
 

brainbuilder/qc/validate_interp_error.py

@@ -685,7 +685,7 @@ def plot_validation_interp_error(
     print(np.mean(df["error_2d"]), "+/-", np.std(df["error_2d"]))
     print(np.mean(df["error"]), "+/-", np.std(df["error"]))
     plt.suptitle("Validation of 2D & 3D Interpolation Error", fontsize=16)
-    plt.savefig(f"{output_dir}/validation_interp_error.png")
+    plt.savefig(f"{output_dir}/validation_interp_error.png", dpi=300)
     print(f"{output_dir}/validation_interp_error.png")
     plt.close()
     plt.cla()
@@ -706,7 +706,7 @@ def plot_validation_interp_error(
     # plt.ylabel('Interpolation Error %')
     # plt.xlabel('Volume Change %')
 
-    plt.savefig(f"{output_dir}/validation_interp_error_vs_volume.png")
+    plt.savefig(f"{output_dir}/validation_interp_error_vs_volume.png", dpi=300)
     print(f"{output_dir}/validation_interp_error_vs_volume.png")
 
     # plt.subplot(2, 2, 4)