detectSpheres.py 45 KB
Newer Older
1
2
3
4
'''
Implementation of tomopack by Stéphane Roux
'''

5
import os
6
import numpy
Benjy Marks's avatar
Benjy Marks committed
7
import projectSphere
Edward Andò's avatar
Edward Andò committed
8
import scipy.ndimage
9
from scipy.spatial import distance
10
11
12
13

import matplotlib.pyplot as plt
import tifffile

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ò's avatar
Edward Andò committed
17
18
19
20
# non-naiive approach
def _PI_projector(f):
    g = f.copy()
    g[f<0.5] = 0.
21
    # g[g.imag < 1e-14] = numpy.round(g[g.imag < 1e-14])
Edward Andò's avatar
Edward Andò committed
22
    # g.real = numpy.round(g.real)
23
    g = numpy.round(numpy.abs(g))
Edward Andò's avatar
Edward Andò committed
24
    return g
25

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
27
# Why not directly look at maxima in the masses, this avoids the
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
37
    f_abs = numpy.abs(f)
38

39
40
41
    # Remove negatives -- however keeping negatives can help sharpen the result of the convolution
    if removeNegatives: f_abs[f<0] = 0

42
    masses = scipy.ndimage.convolve(f_abs,_kernel)
43
44
45

    peaks = masses == scipy.ndimage.maximum_filter(masses,size=(1,3,3))

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()

67
        plt.imshow(masses[0])
68
        plt.colorbar()
69
        plt.title('Local masses (first slice if 3D)')
70
71
        plt.show()

72
        plt.imshow(masses[0]>0.5)
73
        plt.colorbar()
74
        plt.title('Local masses > 0.5 (first slice if 3D)')
75
76
77
78
79
80
        plt.show()

        #plt.imshow(peaks*masses)
        #plt.colorbar()
        #plt.title('Final weighted estimation')
        #plt.show()
81
82
83

    if twoD: return peaks[0]*masses[0]
    else:    return peaks   *masses
84

Edward Andò's avatar
Edward Andò committed
85

Benjy Marks's avatar
Benjy Marks committed
86
def tomopack(radioMM, psiMM, maxIterations=50, l=0.1, epsilon='iterative', kTrustRatio=0.75, GRAPH=False):
87
88
    '''
    'tomopack' FFT-based algorithm for sphere identification by Stéphane Roux.
Edward Andò's avatar
Edward Andò committed
89
90
    the FFT of psiMM as a structuring element to pick up projected spheres in radioMM
    using a parallel projection hypothesis.
91
92
93

    Parameters
    ----------
Edward Andò's avatar
Edward Andò committed
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ò's avatar
Edward Andò committed
100
101
            as radioMM -- this is the structuring element used for FFT recognition.
            Should be same size as radioMM
Edward Andò's avatar
Edward Andò committed
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

111
        epsilon : float (optional), or 'iterative'
Edward Andò's avatar
Edward Andò committed
112
            trust cutoff in psi.
113
114
            If the radiograph is in mm so is this parameter
            Default = 'iterative'
Edward Andò's avatar
Edward Andò committed
115
116
117
118
119
120

        GRAPH : bool (optional)
            VERY noisily make graphs?

    Returns
    -------
121
        f_x : 2D numpy array
Edward Andò's avatar
Edward Andò committed
122
123
            same size as radioMM, high values where we think particles are.
            Consider using indicatorFunctionToDetectorPositions()
124
125

    '''
Edward Andò's avatar
Edward Andò committed
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"
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
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))

141
142

    # This is comparable to Figure 7 in TomoPack
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()
148
149

    ## naiive approach
150
151
    with numpy.errstate(divide='ignore', invalid='ignore'):
        f_k_naiive = radioMM_FFT/psiMM_FFT
152
153
154
155
    f_x_naiive = numpy.fft.ifft2(f_k_naiive)


    # Prepare for iterations
Benjy Marks's avatar
Benjy Marks committed
156

157
    f_x = numpy.zeros_like(radioMM)
158
159
160
    if epsilon == 'iterative':
        epsilon = 1e-5
        k_trust = numpy.abs(psiMM_FFT) > epsilon
Benjy Marks's avatar
Benjy Marks committed
161
        while k_trust.sum()/k_trust.shape[0]/k_trust.shape[1] > kTrustRatio:
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]))
171
172

    count = 0
Edward Andò's avatar
Edward Andò committed
173
    if GRAPH: plt.ion()
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
179
        if GRAPH: plt.clf()
180
181
182

        f_x_old = f_x.copy()
        f_k = numpy.fft.fft2(f_x)
Edward Andò's avatar
Edward Andò committed
183
        # if GRAPH: plt.plot(f_k,'k.')
184
        f_k[k_trust] = f_k_naiive[k_trust]
Edward Andò's avatar
Edward Andò committed
185
        # if GRAPH: plt.plot(f_k,'r.')
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ò's avatar
Edward Andò committed
188
        f_x = f_x + l*(_PI_projector(f_x) - f_x)
189
190
        #print(numpy.amax(numpy.abs(f_x - f_x_old)))

Edward Andò's avatar
Edward Andò committed
191
        if GRAPH:
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),'.')
195
            plt.pause(0.00001)
196
197
198
        count += 1

        if count > maxIterations:
Benjy Marks's avatar
Benjy Marks committed
199
200
            # f_x = numpy.zeros_like(f_x) # NOTE: THIS WILL RETURN NO MATCHES WHEN TOMOPACK DOESN'T CONVERGE
            print('\nKILLED LOOP')
201
            break
Benjy Marks's avatar
Benjy Marks committed
202
203
    # plt.show()
    # print('Took ' + str(count) + ' iterations to complete.')
204

205
206
    if GRAPH:
        #not_k_trust_ind = numpy.where(~k_trust)
207

208
        plt.figure(figsize=[8,6])
209

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()
215

216
217
218
219
        plt.subplot(332)
        plt.title(r'$\psi(x)$')
        plt.imshow(psiMM, vmin=radioMM.min(), vmax=radioMM.max())
        plt.colorbar()
220

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')
226

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()
231

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()
236

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)
248

Edward Andò's avatar
Edward Andò committed
249

250
251
252
253
        plt.subplot(337)
        plt.title(r'Estimated $|f(k)|$')
        plt.imshow(numpy.abs(f_k))
        plt.colorbar()
Edward Andò's avatar
Edward Andò committed
254

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ò's avatar
Edward Andò committed
260

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ò's avatar
Edward Andò committed
289
290
291



292
293
294
        plt.subplots_adjust(hspace=0.5)
        plt.show()
        #plt.savefig('tt3D.png',dpi=200)
Edward Andò's avatar
Edward Andò committed
295
296

    #return pPredictedIJ
297
    return numpy.real(f_x)
Edward Andò's avatar
Edward Andò committed
298

299

Edward Andò's avatar
Edward Andò committed
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

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ò's avatar
Edward Andò committed
308

309
310
311
        debug : bool, optional
            Show debug graphs (especially in the maximum filtering)
            Default = False
Edward Andò's avatar
Edward Andò committed
312

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ò's avatar
Edward Andò committed
318

319
    f_x = _filter_maxima(f_x, debug=debug)
Edward Andò's avatar
Edward Andò committed
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]]
325
        # print(f_x[pos[0],pos[1]])
Edward Andò's avatar
Edward Andò committed
326
327
328
        for i in range(int(val)):
            pPredictedIJ.append([pos[1], pos[0]])
    pPredictedIJ = numpy.array(pPredictedIJ)
329
330

    return pPredictedIJ
331
332


Edward Andò's avatar
Edward Andò committed
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ò's avatar
Edward Andò committed
363
364
    tifffile.imsave(cacheFile, fXseries.astype('float'))
    print(f"saved {cacheFile}")
Edward Andò's avatar
Edward Andò committed
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('<f4'))
        #tifffile.imsave(cachePsiFile, psiXseries.astype('<f4'))
        #print("done.")

Edward Andò's avatar
Edward Andò committed
410
411
    #L_x  = 20 # TODO: SCALING IN X DIRECTION SHOULD BE A FUNCTION OF THE CONE ANGLE
    #L_yz =  2 # TODO: THIS SHOULD BE A FUNCTION OF THE PIXELS PER RADIUS
Edward Andò's avatar
Edward Andò committed
412

Edward Andò's avatar
Edward Andò committed
413
414
415
    #struct = psiXseries[(psiXseries.shape[0])//2 -  L_x:(psiXseries.shape[0])//2 + L_x  + 1,
                        #(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]
Edward Andò's avatar
Edward Andò committed
416

Edward Andò's avatar
Edward Andò committed
417
418
419
420
    #fXconvolvedSeries = scipy.ndimage.convolve(fXseries,struct/struct.sum())
    ##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'))
Edward Andò's avatar
Edward Andò committed
421
422


Edward Andò's avatar
Edward Andò committed
423
    #binaryPeaks = fXconvolvedSeries > massThreshold
Edward Andò's avatar
Edward Andò committed
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ò's avatar
Edward Andò committed
428
429
    fXseriesMaximumFiltered = scipy.ndimage.maximum_filter(fXseries,
                                                              size=(3*numpy.int(numpy.floor(radiusMM/CORxDelta)),
Edward Andò's avatar
Edward Andò committed
430
                                                                    numpy.int(numpy.floor(radiusMM/pixelSizeMM*zoomLevel)),
Edward Andò's avatar
Edward Andò committed
431
432
433
434
435
                                                                    numpy.int(numpy.floor(radiusMM/pixelSizeMM*zoomLevel)))
                                                              )
    allPeaks = fXseries == fXseriesMaximumFiltered
    masses = allPeaks*fXseries

Edward Andò's avatar
Edward Andò committed
436
437
438
439
    if verbose:
        tifffile.imsave(cacheFile[:-4] + '_masses.tif', masses.astype('<f4'))
        tifffile.imsave(cacheFile[:-4] + '_peaks.tif', allPeaks.astype('<f4'))
        tifffile.imsave(cacheFile[:-4] + '_fXseries.tif', fXseries.astype('<f4'))
Edward Andò's avatar
Edward Andò committed
440
        tifffile.imsave(cacheFile[:-4] + '_fXseriesMaximumFiltered.tif', fXseriesMaximumFiltered.astype('<f4'))
Edward Andò's avatar
Edward Andò committed
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484

    if scanFixedNumber:
        # get the indices of all of the peaks, from highest to lowest
        sortedPeakIndices = numpy.argsort(masses, axis=None)[::-1]
        # print(sortedPeakIndices.shape)
        # get just the first scanFixedNumber of those and put them into a scanFixedNumber x 3 array
        peaksCORxPOSnJI = numpy.vstack(numpy.unravel_index(sortedPeakIndices[:scanFixedNumber], masses.shape)).T
        # print(peaksCORxPOSnJI.shape)
    else:
        filteredPeaks = masses > massThreshold
        peaksCORxPOSnJI = numpy.argwhere(filteredPeaks)

        if verbose: tifffile.imsave(cacheFile[:-4] + '_filteredPeaks.tif', filteredPeaks.astype('<f4'))
    # print(peaksCORxPOSnJI)

    ###############################################################
    ### Now we have guesses for all particle according to detector
    ###   (IJ) and position along the X-scanning direction
    ### We're going to convert that to spatial XYZ
    ###############################################################
    print("\nConverting tomopack x-scan to 3D positions\n")
    ## Convert to XYZ in space and mm
    positionsXYZmm = numpy.zeros([peaksCORxPOSnJI.shape[0], 3])

    for i in range(positionsXYZmm.shape[0]):
        # X -- look up which CORx slice the maximum falls in, this could be interpolated instead of rounded
        positionsXYZmm[i,0] = CORxPositions[int(numpy.round(peaksCORxPOSnJI[i,0]))]

        # detector I gives real position Y in mm
        yPosDetMM = -1*(peaksCORxPOSnJI[i,2] - radio.shape[1]/2.0)*pixelSizeMM

        # detector J gives real position Z in mm
        zPosDetMM = -1*(peaksCORxPOSnJI[i,1] - radio.shape[0]/2.0)*pixelSizeMM

        # And now scale down by zoom factor
        # Y
        positionsXYZmm[i,1] = yPosDetMM * ( positionsXYZmm[i,0] / sourceDetectorDistMM )
        # Z
        positionsXYZmm[i,2] = zPosDetMM * ( positionsXYZmm[i,0] / sourceDetectorDistMM )

    print(f"\ntomopackDivergentScanTo3DPositions(): I'm returning {positionsXYZmm.shape[0]} 3D positions.\n")
    return positionsXYZmm


Benjy Marks's avatar
Benjy Marks committed
485
486
def tomopackDivergentScanTo3DPositions(radioMM, radiusMM,
                                       CORxMin=None, CORxMax=None, CORxNumber=100,
Edward Andò's avatar
Edward Andò committed
487
                                       psiXseries = None,
488
                                       massThreshold=0.12,
Benjy Marks's avatar
Benjy Marks committed
489
490
                                       scanPersistenceThresholdRadii=None,
                                       scanFixedNumber=None,
Benjy Marks's avatar
Benjy Marks committed
491
                                       scanPersistenceThreshold=7, maxIterations=50,
Benjy Marks's avatar
Benjy Marks committed
492
493
                                       sourceDetectorDistMM=100, pixelSizeMM=0.1, l=0.2, kTrustRatio=0.7, useCache=True,
                                       numCores=1,
494
                                       blur=0.0,
Benjy Marks's avatar
Benjy Marks committed
495
496
                                       cacheFile='fXseries.tif',
                                       verbose=False):
497
    """
498
499
500
    This function takes in a single divergent projection, and will run tomopack
    generating different psis by varying their position as Centre Of Rotation (COR) in the x-direction,
    from CORxMin to CORxMax in CORxNumber steps.
Edward Andò's avatar
Edward Andò committed
501

502
    The resulting series of indicator functions is analysed and a 3D position guess for
503
    all identified spheres is returned.
Edward Andò's avatar
Edward Andò committed
504

505
506
507
508
509
510
    Parameters
    ----------
        massThreshold : float (optional)
            Threshold for the accepting the result of the convolution of the indicator function.
            Deafult = 0.5

511
512
513
514
        blur : float, optional
            sigma of blur to pass to scipy.ndimage.gaussian_filter to
            blur the radiograph at the end of everything

Edward Andò's avatar
Edward Andò committed
515
        # scanPersistenceThresholdRadii : float (optional)
Benjy Marks's avatar
Benjy Marks committed
516
            # How much +/- in radii in the scanning direction should the particle's indicator function be above massThreshold? If not passed, scanPersistenceThreshold is used
517
518
519
520
521
522

    Returns
    -------
        positionsXYZmm : 2D numpy array
    """
    assert(len(radioMM.shape) == 2), "detectSpheres.tomopackDivergentScanTo3DPositions(): need 2D array"
Edward Andò's avatar
Edward Andò committed
523

524
525
526
527
    if CORxMin is None: CORxMin = CORx-3*radiusMM
    if CORxMax is None: CORxMax = CORx+3*radiusMM

    CORxPositions = numpy.linspace(CORxMin, CORxMax, CORxNumber)
528
    CORxDelta = CORxPositions[1] - CORxPositions[0]
529

Edward Andò's avatar
Edward Andò committed
530
531
532
    loadedCache = False

    if useCache:
533
534
        cachePsiFile = cacheFile[:-4] + '_psi.tif'
        if os.path.isfile(cacheFile) and os.path.isfile(cachePsiFile):
Edward Andò's avatar
Edward Andò committed
535
536
            print("Loading previous indicator functions... ", end="")
            fXseries = tifffile.imread(cacheFile)
537
538
539
540
541
542
            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.")
Edward Andò's avatar
Edward Andò committed
543
544
545
        else:
            print('No cached indicator functions found. Generating them now to cache.')
    if not loadedCache:
546
        fXseries = numpy.zeros((len(CORxPositions), radioMM.shape[0], radioMM.shape[1]))
547
        psiXseries = numpy.zeros_like(fXseries)
548
        psiMMseries = numpy.zeros_like(fXseries)
549
550
551
552
553

        psiRefMM = radioSphere.projectSphere.projectSphereMM(numpy.array([[(CORxMax+CORxMin)/2., 0., 0.]]),
                                                          numpy.array([radiusMM]),
                                                          detectorResolution=radioMM.shape,
                                                          pixelSizeMM=pixelSizeMM,
554
555
                                                          sourceDetectorDistMM=sourceDetectorDistMM,
                                                          blur=blur)
556
557
558
559
560
561
562
563

        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,
564
565
                                                              sourceDetectorDistMM=sourceDetectorDistMM,
                                                              blur=blur)
566

567
            fXseries[posN] = radioSphere.detectSpheres.tomopack(radioMM, psiMM, GRAPH=0, maxIterations=maxIterations, l=l, kTrustRatio=kTrustRatio)
568
            psiXseries[posN] = radioSphere.detectSpheres.tomopack(psiRefMM, psiMM, GRAPH=0, maxIterations=maxIterations, l=l, kTrustRatio=kTrustRatio)
569
            psiMMseries[posN] = psiMM
570

571
572
573
574
    if useCache and not loadedCache:
        print("Saving indicator functions for next time... ", end="")
        tifffile.imsave(cacheFile, fXseries.astype('<f4'))
        tifffile.imsave(cachePsiFile, psiXseries.astype('<f4'))
575
        tifffile.imsave(cachePsiFile[:-4]+'_MM.tif', psiMMseries.astype('<f4'))
576
        print("done.")
577

578

579
    use_old_method = False
580
581
582
583
584
585
586
    if use_old_method:
        ### Step 1. Maximum filtering... Trying only in the detector direction, in order not to bridge over the bad
        #           COR x positions where there are weird resonances
        ### Step 2. Lump together objects in the space of IJ and X scan
        binaryPeaks = scipy.ndimage.convolve(fXseries,_kernel) > massThreshold
        structure = numpy.ones([5,3,3])
        binaryPeaks = scipy.ndimage.binary_closing(binaryPeaks,structure=structure) # remove any gaps
587
        if verbose: tifffile.imsave('closed_peaks.tif', binaryPeaks.astype('<f4'))
588
    else:
Edward Andò's avatar
Edward Andò committed
589
590
        L_x  = 20 # TODO: SCALING IN X DIRECTION SHOULD BE A FUNCTION OF THE CONE ANGLE
        L_yz =  2 # TODO: THIS SHOULD BE A FUNCTION OF THE PIXELS PER RADIUS
591

592
593
594
595
596
        struct = psiXseries[(psiXseries.shape[0])//2 -  L_x:(psiXseries.shape[0])//2 + L_x  + 1,
                            (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())
597
598
599
600
601
        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'))


602
        binaryPeaks = fXconvolvedSeries > massThreshold
603

604
    get_peaks_from_persistence = False
605
    if get_peaks_from_persistence:
Edward Andò's avatar
Edward Andò committed
606

607
608
        labelledPeaks = scipy.ndimage.label(binaryPeaks)[0]
        labelledPeaksBB = scipy.ndimage.find_objects(labelledPeaks)
Edward Andò's avatar
Edward Andò committed
609

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))
613

614
615
616
617
        try:
            import spam.label
            #print("Doing it with spam")
            BB = spam.label.boundingBoxes(labelledPeaks)
Edward Andò's avatar
Edward Andò committed
618
619

            # save x-y-z positions where x is the slice number in the scan direction and y, z pixel positions on detector
620
            peaksCORxPOSnJI[:,0:3] = spam.label.centresOfMass(labelledPeaks, boundingBoxes=BB)[1:]
Edward Andò's avatar
Edward Andò committed
621
622

            # also add in how many slices this is alive for
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:
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
663
        masses = allPeaks*fXconvolvedSeries
664
        if verbose:
665
666
            tifffile.imsave(cacheFile[:-4] + '_masses.tif', masses.astype('<f4'))
            tifffile.imsave(cacheFile[:-4] + '_peaks.tif', allPeaks.astype('<f4'))
667
            tifffile.imsave(cacheFile[:-4] + '_fXconvolvedSeries.tif', fXconvolvedSeries.astype('<f4'))
668
            tifffile.imsave(cacheFile[:-4] + '_fXconvolvedMaximumFiltered.tif', fXconvolvedMaximumFiltered.astype('<f4'))
669
670
671
672

        if scanFixedNumber:
            # get the indices of all of the peaks, from highest to lowest
            sortedPeakIndices = numpy.argsort(masses, axis=None)[::-1]
673
            # print(sortedPeakIndices.shape)
674
675
            # get just the first scanFixedNumber of those and put them into a scanFixedNumber x 3 array
            peaksCORxPOSnJI = numpy.vstack(numpy.unravel_index(sortedPeakIndices[:scanFixedNumber], masses.shape)).T
676
            # print(peaksCORxPOSnJI.shape)
677
678
679
        else:
            filteredPeaks = masses > massThreshold
            peaksCORxPOSnJI = numpy.argwhere(filteredPeaks)
680
681

            if verbose: tifffile.imsave(cacheFile[:-4] + '_filteredPeaks.tif', filteredPeaks.astype('<f4'))
682
        # print(peaksCORxPOSnJI)
683

684
685
686



687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
    ###############################################################
    ### Now we have guesses for all particle according to detector
    ###   (IJ) and position along the X-scanning direction
    ### We're going to convert that to spatial XYZ
    ###############################################################
    print("\nConverting tomopack x-scan to 3D positions\n")
    ## Convert to XYZ in space and mm
    positionsXYZmm = numpy.zeros([peaksCORxPOSnJI.shape[0], 3])

    for i in range(positionsXYZmm.shape[0]):
        # X -- look up which CORx slice the maximum falls in, this could be interpolated instead of rounded
        positionsXYZmm[i,0] = CORxPositions[int(numpy.round(peaksCORxPOSnJI[i,0]))]

        # detector I gives real position Y in mm
        yPosDetMM = -1*(peaksCORxPOSnJI[i,2] - radioMM.shape[1]/2.0)*pixelSizeMM

        # detector J gives real position Z in mm
        zPosDetMM = -1*(peaksCORxPOSnJI[i,1] - radioMM.shape[0]/2.0)*pixelSizeMM

        # And now scale down by zoom factor
        # Y
        positionsXYZmm[i,1] = yPosDetMM * ( positionsXYZmm[i,0] / sourceDetectorDistMM )
        # Z
        positionsXYZmm[i,2] = zPosDetMM * ( positionsXYZmm[i,0] / sourceDetectorDistMM )

    print(f"\ntomopackDivergentScanTo3DPositions(): I'm returning {positionsXYZmm.shape[0]} 3D positions.\n")
    return positionsXYZmm
714

715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
def removeParticle(positionsXYZmm,residual,radioMM,radiiMM,pixelSizeDetectorMM,zoomLevel,sourceObjectDistMM,verbose,GRAPH):
    if verbose: print('Removing a particle')
    # find the location of the highest peak on the detector panel. This is presumably the centroid.
    residualMMPeakIndices = numpy.unravel_index(numpy.argmax(residual, axis=None), residual.shape)
    if verbose: print(f'residualMMPeakIndices: {residualMMPeakIndices}')

    # define unit vector between source and peak location
    yPosDetMM = -1*(residualMMPeakIndices[1] - radioMM.shape[1]/2.0)*pixelSizeDetectorMM
    zPosDetMM = -1*(residualMMPeakIndices[0] - radioMM.shape[0]/2.0)*pixelSizeDetectorMM
    magnitude = numpy.sqrt(zoomLevel**2*sourceObjectDistMM**2 + yPosDetMM**2 + zPosDetMM**2)
    s = numpy.array([
        zoomLevel*sourceObjectDistMM/magnitude,
        yPosDetMM/magnitude,
        zPosDetMM/magnitude])
    if verbose: print(f's: {s}')

    # find distance of every particle from the line defined by this unit vector
    distances = numpy.linalg.norm(numpy.cross(positionsXYZmm,s),axis=1)
    if verbose: print(f'distances: {distances}')

    # remove the particle closest to the line
    closestParticleIndex = numpy.argmin(distances)
    positionsXYZmm = numpy.delete(positionsXYZmm, closestParticleIndex, axis=0)

    p_f_x = radioSphere.projectSphere.projectSphereMM(positionsXYZmm,
                                                      radiiMM[0]*numpy.ones(len(positionsXYZmm)),
                                                      sourceDetectorDistMM=zoomLevel*sourceObjectDistMM,
                                                      pixelSizeMM=pixelSizeDetectorMM,
                                                      detectorResolution=radioMM.shape)
    residual = p_f_x - radioMM

    if GRAPH:
        plt.imshow(residual)
        plt.show()

    return positionsXYZmm, residual

752
def addParticle(*args, **kwargs):
Benjy Marks's avatar
Benjy Marks committed
753
754
    # return addParticleRaster(*args, **kwargs)
    return addParticleSensitivity(*args, **kwargs)
755

Edward Andò's avatar
Edward Andò committed
756

757
def addParticleRaster(positionsXYZmm,residual,radioMM,radiiMM,pixelSizeDetectorMM,zoomLevel,sourceObjectDistMM,CORxMin, CORxMax,CORxNumber,verbose,GRAPH):
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
    if verbose: print('Adding a particle')
    # find the location of the highest peak on the detector panel. This is presumably the centroid.
    residualMMPeakIndices = numpy.unravel_index(numpy.argmin(residual, axis=None), residual.shape)
    if verbose: print(f'residualMMPeakIndices: {residualMMPeakIndices}')

    # define unit vector between source and peak location
    yPosDetMM = -1*(residualMMPeakIndices[1] - radioMM.shape[1]/2.0)*pixelSizeDetectorMM
    zPosDetMM = -1*(residualMMPeakIndices[0] - radioMM.shape[0]/2.0)*pixelSizeDetectorMM
    magnitude = numpy.sqrt(zoomLevel**2*sourceObjectDistMM**2 + yPosDetMM**2 + zPosDetMM**2)
    s = numpy.array([
        zoomLevel*sourceObjectDistMM/magnitude,
        yPosDetMM/magnitude,
        zPosDetMM/magnitude])
    if verbose: print(f's: {s}')

773
    # trying to find an optimal solution by doing a raster scan and looking for minimal residual
774
    x_test = numpy.linspace(CORxMin,CORxMax,CORxNumber)
775
776
    best_index = 0
    best_residual = numpy.inf
Benjy Marks's avatar
Benjy Marks committed
777
778
779
780
781
    ref_projection = radioSphere.projectSphere.projectSphereMM(positionsXYZmm,
                                                      radiiMM[0]*numpy.ones(len(positionsXYZmm)),
                                                      sourceDetectorDistMM=zoomLevel*sourceObjectDistMM,
                                                      pixelSizeMM=pixelSizeDetectorMM,
                                                      detectorResolution=radioMM.shape)
782

Benjy Marks's avatar
Benjy Marks committed
783
784
785
786
787
788
    limits = radioSphere.projectSphere.singleSphereToDetectorPixelRange(s*x_test[0],
                                                                        radiiMM[0],
                                                                        radiusMargin=0.1,
                                                                        sourceDetectorDistMM=zoomLevel*sourceObjectDistMM,
                                                                        pixelSizeMM=pixelSizeDetectorMM,
                                                                        detectorResolution=radioMM.shape)
Benjy Marks's avatar
Benjy Marks committed
789
790
791

    ref_projection_crop = ref_projection[limits[0,0]:limits[0,1], limits[1,0]:limits[1,1]]
    radioMM_crop = radioMM[limits[0,0]:limits[0,1], limits[1,0]:limits[1,1]]
Benjy Marks's avatar
Benjy Marks committed
792
793
794
    for i,x in enumerate(x_test):
        single_particle_projection = radioSphere.projectSphere.projectSphereMM(numpy.expand_dims(s*x,axis=0),
                                                          numpy.expand_dims(radiiMM[0],axis=0),
795
796
                                                          sourceDetectorDistMM=zoomLevel*sourceObjectDistMM,
                                                          pixelSizeMM=pixelSizeDetectorMM,
Benjy Marks's avatar
Benjy Marks committed
797
798
799
800
                                                          detectorResolution=radioMM.shape,
                                                          ROIcentreMM=s*x_test[0],
                                                          ROIradiusMM=radiiMM[0])

Benjy Marks's avatar
Benjy Marks committed
801
        residual = ref_projection_crop + single_particle_projection - radioMM_crop
802
        # print(i, (residual**2).sum(), best_residual, best_index)
Benjy Marks's avatar
Benjy Marks committed
803

804
805
806
        if (residual**2).sum() < best_residual:
            best_index = i
            best_residual = (residual**2).sum()
Benjy Marks's avatar
Benjy Marks committed
807
808
809
810
811

        # plt.ion()
        # plt.title(i)
        # plt.imshow(residual)
        # plt.pause(0.001)
Benjy Marks's avatar
Benjy Marks committed
812
    bestPositionXYZmm = s*x_test[best_index]
813
814
815

    print(f'Best location at {best_index}-th x value: {x_test[best_index]}')

Benjy Marks's avatar
Benjy Marks committed
816
817
    optimise = True
    if optimise:
Edward Andò's avatar
Edward Andò committed
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
        import radiosphere.optimsePositions
        bestPositionXYZmm = radioSphere.optimisePositions.optimiseSensitivityFields(radioMM,
                                                                                    numpy.expand_dims(bestPositionXYZmm,axis=0), # try the middle of the sample
                                                                                    numpy.expand_dims(radiiMM[0],axis=0),
                                                                                    perturbationMM=(0.01, 0.01, 0.01),
                                                                                    # perturbationMM=(0.5, 0.25, 0.25),
                                                                                    # perturbationMM=(1, 0.5, 0.5),
                                                                                    # perturbationMM=(3, 1, 1),
                                                                                    minDeltaMM=0.0001,
                                                                                    iterationsMax=500,
                                                                                    sourceDetectorDistMM=zoomLevel*sourceObjectDistMM,
                                                                                    pixelSizeMM=pixelSizeDetectorMM,
                                                                                    detectorResolution=radioMM.shape,
                                                                                    verbose=False,
                                                                                    # GRAPH=True,
                                                                                    # NDDEM_output=True
                                                                                    )
Benjy Marks's avatar
Benjy Marks committed
835
836

    positionsXYZmm = numpy.vstack([positionsXYZmm, bestPositionXYZmm])
837
838
839
840
841
842
843
    p_f_x = radioSphere.projectSphere.projectSphereMM(positionsXYZmm,
                                                      radiiMM[0]*numpy.ones(len(positionsXYZmm)),
                                                      sourceDetectorDistMM=zoomLevel*sourceObjectDistMM,
                                                      pixelSizeMM=pixelSizeDetectorMM,
                                                      detectorResolution=radioMM.shape)
    residual = p_f_x - radioMM

844
845
846
    if GRAPH:
        plt.imshow(residual)
        plt.show()
847

848
    return positionsXYZmm, residual
849

Edward Andò's avatar
Edward Andò committed
850

851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
def addParticleSensitivity(positionsXYZmm,residual,radioMM,radiiMM,pixelSizeDetectorMM,zoomLevel,sourceObjectDistMM,CORxMin, CORxMax,CORxNumber,verbose,GRAPH):
    if verbose: print('Adding a particle')
    # find the location of the highest peak on the detector panel. This is presumably the centroid.
    residualMMPeakIndices = numpy.unravel_index(numpy.argmin(residual, axis=None), residual.shape)
    if verbose: print(f'residualMMPeakIndices: {residualMMPeakIndices}')

    # define unit vector between source and peak location
    yPosDetMM = -1*(residualMMPeakIndices[1] - radioMM.shape[1]/2.0)*pixelSizeDetectorMM
    zPosDetMM = -1*(residualMMPeakIndices[0] - radioMM.shape[0]/2.0)*pixelSizeDetectorMM
    magnitude = numpy.sqrt(zoomLevel**2*sourceObjectDistMM**2 + yPosDetMM**2 + zPosDetMM**2)
    s = numpy.array([
        zoomLevel*sourceObjectDistMM/magnitude,
        yPosDetMM/magnitude,
        zPosDetMM/magnitude])
    if verbose: print(f's: {s}')

    # trying to find an optimal solution by doing a raster scan and looking for minimal residual
Edward Andò's avatar
Edward Andò committed
868
    positionXYZmmOpt = radioSphere.optimisePositions.optimiseSensitivityFields(   radioMM,
869
870
                            numpy.expand_dims(s*sourceObjectDistMM,axis=0), # try the middle of the sample
                            numpy.expand_dims(radiiMM[0],axis=0),
Benjy Marks's avatar
Benjy Marks committed
871
                            perturbationMM=(0.001, 0.001, 0.001),
872
873
874
875
876
877
                            minDeltaMM=0.0001,
                            iterationsMax=500,
                            sourceDetectorDistMM=zoomLevel*sourceObjectDistMM,
                            pixelSizeMM=pixelSizeDetectorMM,
                            detectorResolution=radioMM.shape,
                            verbose=False,
Benjy Marks's avatar
Benjy Marks committed
878
                            NDDEM_output=True
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
                            )



    positionsXYZmm = numpy.vstack([positionsXYZmm, positionXYZmmOpt])
    p_f_x = radioSphere.projectSphere.projectSphereMM(positionsXYZmm,
                                                      radiiMM[0]*numpy.ones(len(positionsXYZmm)),
                                                      sourceDetectorDistMM=zoomLevel*sourceObjectDistMM,
                                                      pixelSizeMM=pixelSizeDetectorMM,
                                                      detectorResolution=radioMM.shape)
    residual = p_f_x - radioMM

    if GRAPH:
        plt.imshow(residual)
        plt.show()

    return positionsXYZmm, residual

Edward Andò's avatar
Edward Andò committed
897

898
def cleanDivergentScan(positionsXYZmm,radioMM,radiiMM,zoomLevel,sourceObjectDistMM,pixelSizeDetectorMM,CORxMin,CORxMax, CORxNumber,verbose=False,GRAPH=False):
899
900
901
902
903
    p_f_x = radioSphere.projectSphere.projectSphereMM(positionsXYZmm,
                                                      radiiMM[0]*numpy.ones(len(positionsXYZmm)),
                                                      sourceDetectorDistMM=zoomLevel*sourceObjectDistMM,
                                                      pixelSizeMM=pixelSizeDetectorMM,
                                                      detectorResolution=radioMM.shape)
904

905
906
907
908
909
910
911
912
913
914
    residual = p_f_x - radioMM
    print(f'Maximum |residual| {numpy.abs(residual).max()}')
    # Found an extra particle that shouldn't be there
    removed = 0
    added = 0
    # remove any overlaps
    while residual.max() > 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)
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)
923
        added +=1
Benjy Marks's avatar
Benjy Marks committed
924
925

    print(f'Removed {removed} particles and added {added}.')
926
    return positionsXYZmm
927

Edward Andò's avatar
Edward Andò committed
928

Benjy Marks's avatar
Benjy Marks committed
929
930
def calculateErrors(posXYZa,posXYZb,radiiMM,verbose=True):
    distances = distance.cdist(posXYZa,posXYZb)
931
932
    distance_threshold = radiiMM[0]

Benjy Marks's avatar
Benjy Marks committed
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ò's avatar
Edward Andò committed
935
    if number_lost > 0: print(f"detectSpheres.calculateErrors(): number_lost = {number_lost}")
936

Benjy Marks's avatar
Benjy Marks committed
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
940

Benjy Marks's avatar
Benjy Marks committed
941
942
    err_mean    = numpy.mean(min_valid_errors)
    err_std     = numpy.std( min_valid_errors)
943
944

    return err_mean, err_std, number_lost