Changed mask behaviour: excluding from smoothing.

src/tirl/costs/mind.py

@@ -194,7 +194,7 @@ class CostMIND(CostMSD):
 
         vc = img.domain.get_voxel_coordinates()
         if img.mask is not None:
-            masked = np.flatnonzero(np.isclose(img.mask, 0, atol=1e-3))
+            masked = np.flatnonzero(np.isclose(img.mask, 0, atol=1e-2))
         else:
             masked = []
 
@@ -208,6 +208,11 @@ class CostMIND(CostMSD):
                 if vi:
                     img_shifted.data[...] = np.roll(img_shifted.data, vi, axis)
 
+            if img.mask is not None:
+                ssdmask = img.mask.copy()
+            else:
+                ssdmask = np.ones(img.vshape, dtype=np.float32)
+
             # Calculate voxelwise squared differences over the entire image
             ssd = (img_shifted.data - img.data) ** 2
             del img_shifted
@@ -222,23 +227,39 @@ class CostMIND(CostMSD):
             else:
                 bindices = np.unique(np.concatenate(bindices))
             ssd.flat[bindices] = 0
+            ssdmask.flat[bindices] = 0
 
             # Treat mask edges the same as image borders
+            mindices = []
             if self.metaparameters.get("ignore_masked_edges", False):
                 tc = np.delete(test_coordinates, bindices, axis=0)
                 vcc = np.delete(vc, bindices, axis=0)
+                masked = np.asarray(masked)
                 mindices = np.ravel_multi_index(tc.T, img.vshape)
                 mindices = np.ravel_multi_index(
-                    vcc[np.in1d(mindices, masked)].T, img.vshape)
+                    vcc[np.in1d(mindices, masked), :].T, img.vshape)
                 ssd.flat[mindices] = 0
+                ssdmask.flat[mindices] = 0
+                ssd.flat[masked] = 0
+                ssdmask.flat[masked] = 0
 
             # Convert distances to patch-wise distance
             # This step ensures robust estimation of the distance
             # from the neighbour. Do not smooth across tensor values.
             sigma_nd = (sigma,) * ndim + (0,) * img.tdim
-            out.tensors[i] = gaussian_filter(
-                ssd, sigma=sigma_nd, truncate=truncate)
-            del ssd
+            if self.metaparameters.get("ignore_masked_regions", False):
+                tmp = \
+                    gaussian_filter(ssd, sigma=sigma_nd, truncate=truncate) / \
+                    gaussian_filter(ssdmask, sigma=sigma_nd, truncate=truncate)
+            else:
+                tmp = gaussian_filter(ssd, sigma=sigma_nd, truncate=truncate)
+            tmp[~np.isfinite(tmp)] = 0
+            tmp.flat[bindices] = 0
+            if self.metaparameters.get("ignore_masked_edges", False):
+                tmp.flat[mindices] = 0
+            out.tensors[i] = tmp
+
+            del ssd, tmp
 
         # Calculate MIND descriptor from patch distances
         # This is essentially making MIND vectors universally comparable