Fixed precision for NumPy output.

src/tirl/scripts/mnd/find_sites.py

@@ -102,6 +102,8 @@ OPTIONS_NONLINEAR = {
     "visualise": False
 }
 
+np.set_printoptions(precision=4)
+
 
 # IMPLEMENTATION
 
@@ -350,7 +352,7 @@ def _worker_func(job, outputdir, p):
         i.domain.chain[-1] = warp.resample(warp.domain.resize(0.25))
 
     # Identify insertion sites
-    sites = identify_sites(i, s, p)
+    sites = identify_sites(i, s, p, index)
 
     # Save voxel-to-voxel transformation chains
     # (they are saved because they are hard to pickle at the moment)
@@ -475,7 +477,7 @@ def diffreg2d(fixed, moving, logger=None, **options):
         moving.resample(1, copy=False)
 
 
-def identify_sites(intact, sampled, p):
+def identify_sites(intact, sampled, p, index):
     """
     Finds potential insertion sites given a registered pair of photographs
     taken before and after the tissue block sampling. At any one stage, there
@@ -509,18 +511,44 @@ def identify_sites(intact, sampled, p):
         plt.imshow(binary)
         plt.show()
 
-    # Streak filter
+    # Generate distance map
     difference = median(difference, np.ones((10, 10)))
     distmap = edt(difference)
-    binary = np.ones_like(distmap, dtype=np.int64)
-    binary[distmap < p.width / res] = 0
 
-    if p.debug:
-        plt.imshow(binary)
-        plt.show()
+    # How many blocks are expected? (0: data driven)
+    if p.n:
+        n = p.n[index] if len(p.n) > 1 else p.n[0]
+    else:
+        n = 0
+
+    # Use k-Means clustering if a certain number of blocks are expected
+    if n > 0:
+        # Streak filtering
+        distmap[distmap < p.width / res] = 0  # streak filter on distance map
+        if p.debug:
+            plt.imshow(distmap)
+            plt.show()
+
+        # Fit a fixed number of region centroids
+        from sklearn.cluster import KMeans
+        coords = np.stack(np.meshgrid(
+            *tuple(np.arange(dim) for dim in distmap.shape),
+            indexing="ij"), axis=-1).reshape((-1, 2))
+        km = KMeans(n_clusters=n)
+        km.fit(X=coords[distmap.ravel() > 0, :])
+        sites = km.cluster_centers_.tolist()
+
+    # Get sites from binary segmentation
+    else:
+        # Streak filtering
+        binary = np.ones_like(distmap, dtype=np.int64)
+        binary[distmap < p.width / res] = 0
+        if p.debug:
+            plt.imshow(binary)
+            plt.show()
+        # Extract region centroids
+        sites = [region.centroid for region in regionprops(label(binary))]
 
-    # Extract region centroids
-    sites = [region.centroid for region in regionprops(label(binary))]
     if sites:
         return np.atleast_2d(sites)
     else:
@@ -560,7 +588,10 @@ def find_sites(args):
                     f"{' '.join(sys.argv)}")
 
     # Sort images based on foreground area
-    sorted_images = sort_images(*args.images)
+    if p.nosort:
+        sorted_images = args.images
+    else:
+        sorted_images = sort_images(*args.images)
 
     # Perform parallel pairwise rigid registrations to find sites
     # and conclude run
@@ -571,7 +602,7 @@ def find_sites(args):
         raise
 
     if sites is not None:
-        np.savetxt(args.sites, sites)
+        np.savetxt(args.sites, sites, fmt='%1.3f')
         logger.fatal(f"Insertion site search complete. "
                      f"Found {sites.shape[0]} site(s).")
     else:
@@ -613,6 +644,14 @@ def create_cli():
                         required=False,
                         help="Only detect blocks that are larger than x mm2. "
                              "Default: 100 mm2")
+    parser.add_argument("-n", metavar="N", type=int, default=0, nargs="+",
+                        required=False,
+                        help="Look for N0, N1, ... number of sites between "
+                             "consecutive image pairs. Default: 0 (not fixed)")
+    parser.add_argument("--nosort", default=False, action="store_true",
+                        required=False,
+                        help="Do not sort images. Handle them in the order as "
+                             "they are specified.")
 
     parser.add_argument("--threshold", type=float, default=0.1, metavar="x",
                         required=False, help="Fractional lower intensity "