fix: do not fit intercept, improve tests

Author: oesteban

sdcflows/interfaces/bspline.py

@@ -140,15 +140,15 @@ def _run_interface(self, runtime):
 
         # Load in the fieldmap
         fmapnii = nb.load(self.inputs.in_data)
-        fmapnii = nb.as_closest_canonical(fmapnii)
+        # fmapnii = nb.as_closest_canonical(fmapnii)
         zooms = fmapnii.header.get_zooms()
 
         # Get a mask (or define on the spot to cover the full extent)
         masknii = (
             nb.load(self.inputs.in_mask) if isdefined(self.inputs.in_mask) else None
         )
-        if masknii is not None:
-            masknii = nb.as_closest_canonical(masknii)
+        # if masknii is not None:
+        #     masknii = nb.as_closest_canonical(masknii)
 
         # Determine the shape of bspline coefficients
         # This should not change with resizing, so do it first
@@ -212,7 +212,7 @@ def _run_interface(self, runtime):
         )
 
         # Fit the model
-        model = lm.Ridge(alpha=self.inputs.ridge_alpha, fit_intercept=True)
+        model = lm.Ridge(alpha=self.inputs.ridge_alpha, fit_intercept=False)
         for attempt in range(3):
             model.fit(colmat, data.reshape(-1))
             extreme = np.abs(model.coef_).max()
@@ -248,11 +248,11 @@ def _run_interface(self, runtime):
         # Interpolating in the original grid will require a new collocation matrix
         if need_resize:
             fmapnii = nb.load(self.inputs.in_data)
-            fmapnii = nb.as_closest_canonical(fmapnii)
+            # fmapnii = nb.as_closest_canonical(fmapnii)
             data = fmapnii.get_fdata(dtype="float32") - center
             if masknii is not None:
                 masknii = nb.load(self.inputs.in_mask)
-                masknii = nb.as_closest_canonical(masknii)
+                # masknii = nb.as_closest_canonical(masknii)
                 mask = np.asanyarray(masknii.dataobj) > 1e-4
             else:
                 mask = np.ones_like(fmapnii.dataobj, dtype=bool)
@@ -268,14 +268,20 @@ def _run_interface(self, runtime):
         # Store interpolated field
         hdr = fmapnii.header.copy()
         hdr.set_data_dtype("float32")
-        fmapnii.__class__(interp_data, fmapnii.affine, hdr).to_filename(out_name)
+        outnii = fmapnii.__class__(interp_data, fmapnii.affine, hdr)
+        outnii.header["cal_max"] = np.abs(outnii.dataobj).max()
+        outnii.header["cal_min"] = -outnii.header["cal_max"]
+        outnii.to_filename(out_name)
         self._results["out_field"] = out_name
 
         # Write out fitting-error map
         self._results["out_error"] = out_name.replace("_field.", "_error.")
-        fmapnii.__class__(
-            data * mask - interp_data + model.intercept_, fmapnii.affine, fmapnii.header
-        ).to_filename(self._results["out_error"])
+        errornii = fmapnii.__class__(
+            (data - interp_data) * mask, fmapnii.affine, fmapnii.header
+        )
+        errornii.header["cal_min"] = 0
+        errornii.header["cal_max"] = np.max(errornii.dataobj)
+        errornii.to_filename(self._results["out_error"])
 
         if not self.inputs.extrapolate:
             return runtime

sdcflows/interfaces/tests/test_bspline.py

@@ -40,41 +40,64 @@
 @pytest.mark.parametrize("testnum", range(100))
 def test_bsplines(tmp_path, testnum):
     """Test idempotency of B-Splines interpolation + approximation."""
-    targetshape = (10, 12, 9)
+    targetshape = (50, 50, 30)
 
     # Generate an oblique affine matrix for the target - it will be a common case.
     targetaff = nb.affines.from_matvec(
-        nb.eulerangles.euler2mat(x=0.9, y=0.001, z=0.001) @ np.diag((2, 3, 4)),
+        nb.eulerangles.euler2mat(x=-0.9, y=0.001, z=0.001) @ np.diag((2, 2, 2.4)),
     )
 
     # Intendedly mis-centered (exercise we may not have volume-centered NIfTIs)
     targetaff[:3, 3] = nb.affines.apply_affine(
         targetaff, 0.5 * (np.array(targetshape) - 3)
     )
 
+    mask = np.zeros(targetshape)
+    mask[10:-10, 10:-10, 6:-6] = 1
     # Generate some target grid
-    targetnii = nb.Nifti1Image(np.ones(targetshape), targetaff, None)
-    targetnii.to_filename(tmp_path / "target.nii.gz")
+    targetnii = nb.Nifti1Image(mask, targetaff, None)
+    targetnii.header.set_qform(targetaff, code=1)
+    targetnii.header.set_sform(targetaff, code=1)
+    targetnii.to_filename(tmp_path / "mask.nii.gz")
 
     # Generate random coefficients
-    gridnii = bspline_grid(targetnii, control_zooms_mm=(4, 6, 8))
-    coeff = (rng.random(size=gridnii.shape) - 0.5) * 500
+    gridnii = bspline_grid(targetnii, control_zooms_mm=(40, 40, 16))
+    coeff = (rng.standard_normal(size=gridnii.shape)) * 100
     coeffnii = nb.Nifti1Image(coeff.astype("float32"), gridnii.affine, gridnii.header)
+    coeffnii.header["cal_max"] = np.abs(coeff).max()
+    coeffnii.header["cal_min"] = -coeffnii.header["cal_max"]
+    coeffnii.header.set_qform(gridnii.affine, code=1)
+    coeffnii.header.set_sform(gridnii.affine, code=1)
     coeffnii.to_filename(tmp_path / "coeffs.nii.gz")
 
     os.chdir(tmp_path)
     # Check that we can interpolate the coefficients on a target
     test1 = ApplyCoeffsField(
-        in_data=str(tmp_path / "target.nii.gz"),
+        in_data=str(tmp_path / "mask.nii.gz"),
         in_coeff=str(tmp_path / "coeffs.nii.gz"),
         pe_dir="j-",
         ro_time=1.0,
     ).run()
 
+    fieldnii = nb.load(test1.outputs.out_field)
+    fielddata = fieldnii.get_fdata()
+    fielddata -= np.median(fielddata)
+    fielddata = 200 * fielddata / np.abs(fielddata).max()
+
+    fieldnii.header["cal_max"] = np.abs(fielddata).max()
+    fieldnii.header["cal_min"] = -fieldnii.header["cal_max"]
+    fieldnii.header.set_qform(targetaff, code=1)
+    fieldnii.header.set_sform(targetaff, code=1)
+
+    nb.Nifti1Image(fielddata, targetaff, fieldnii.header).to_filename(
+        tmp_path / "testfield.nii.gz",
+    )
+
     # Approximate the interpolated target
     test2 = BSplineApprox(
-        in_data=test1.outputs.out_field,
-        bs_spacing=[(4, 6, 8)],
+        in_data=str(tmp_path / "testfield.nii.gz"),
+        # in_mask=str(tmp_path / "mask.nii.gz"),
+        bs_spacing=[(40, 40, 16)],
         zooms_min=0,
         recenter=False,
         ridge_alpha=1e-4,