detectSpheres.py 45 KB
 Edward Andò committed Oct 16, 2020 1 2 3 4 ''' Implementation of tomopack by Stéphane Roux '''  Benjy Marks committed Nov 30, 2020 5 import os  Edward Andò committed Oct 16, 2020 6 import numpy  Benjy Marks committed Oct 28, 2021 7 import projectSphere  Edward Andò committed Nov 02, 2020 8 import scipy.ndimage  Benjy Marks committed Jan 21, 2021 9 from scipy.spatial import distance  Edward Andò committed Oct 16, 2020 10 11 12 13  import matplotlib.pyplot as plt import tifffile  Benjy Marks committed Nov 30, 2020 14 15 16 # _kernel = numpy.ones((1,3,3))/9. # square kernel _kernel = numpy.array([[[1,2,1],[2,4,2],[1,2,1]]])/16. # gaussian kernel  Edward Andò committed Oct 26, 2020 17 18 19 20 # non-naiive approach def _PI_projector(f): g = f.copy() g[f<0.5] = 0.  Benjy Marks committed Oct 26, 2020 21  # g[g.imag < 1e-14] = numpy.round(g[g.imag < 1e-14])  Edward Andò committed Oct 26, 2020 22  # g.real = numpy.round(g.real)  Benjy Marks committed Oct 26, 2020 23  g = numpy.round(numpy.abs(g))  Edward Andò committed Oct 26, 2020 24  return g  Edward Andò committed Oct 19, 2020 25   Edward Andò committed Nov 02, 2020 26 # Bottom of page 8 in TomoPack.pdf --- just take peaks in a 3x3 pixel area, and weight by all of the mass in that area  Benjy Marks committed Nov 30, 2020 27 # Why not directly look at maxima in the masses, this avoids the  Edward Andò committed Nov 23, 2020 28 29 30 31 32 33 34 35 36 def _filter_maxima(f, debug=False, removeNegatives=False): # Detect unpadded 2D image first: if len(f.shape) == 2: f = f[numpy.newaxis, ...] twoD = True else: twoD = False # This is making our complex result into a scalar result  Edward Andò committed Nov 02, 2020 37  f_abs = numpy.abs(f)  Benjy Marks committed Nov 30, 2020 38   Edward Andò committed Nov 23, 2020 39 40 41  # Remove negatives -- however keeping negatives can help sharpen the result of the convolution if removeNegatives: f_abs[f<0] = 0  Benjy Marks committed Nov 30, 2020 42  masses = scipy.ndimage.convolve(f_abs,_kernel)  Edward Andò committed Nov 23, 2020 43 44 45  peaks = masses == scipy.ndimage.maximum_filter(masses,size=(1,3,3))  Edward Andò committed Nov 02, 2020 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66  if debug: #plt.imshow(numpy.real(f)) #plt.title('real part of original estimation') #plt.colorbar() #plt.show() #plt.imshow(f_abs) #plt.title('positive real component') #plt.colorbar() #plt.show() #plt.imshow(g) #plt.colorbar() #plt.title('Local filtered maxima') #plt.show() #plt.imshow(peaks) #plt.colorbar() #plt.title('Location of local peaks') #plt.show()  Edward Andò committed Nov 23, 2020 67  plt.imshow(masses[0])  Edward Andò committed Nov 02, 2020 68  plt.colorbar()  Edward Andò committed Nov 23, 2020 69  plt.title('Local masses (first slice if 3D)')  Edward Andò committed Nov 02, 2020 70 71  plt.show()  Edward Andò committed Nov 23, 2020 72  plt.imshow(masses[0]>0.5)  Edward Andò committed Nov 02, 2020 73  plt.colorbar()  Edward Andò committed Nov 23, 2020 74  plt.title('Local masses > 0.5 (first slice if 3D)')  Edward Andò committed Nov 02, 2020 75 76 77 78 79 80  plt.show() #plt.imshow(peaks*masses) #plt.colorbar() #plt.title('Final weighted estimation') #plt.show()  Edward Andò committed Nov 23, 2020 81 82 83  if twoD: return peaks[0]*masses[0] else: return peaks *masses  Edward Andò committed Nov 02, 2020 84   Edward Andò committed May 10, 2021 85   Benjy Marks committed Dec 14, 2020 86 def tomopack(radioMM, psiMM, maxIterations=50, l=0.1, epsilon='iterative', kTrustRatio=0.75, GRAPH=False):  Edward Andò committed Oct 19, 2020 87 88  ''' 'tomopack' FFT-based algorithm for sphere identification by Stéphane Roux.  Edward Andò committed Oct 22, 2020 89 90  the FFT of psiMM as a structuring element to pick up projected spheres in radioMM using a parallel projection hypothesis.  Edward Andò committed Oct 19, 2020 91 92 93  Parameters ----------  Edward Andò committed Oct 22, 2020 94 95 96 97 98 99  radioMM : 2D numpy array of floats radiography containing "distance" projections in mm This is typically the result of mu*log(I/I_0) psiMM : 2D numpy array of floats A (synthetic) projection of a single sphere in approx the same settings,  Edward Andò committed Oct 26, 2020 100 101  as radioMM -- this is the structuring element used for FFT recognition. Should be same size as radioMM  Edward Andò committed Oct 22, 2020 102 103 104 105 106 107 108 109 110  maxIterations : int (optional) Number of iterations to run the detection for. Default = 50 l : float (optional) "lambda" parameter which controls relaxation in the iterations. Default = 0.1  Edward Andò committed Nov 02, 2020 111  epsilon : float (optional), or 'iterative'  Edward Andò committed Oct 22, 2020 112  trust cutoff in psi.  Edward Andò committed Nov 02, 2020 113 114  If the radiograph is in mm so is this parameter Default = 'iterative'  Edward Andò committed Oct 22, 2020 115 116 117 118 119 120  GRAPH : bool (optional) VERY noisily make graphs? Returns -------  Edward Andò committed Nov 02, 2020 121  f_x : 2D numpy array  Edward Andò committed Oct 26, 2020 122 123  same size as radioMM, high values where we think particles are. Consider using indicatorFunctionToDetectorPositions()  Edward Andò committed Oct 19, 2020 124 125  '''  Edward Andò committed Oct 26, 2020 126 127 128 129 130  #assert(psi.sum() < radioMM.sum()), "detectSpheres.tomopack(): sum(psi) > sum(radioMM), doesn't look right" assert(len(psiMM.shape) == 2), "detectSpheres.tomopack(): psi and radioMM should be same size" assert(len(psiMM.shape) == len(radioMM.shape)), "detectSpheres.tomopack(): psi and radioMM should be same size" assert(psiMM.shape[0] == radioMM.shape[0]), "detectSpheres.tomopack(): psi and radioMM should be same size" assert(psiMM.shape[1] == radioMM.shape[1]), "detectSpheres.tomopack(): psi and radioMM should be same size"  Edward Andò committed Oct 19, 2020 131 132 133 134 135  radioMM_FFT = numpy.fft.fft2(radioMM) # define the psi function and its FFT -- the structuring element for SR # Start with the projection of a centred sphere  Olga Stamati committed Oct 26, 2021 136 137 138 139 140  #psiMM_FFT = numpy.fft.fft2(psiMM) # 2021-10-26 OS: phase shift of psi so that zero phase angle is positioned at the centre of the detector # avoids the HACK of shifting f_x later on line 187 psiMM_FFT = numpy.fft.fft2(numpy.fft.fftshift(psiMM))  Edward Andò committed Oct 19, 2020 141 142  # This is comparable to Figure 7 in TomoPack  Benjy Marks committed Oct 26, 2020 143 144 145 146 147  # if GRAPH: # plt.imshow(1./numpy.abs(psiMM_FFT), vmin=0, vmax=1, cmap='hot') # plt.title(r"$1/|\psi_{FFT}|$") # plt.colorbar() # plt.show()  Edward Andò committed Oct 19, 2020 148 149  ## naiive approach  Benjy Marks committed Nov 30, 2020 150 151  with numpy.errstate(divide='ignore', invalid='ignore'): f_k_naiive = radioMM_FFT/psiMM_FFT  Edward Andò committed Oct 19, 2020 152 153 154 155  f_x_naiive = numpy.fft.ifft2(f_k_naiive) # Prepare for iterations  Benjy Marks committed Dec 14, 2020 156   Edward Andò committed Oct 19, 2020 157  f_x = numpy.zeros_like(radioMM)  Benjy Marks committed Oct 26, 2020 158 159 160  if epsilon == 'iterative': epsilon = 1e-5 k_trust = numpy.abs(psiMM_FFT) > epsilon  Benjy Marks committed Dec 14, 2020 161  while k_trust.sum()/k_trust.shape[0]/k_trust.shape[1] > kTrustRatio:  Benjy Marks committed Oct 26, 2020 162 163 164 165 166 167 168 169 170  # print(epsilon) epsilon *= 1.5 k_trust = numpy.abs(psiMM_FFT) > epsilon f_x_old = epsilon*numpy.ones_like(radioMM) else: f_x_old = epsilon*numpy.ones_like(radioMM) k_trust = numpy.abs(psiMM_FFT) > epsilon if GRAPH: print("INFO: k_trust maintains {}% of wavenumubers".format(100*k_trust.sum()/k_trust.shape[0]/k_trust.shape[1]))  Edward Andò committed Oct 19, 2020 171 172  count = 0  Edward Andò committed Oct 22, 2020 173  if GRAPH: plt.ion()  Edward Andò committed Oct 19, 2020 174 175 176 177 178  # Start iterations as per Algorithm 1 in TomoPack # Objective: get a good indictor function f_x while numpy.linalg.norm(f_x - f_x_old) > 1e-8: # NOTE: Using a different epsilon here  Edward Andò committed Nov 02, 2020 179  if GRAPH: plt.clf()  Edward Andò committed Oct 19, 2020 180 181 182  f_x_old = f_x.copy() f_k = numpy.fft.fft2(f_x)  Edward Andò committed Oct 22, 2020 183  # if GRAPH: plt.plot(f_k,'k.')  Edward Andò committed Oct 19, 2020 184  f_k[k_trust] = f_k_naiive[k_trust]  Edward Andò committed Oct 22, 2020 185  # if GRAPH: plt.plot(f_k,'r.')  Olga Stamati committed Oct 26, 2021 186 187  f_x = numpy.fft.ifft2(f_k) #f_x = numpy.fft.fftshift(numpy.fft.ifft2(f_k)) # importing HACK from 1D version  Edward Andò committed Oct 26, 2020 188  f_x = f_x + l*(_PI_projector(f_x) - f_x)  Edward Andò committed Oct 19, 2020 189 190  #print(numpy.amax(numpy.abs(f_x - f_x_old)))  Edward Andò committed Oct 22, 2020 191  if GRAPH:  Edward Andò committed Oct 19, 2020 192 193 194  plt.imshow(numpy.abs(f_x), vmin=0)#, vmax=1) plt.colorbar() # plt.plot(numpy.real(f_x),numpy.imag(f_x),'.')  Benjy Marks committed Oct 26, 2020 195  plt.pause(0.00001)  Edward Andò committed Oct 19, 2020 196 197 198  count += 1 if count > maxIterations:  Benjy Marks committed Dec 14, 2020 199 200  # f_x = numpy.zeros_like(f_x) # NOTE: THIS WILL RETURN NO MATCHES WHEN TOMOPACK DOESN'T CONVERGE print('\nKILLED LOOP')  Edward Andò committed Oct 19, 2020 201  break  Benjy Marks committed Oct 22, 2020 202 203  # plt.show() # print('Took ' + str(count) + ' iterations to complete.')  Edward Andò committed Oct 19, 2020 204   Benjy Marks committed Oct 26, 2020 205 206  if GRAPH: #not_k_trust_ind = numpy.where(~k_trust)  Edward Andò committed Oct 19, 2020 207   Benjy Marks committed Oct 26, 2020 208  plt.figure(figsize=[8,6])  Edward Andò committed Oct 19, 2020 209   Benjy Marks committed Oct 26, 2020 210 211 212 213 214  plt.subplot(331) plt.title(r'$p(x)$') plt.imshow(radioMM) # plt.plot(pPredictedIJ[:,0], pPredictedIJ[:,1], 'rx') plt.colorbar()  Edward Andò committed Oct 19, 2020 215   Benjy Marks committed Oct 26, 2020 216 217 218 219  plt.subplot(332) plt.title(r'$\psi(x)$') plt.imshow(psiMM, vmin=radioMM.min(), vmax=radioMM.max()) plt.colorbar()  Edward Andò committed Oct 19, 2020 220   Benjy Marks committed Oct 26, 2020 221 222 223 224 225  ax3 = plt.subplot(333) plt.title(r'True $f(x)$') # HACK ky-flat image plt.imshow(numpy.zeros_like(radioMM)) #ax3.set_aspect('equal', 'datalim')  Edward Andò committed Oct 19, 2020 226   Benjy Marks committed Oct 26, 2020 227 228 229 230  plt.subplot(334) plt.title(r'$|\tilde{p}(k)|$') plt.imshow(numpy.abs(radioMM_FFT))#, vmin=0, vmax=10) plt.colorbar()  Edward Andò committed Oct 19, 2020 231   Benjy Marks committed Oct 26, 2020 232 233 234 235  plt.subplot(335) plt.title(r'$|\tilde{\psi}(k)|$') plt.imshow(numpy.abs(psiMM_FFT), vmin=0, vmax=10) plt.colorbar()  Edward Andò committed Oct 19, 2020 236   Benjy Marks committed Oct 26, 2020 237 238 239 240 241 242 243 244 245 246 247  one_on_psi_plot = 1./numpy.abs(psiMM_FFT) one_on_psi_plot[~k_trust] = numpy.nan plt.subplot(336) plt.title(r'$1/|\psi(k)|$') plt.imshow(one_on_psi_plot) # plt.semilogy(numpy.arange(len(psi_FFT))[k_trust],1./numpy.abs(psi_FFT)[k_trust],'g.') #plt.semilogy(numpy.arange(len(psi_FFT))[~k_trust],1./numpy.abs(psi_FFT)[~k_trust],'r.') #plt.semilogy([0,len(psi_FFT)],[1./epsilon,1./epsilon],'k--') # plt.plot(not_k_trust_ind[0], not_k_trust_ind[1], 'r.') plt.colorbar() # plt.ylim(0,1./epsilon)  Edward Andò committed Oct 19, 2020 248   Edward Andò committed Oct 26, 2020 249   Benjy Marks committed Oct 26, 2020 250 251 252 253  plt.subplot(337) plt.title(r'Estimated $|f(k)|$') plt.imshow(numpy.abs(f_k)) plt.colorbar()  Edward Andò committed Oct 26, 2020 254   Benjy Marks committed Oct 26, 2020 255 256 257 258 259  plt.subplot(338) plt.title(r'Estimated $f(x)$') plt.imshow(numpy.real(f_x)) # plt.plot(pPredictedIJ[:,0], pPredictedIJ[:,1], 'o', mec='r', mfc='None', mew=0.1) plt.colorbar()  Edward Andò committed Oct 26, 2020 260   Benjy Marks committed Oct 26, 2020 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288  #plt.subplot(339) #plt.title(r'Residual $p(x) - p(f(x))$') #p_f_x = numpy.zeros_like(radioMM, dtype=float) ### Convert to XYZ in space and mm #positionsXYZmm = numpy.zeros([pPredictedIJ.shape[0], 3]) ## Y #positionsXYZmm[:,1] = -1*(pPredictedIJ[:,0] - radioMM.shape[0]/2.0)*pixelSizeMM ## Z #positionsXYZmm[:,2] = -1*(pPredictedIJ[:,1] - radioMM.shape[1]/2.0)*pixelSizeMM #radii = numpy.ones([pPredictedIJ.shape[0]]) #radii *= radiusMM #print("\n\n\n", pPredictedIJ, "\n\n\n") #print("\n\n\n", positionsXYZmm, "\n\n\n") #p_f_x = radioSphere.projectSphere.projectSphereMM(positionsXYZmm, #radii, #sourceDetectorDistMM=numpy.inf, #pixelSizeMM=pixelSizeMM, #detectorResolution=detectorResolution) #residual = p_f_x #residual = p_f_x - radioMM #vmax = numpy.abs(residual).max() #plt.imshow(residual, vmin=-vmax, vmax=vmax, cmap='coolwarm') #plt.colorbar()  Edward Andò committed Oct 26, 2020 289 290 291   Benjy Marks committed Oct 26, 2020 292 293 294  plt.subplots_adjust(hspace=0.5) plt.show() #plt.savefig('tt3D.png',dpi=200)  Edward Andò committed Oct 26, 2020 295 296  #return pPredictedIJ  Benjy Marks committed Oct 26, 2020 297  return numpy.real(f_x)  Edward Andò committed Oct 26, 2020 298   Edward Andò committed Nov 23, 2020 299   Edward Andò committed Oct 26, 2020 300 301 302 303 def indicatorFunctionToDetectorPositions(f_x, debug=False, particlesPresenceThreshold=0.9): ''' Takes an approximation of the indicator function f_x from tomopack and returns positions on the detector  Edward Andò committed Nov 02, 2020 304 305 306 307  Parameters ---------- f_x : 2D numpy array of floats The value of the indicator function for each pixel on the detector  Edward Andò committed Nov 02, 2020 308   Edward Andò committed Nov 02, 2020 309 310 311  debug : bool, optional Show debug graphs (especially in the maximum filtering) Default = False  Edward Andò committed Nov 02, 2020 312   Edward Andò committed Nov 02, 2020 313 314 315 316 317  particlesPresence : float, optional Threshold for accepting a particle Default = 0.9 ''' assert(len(f_x.shape) == 2), "detectSpheres.indicatorFunctionToDetectorPositions(): need 2D array"  Edward Andò committed Nov 02, 2020 318   Edward Andò committed Nov 02, 2020 319  f_x = _filter_maxima(f_x, debug=debug)  Edward Andò committed Oct 26, 2020 320 321 322 323 324  particlesPresence = numpy.round(numpy.real(_PI_projector(f_x))) pPredictedIJ = [] for pos in numpy.array(numpy.where(particlesPresence > particlesPresenceThreshold)).T: val = particlesPresence[pos[0], pos[1]]  Benjy Marks committed Oct 26, 2020 325  # print(f_x[pos[0],pos[1]])  Edward Andò committed Oct 26, 2020 326 327 328  for i in range(int(val)): pPredictedIJ.append([pos[1], pos[0]]) pPredictedIJ = numpy.array(pPredictedIJ)  Edward Andò committed Oct 19, 2020 329 330  return pPredictedIJ  Edward Andò committed Nov 23, 2020 331 332   Edward Andò committed May 10, 2021 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 def psiSeriesScanTo3DPositions(radio, # not (necessarily) MM psiXseries, # Obviously in the same units as radio above, please radiusMM, # can be removed? CORxPositions = None, massThreshold=0.12, scanPersistenceThresholdRadii=None, 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): # it is our objective to fill in fx series fXseries = numpy.zeros_like(psiXseries) if CORxPositions is None: print("xPositions is not passed, just putting 1, 2, 3...") CORxPositions = numpy.arange(psiXseries.shape[0]) for posN, CORxPos in enumerate(CORxPositions): ### "Structuring Element" print("\t{}/{} CORxPos = {:0.2f} mm".format(posN+1, len(CORxPositions), CORxPos), end='\r') fXseries[posN] = radioSphere.detectSpheres.tomopack(radio, psiXseries[posN], GRAPH=0, maxIterations=maxIterations, l=l, kTrustRatio=kTrustRatio)  Edward Andò committed May 18, 2021 363 364  tifffile.imsave(cacheFile, fXseries.astype('float')) print(f"saved {cacheFile}")  Edward Andò committed May 10, 2021 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409  #loadedCache = False #if useCache: #cachePsiFile = cacheFile[:-4] + '_psi.tif' #if os.path.isfile(cacheFile) and os.path.isfile(cachePsiFile): #print("Loading previous indicator functions... ", end="") #fXseries = tifffile.imread(cacheFile) #psiXseries = tifffile.imread(cachePsiFile) #if ( fXseries.shape[0] == CORxNumber ) and ( fXseries.shape[1] == radioMM.shape[0] ) and ( fXseries.shape[2] == radioMM.shape[1] ): #print("done.") #loadedCache = True #else: #print("cached file had wrong dimensions. Generating new cache file.") #else: #print('No cached indicator functions found. Generating them now to cache.') #if not loadedCache: #fXseries = numpy.zeros((len(CORxPositions), radioMM.shape[0], radioMM.shape[1])) #psiXseries = numpy.zeros_like(fXseries) #psiRefMM = radioSphere.projectSphere.projectSphereMM(numpy.array([[(CORxMax+CORxMin)/2., 0., 0.]]), #numpy.array([radiusMM]), #detectorResolution=radioMM.shape, #pixelSizeMM=pixelSizeMM, #sourceDetectorDistMM=sourceDetectorDistMM, #blur=blur) #for posN, CORxPos in enumerate(CORxPositions): #### "Structuring Element" #print("\t{}/{} CORxPos = {:0.2f}mm".format(posN+1, len(CORxPositions), CORxPos), end='\r') #psiMM = radioSphere.projectSphere.projectSphereMM(numpy.array([[CORxPos, 0., 0.]]), #numpy.array([radiusMM]), #detectorResolution=radioMM.shape, #pixelSizeMM=pixelSizeMM, #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) #if useCache and not loadedCache: #print("Saving indicator functions for next time... ", end="") #tifffile.imsave(cacheFile, fXseries.astype(' massThreshold  Edward Andò committed May 10, 2021 424 425 426 427  zoomLevel = sourceDetectorDistMM/((CORxPositions[0] + CORxPositions[-1])/2) CORxDelta = numpy.abs(CORxPositions[0]-CORxPositions[1]) # Look in a volume of +/- half a radius in all directions for the highest value (+/- 1 radius keeps overlapping and causing issues, half a radius doesn't overlap particles, but still contains one clean peak)  Edward Andò committed Jun 08, 2021 428 429  fXseriesMaximumFiltered = scipy.ndimage.maximum_filter(fXseries, size=(3*numpy.int(numpy.floor(radiusMM/CORxDelta)),  Edward Andò committed May 10, 2021 430  numpy.int(numpy.floor(radiusMM/pixelSizeMM*zoomLevel)),  Edward Andò committed Jun 08, 2021 431 432 433 434 435  numpy.int(numpy.floor(radiusMM/pixelSizeMM*zoomLevel))) ) allPeaks = fXseries == fXseriesMaximumFiltered masses = allPeaks*fXseries  Edward Andò committed May 10, 2021 436 437 438 439  if verbose: tifffile.imsave(cacheFile[:-4] + '_masses.tif', masses.astype(' massThreshold peaksCORxPOSnJI = numpy.argwhere(filteredPeaks) if verbose: tifffile.imsave(cacheFile[:-4] + '_filteredPeaks.tif', filteredPeaks.astype(' massThreshold structure = numpy.ones([5,3,3]) binaryPeaks = scipy.ndimage.binary_closing(binaryPeaks,structure=structure) # remove any gaps  Benjy Marks committed Jan 18, 2021 587  if verbose: tifffile.imsave('closed_peaks.tif', binaryPeaks.astype(' massThreshold  Benjy Marks committed Nov 30, 2020 603   Benjy Marks committed Jan 14, 2021 604  get_peaks_from_persistence = False  Benjy Marks committed Dec 24, 2020 605  if get_peaks_from_persistence:  Edward Andò committed Dec 08, 2020 606   Benjy Marks committed Dec 24, 2020 607 608  labelledPeaks = scipy.ndimage.label(binaryPeaks)[0] labelledPeaksBB = scipy.ndimage.find_objects(labelledPeaks)  Edward Andò committed Dec 08, 2020 609   Benjy Marks committed Dec 24, 2020 610 611 612  print(f"Found {labelledPeaks.max()} peaks across the indicator function series.") # Compute centres of mass and persistence in the scanning direction peaksCORxPOSnJI = numpy.zeros((labelledPeaks.max(), 4))  Edward Andò committed Nov 23, 2020 613   Benjy Marks committed Dec 24, 2020 614 615 616 617  try: import spam.label #print("Doing it with spam") BB = spam.label.boundingBoxes(labelledPeaks)  Edward Andò committed Dec 08, 2020 618 619  # save x-y-z positions where x is the slice number in the scan direction and y, z pixel positions on detector  Benjy Marks committed Dec 24, 2020 620  peaksCORxPOSnJI[:,0:3] = spam.label.centresOfMass(labelledPeaks, boundingBoxes=BB)[1:]  Edward Andò committed Dec 08, 2020 621 622  # also add in how many slices this is alive for  Benjy Marks committed Dec 24, 2020 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653  peaksCORxPOSnJI[:,-1] = BB[1:,1] - BB[1:,0] + 1 #print("Done it with spam") except: print("Not doing peak label analysis with spam (not found) this will be very slow") for label in range(1,labelledPeaks.max()+1): # save x-y-z positions where x is the slice number in the scan direction and y, z pixel positions on detector peaksCORxPOSnJI[label-1][0:3] = scipy.ndimage.center_of_mass(labelledPeaks==label) # also add in how many slices this is alive for peaksCORxPOSnJI[label-1][-1] = labelledPeaksBB[label-1][0].stop - labelledPeaksBB[label-1][0].start + 1 print("...it was slow") if verbose: print("tomopackDivergentScanTo3DPositions(): object persistence in scan direction:\n", peaksCORxPOSnJI[:, -1]) # Apply persistence threshold: persistenceMask = numpy.zeros(labelledPeaks.max(), dtype=bool) CORxSpacing = CORxPositions[1] - CORxPositions[0] if scanPersistenceThresholdRadii is not None: scanPersistenceThreshold = int(numpy.floor((scanPersistenceThresholdRadii*radiusMM)/CORxSpacing)) print(f"Calculated scan persistence threshold for {scanPersistenceThresholdRadii}\% particle radii: {scanPersistenceThreshold}") if scanFixedNumber: # print('Scanning for a fixed number of particles is not implemented!') sortedPeakIndices = numpy.argsort(peaksCORxPOSnJI[:, -1])[::-1] for i in range(scanFixedNumber): persistenceMask[sortedPeakIndices[i]] = True else: # use threshold value persistenceMask[peaksCORxPOSnJI[:, -1] > scanPersistenceThreshold] = True # Apply mask and remove 4th column: peaksCORxPOSnJI = peaksCORxPOSnJI[persistenceMask, 0:3] else:  Benjy Marks committed Jan 14, 2021 654 655 656 657 658 659 660 661 662  zoomLevel = sourceDetectorDistMM/((CORxMin + CORxMax)/2) # Look in a volume of +/- half a radius in all directions for the highest value (+/- 1 radius keeps overlapping and causing issues, half a radius doesn't overlap particles, but still contains one clean peak) fXconvolvedMaximumFiltered = scipy.ndimage.maximum_filter(fXconvolvedSeries, size=(numpy.int(numpy.floor(radiusMM/CORxDelta)), numpy.int(numpy.floor(radiusMM/pixelSizeMM*zoomLevel)), numpy.int(numpy.floor(radiusMM/pixelSizeMM*zoomLevel)) ) ) allPeaks = fXconvolvedSeries == fXconvolvedMaximumFiltered  Benjy Marks committed Dec 24, 2020 663  masses = allPeaks*fXconvolvedSeries  Benjy Marks committed Jan 14, 2021 664  if verbose:  Benjy Marks committed Jan 18, 2021 665 666  tifffile.imsave(cacheFile[:-4] + '_masses.tif', masses.astype(' massThreshold peaksCORxPOSnJI = numpy.argwhere(filteredPeaks)  Benjy Marks committed Jan 14, 2021 680 681  if verbose: tifffile.imsave(cacheFile[:-4] + '_filteredPeaks.tif', filteredPeaks.astype(' 0.9: positionsXYZmm, residual = removeParticle(positionsXYZmm,residual,radioMM,radiiMM,pixelSizeDetectorMM,zoomLevel,sourceObjectDistMM,verbose,GRAPH) removed += 1 while residual.min() < -0.9:  915  positionsXYZmm, residual = addParticle(positionsXYZmm,residual,radioMM,radiiMM,pixelSizeDetectorMM,zoomLevel,sourceObjectDistMM,CORxMin, CORxMax,CORxNumber,verbose,GRAPH)  Benjy Marks committed Jan 14, 2021 916 917 918 919 920 921  added +=1 # now we clean up the cleaning up to make sure we have the right number of particles while added > removed: positionsXYZmm, residual = removeParticle(positionsXYZmm,residual,radioMM,radiiMM,pixelSizeDetectorMM,zoomLevel,sourceObjectDistMM,verbose,GRAPH) removed += 1 while removed > added:  922  positionsXYZmm, residual = addParticle(positionsXYZmm,residual,radioMM,radiiMM,pixelSizeDetectorMM,zoomLevel,sourceObjectDistMM,CORxMin, CORxMax,CORxNumber,verbose,GRAPH)  Benjy Marks committed Jan 14, 2021 923  added +=1  Benjy Marks committed Jan 21, 2021 924 925  print(f'Removed {removed} particles and added {added}.')  Benjy Marks committed Dec 15, 2020 926  return positionsXYZmm  Benjy Marks committed Jan 21, 2021 927   Edward Andò committed May 10, 2021 928   Benjy Marks committed Jan 21, 2021 929 930 def calculateErrors(posXYZa,posXYZb,radiiMM,verbose=True): distances = distance.cdist(posXYZa,posXYZb)  Benjy Marks committed Jan 21, 2021 931 932  distance_threshold = radiiMM[0]  Benjy Marks committed Jan 21, 2021 933 934  # an empty column in the distance matrix is a lost particle number_lost = numpy.sum((numpy.sum(distances <= distance_threshold, axis=0) == 0))  Edward Andò committed Feb 01, 2021 935  if number_lost > 0: print(f"detectSpheres.calculateErrors(): number_lost = {number_lost}")  Benjy Marks committed Jan 21, 2021 936   Benjy Marks committed Jan 21, 2021 937 938 939  # get errors min_errors = numpy.min(distances, axis=0) # best match for each column min_valid_errors = min_errors[min_errors < distance_threshold] # just where we found a particle  Benjy Marks committed Jan 21, 2021 940   Benjy Marks committed Jan 21, 2021 941 942  err_mean = numpy.mean(min_valid_errors) err_std = numpy.std( min_valid_errors)  Benjy Marks committed Jan 21, 2021 943 944  return err_mean, err_std, number_lost