From 7facb16b585026e291d22bbd14a4c7c13817a9ac Mon Sep 17 00:00:00 2001
From: seanmcculloch <sean.mcculloch@alleninstitute.org>
Date: Fri, 6 Dec 2024 14:20:03 -0800
Subject: [PATCH 1/2] add nan_to_num before cast

---
 bigstream/transform.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/bigstream/transform.py b/bigstream/transform.py
index 89461ef..1100296 100644
--- a/bigstream/transform.py
+++ b/bigstream/transform.py
@@ -125,7 +125,7 @@ def apply_transform(
 
     # execute, return as numpy array
     resampled = resampler.Execute(mov)
-    return sitk.GetArrayViewFromImage(resampled).astype(dtype)
+    return np.nan_to_num(sitk.GetArrayViewFromImage(resampled),nan=0.0,posinf=0.0,neginf=0.0).astype(dtype)
 
 
 def apply_transform_to_coordinates(

From d05e45872f6c83fe8485ac9d348db3ef5d95cb8f Mon Sep 17 00:00:00 2001
From: seanmcculloch <sean.mcculloch@alleninstitute.org>
Date: Tue, 10 Dec 2024 16:38:14 -0800
Subject: [PATCH 2/2] chunked interest points

---
 bigstream/features.py | 69 +++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 66 insertions(+), 3 deletions(-)

diff --git a/bigstream/features.py b/bigstream/features.py
index f63c88a..dda3973 100644
--- a/bigstream/features.py
+++ b/bigstream/features.py
@@ -1,8 +1,72 @@
 import numpy as np
 from fishspot.filter import white_tophat, apply_foreground_mask
-from fishspot.detect import detect_spots_log
+
 from scipy.stats.mstats import winsorize
 
+import numpy as np
+from skimage.feature import blob_log
+
+
+# TODO: potential major improvement - after finding coordinates
+#       with LoG filter, template match the PSF to the region around
+#       each detected point (Fourier Phase Correlation maybe?).
+#       upsample the PSF and data to achieve subvoxel accuracy.
+
+
+def detect_spots_log(
+    image,
+    min_radius,
+    max_radius,
+    num_sigma=5,
+    **kwargs,
+):
+    """
+    """
+
+    # ensure iterable radii
+    if not isinstance(min_radius, (tuple, list, np.ndarray)):
+        min_radius = (min_radius,)*image.ndim
+    if not isinstance(max_radius, (tuple, list, np.ndarray)):
+        max_radius = (max_radius,)*image.ndim
+
+    # set given arguments
+    kwargs['min_sigma'] = np.array(min_radius) / np.sqrt(image.ndim)
+    kwargs['max_sigma'] = np.array(max_radius) / np.sqrt(image.ndim)
+    kwargs['num_sigma'] = num_sigma
+
+    # set additional defaults
+    if 'threshold' not in kwargs or kwargs['threshold'] is None:
+        kwargs['threshold'] = None
+        kwargs['threshold_rel'] = 0.1
+
+    # run
+    #return blob_log(image, **kwargs)
+    return chunked_blob_log(image, **kwargs)
+
+
+def chunked_blob_log(image, **kwargs): #sigma_list, chunk_size=(128, 128, 128), overlap=(32, 32, 32)):
+    chunk_size=(256, 256, 256)
+    overlap=(64, 64, 64)
+    z_chunks, y_chunks, x_chunks = [
+        range(0, dim, chunk_size[i] - overlap[i])
+        for i, dim in enumerate(image.shape)
+    ]
+    blobs = []
+    
+    for z_start in z_chunks:
+        for y_start in y_chunks:
+            for x_start in x_chunks:
+                z_end = min(z_start + chunk_size[0], image.shape[0])
+                y_end = min(y_start + chunk_size[1], image.shape[1])
+                x_end = min(x_start + chunk_size[2], image.shape[2])
+                
+                chunk = image[z_start:z_end, y_start:y_end, x_start:x_end]
+                chunk_blobs = blob_log(chunk, **kwargs)
+                chunk_blobs[:, :3] += [z_start, y_start, x_start]  # Adjust positions
+                blobs.extend(chunk_blobs)
+    
+    return np.array(blobs)
+
 
 def blob_detection(
     image,
@@ -25,8 +89,7 @@ def blob_detection(
     max_blob_radius : scalar float
         The largest size blob you want to find in voxel units
     **kwargs : any additional kwargs
-        Passed to fishspot.detect_spots_log
-
+        Passed to fishspot.detect_spots_log∂
     Returns
     -------
     blob_coordinates_and_intensities : nd-array Nx4