fourier convolution of fx series with psi series SE in divergent beam

tools/detectSpheres.py

 '''
-Implementation of tomopack by Stéphane Roux
+Tool for sphere detection in parallel and divergent beam
 '''
 
 import os
 import numpy
 import radioSphere.projectSphere
 import scipy.ndimage
+import scipy.signal
 from scipy.spatial import distance
 
 import matplotlib.pyplot as plt
@@ -437,7 +438,6 @@ def tomopackDivergentScanTo3DPositions(radioMM, radiusMM,
                                        scanFixedNumber=None,
                                        scanPersistenceThreshold=7, maxIterations=50,
                                        sourceDetectorDistMM=100, pixelSizeMM=0.1, l=0.2, kTrustRatio=0.7, useCache=True,
-                                       numCores=1,
                                        blur=0.0,
                                        cacheFile='fXseries.tif',
                                        verbose=False):
@@ -511,8 +511,8 @@ def tomopackDivergentScanTo3DPositions(radioMM, radiusMM,
                                                               sourceDetectorDistMM=sourceDetectorDistMM,
                                                               blur=blur)
 
-            fXseries[posN] = radioSphere.detectSpheres.tomopack(radioMM, psiMM, GRAPH=0, maxIterations=maxIterations, l=l, kTrustRatio=kTrustRatio)
-            psiXseries[posN] = radioSphere.detectSpheres.tomopack(psiRefMM, psiMM, GRAPH=0, maxIterations=maxIterations, l=l, kTrustRatio=kTrustRatio)
+            fXseries[posN] = radioSphere.detectSpheres.tomopack(radioMM, psiMM, maxIterations=maxIterations, l=l, kTrustRatio=kTrustRatio)
+            psiXseries[posN] = radioSphere.detectSpheres.tomopack(psiRefMM, psiMM, maxIterations=maxIterations, l=l, kTrustRatio=kTrustRatio)
             psiMMseries[posN] = psiMM
 
     if useCache and not loadedCache:
@@ -540,14 +540,14 @@ def tomopackDivergentScanTo3DPositions(radioMM, radiusMM,
                             (psiXseries.shape[1])//2 - L_yz:(psiXseries.shape[1])//2 + L_yz + 1,
                             (psiXseries.shape[2])//2 - L_yz:(psiXseries.shape[2])//2 + L_yz + 1]
 
-        fXconvolvedSeries = scipy.ndimage.convolve(fXseries,struct/struct.sum())
+        #fXconvolvedSeries = scipy.ndimage.convolve(fXseries,struct/struct.sum())
+        # OS 2021-11-19: Let's do a fourier convolution which should be way faster
+        fXconvolvedSeries = scipy.signal.convolve(fXseries, struct/struct.sum(), mode="same")
         if useCache and not loadedCache:
             tifffile.imsave(f'{cacheFile[:-4]}_struct.tif', struct.astype('<f4'))
             tifffile.imsave(f'{cacheFile[:-4]}_fXconvolvedSeries.tif', fXconvolvedSeries.astype('<f4'))
 
 
-        binaryPeaks = fXconvolvedSeries > massThreshold
-
     get_peaks_from_persistence = False
     if get_peaks_from_persistence: