Skip to content

GitLab

Commit af16e5db authored by Olga Stamati's avatar Olga Stamati
Browse files

clean and document tomopack()

parent 6af2860b
Pipeline #79967 passed with stage
in 2 minutes and 8 seconds
Hide whitespace changes
Inline Side-by-side
tools/detectSpheres.py
View file @ af16e5db
......@@ -15,15 +15,14 @@ import tifffile
#_kernel = numpy.array([[[1,2,1],[2,4,2],[1,2,1]]])/16. # gaussian kernel
_kernel = numpy.array([[1,2,1], [2,4,2], [1,2,1]])/16. # gaussian kernel 2D
# non-naiive approach
def _PI_projector(f):
# helper projector function operating on fx (approximation of "indicator" function)
def _PIprojector(f):
g = f.copy()
g[f<0.5] = 0.
# g[g.imag < 1e-14] = numpy.round(g[g.imag < 1e-14])
# g.real = numpy.round(g.real)
g = numpy.round(numpy.abs(g))
g[f<0.] = 0. # discard negative values
g = numpy.round(g) # round to the nearest integer
return g
# 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
# Why not directly look at maxima in the masses, this avoids the
def _filter_maxima(f, debug=False, removeNegatives=False):
......@@ -75,7 +74,7 @@ def _filter_maxima(f, debug=False, removeNegatives=False):
return masses
def tomopack(radioMM, psiMM, maxIterations=50, l=0.1, epsilon='iterative', kTrustRatio=0.75, GRAPH=False):
def tomopack(radioMM, psiMM, maxIterations=50, l=0.5, epsilon='iterative', kTrustRatio=0.75, verbose=False, graphShow=False):
'''
'tomopack' FFT-based algorithm for sphere identification by Stéphane Roux.
the FFT of psiMM as a structuring element to pick up projected spheres in radioMM
......@@ -84,209 +83,146 @@ def tomopack(radioMM, psiMM, maxIterations=50, l=0.1, epsilon='iterative', kTrus
Parameters
----------
radioMM : 2D numpy array of floats
radiography containing "distance" projections in mm
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,
as radioMM -- this is the structuring element used for FFT recognition.
A (synthetic) projection in mm of a single sphere in the centre of the detector
This is the structuring element used for FFT recognition.
Should be same size as radioMM
maxIterations : int (optional)
Number of iterations to run the detection for.
Number of iterations to run the detection for
Default = 50
l : float (optional)
"lambda" parameter which controls relaxation in the iterations.
Default = 0.1
"lambda" parameter which controls relaxation in the iterations
Default = 0.5
epsilon : float (optional), or 'iterative'
trust cutoff in psi.
Trust cutoff wavelengths in psi.
If the radiograph is in mm so is this parameter
Default = 'iterative'
Default = 'iterative' -- updating until `kTrustRatio` (below) is reached
GRAPH : bool (optional)
VERY noisily make graphs?
kTrustRatio : float (optional)
Ratio of cutoff psi wavelengths
Default = 0.75
verbose : bool, optional
Get to know what the function is really thinking
Default = False
graphShow : bool (optional)
Show evolution of fx during iterations, and final one together with radio and psi
Default = False
Returns
-------
f_x : 2D numpy array
same size as radioMM, high values where we think particles are.
Consider using indicatorFunctionToDetectorPositions()
fx : 2D numpy array
The approximation of the `indicator` function
Same size as radioMM, with high values where we think spheres centres are
Consider using `indicatorFunctionToDetectorPositions()` to get the IJ detector coordinates of the centres
'''
#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(radioMM.shape) == 2), "detectSpheres.tomopack(): radio should be a 2D array"
assert(len(psiMM.shape) == 2), "detectSpheres.tomopack(): psi should be a 2D array"
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"
radioMM_FFT = numpy.fft.fft2(radioMM)
# compute FFT of input radio
radioMMFFT = numpy.fft.fft2(radioMM)
# define the psi function and its FFT -- the structuring element for SR
# Start with the projection of a centred sphere
# compute FFT of input psi (structuring element or shape function)
#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))
# avoids the HACK of shifting fx later
psiMMFFT = numpy.fft.fft2(numpy.fft.fftshift(psiMM))
# This is comparable to Figure 7 in TomoPack
# 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()
## naiive approach
# Naiive approach to compute fx by a direct division and deconvolution
with numpy.errstate(divide='ignore', invalid='ignore'):
f_k_naiive = radioMM_FFT/psiMM_FFT
f_x_naiive = numpy.fft.ifft2(f_k_naiive)
fkNaiive = radioMMFFT/psiMMFFT
fxNaiive = numpy.fft.ifft2(fkNaiive)
# Prepare for iterations
f_x = numpy.zeros_like(radioMM)
fx = numpy.zeros_like(radioMM)
if epsilon == 'iterative':
epsilon = 1e-5
k_trust = numpy.abs(psiMM_FFT) > epsilon
while k_trust.sum()/k_trust.shape[0]/k_trust.shape[1] > kTrustRatio:
# print(epsilon)
kTrust = numpy.abs(psiMMFFT) > epsilon
while kTrust.sum()/kTrust.shape[0]/kTrust.shape[1] > kTrustRatio:
epsilon *= 1.5
k_trust = numpy.abs(psiMM_FFT) > epsilon
f_x_old = epsilon*numpy.ones_like(radioMM)
kTrust = numpy.abs(psiMMFFT) > epsilon
fxOld = 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]))
count = 0
if GRAPH: plt.ion()
fxOld = epsilon*numpy.ones_like(radioMM)
kTrust = numpy.abs(psiMMFFT) > epsilon
if verbose: print(f"kTrust maintains {100*kTrust.sum()/kTrust.shape[0]/kTrust.shape[1]}% of wavenumubers")
it = 0
# 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
if GRAPH: plt.clf()
f_x_old = f_x.copy()
f_k = numpy.fft.fft2(f_x)
# if GRAPH: plt.plot(f_k,'k.')
f_k[k_trust] = f_k_naiive[k_trust]
# if GRAPH: plt.plot(f_k,'r.')
f_x = numpy.fft.ifft2(f_k)
#f_x = numpy.fft.fftshift(numpy.fft.ifft2(f_k)) # importing HACK from 1D version
f_x = f_x + l*(_PI_projector(f_x) - f_x)
#print(numpy.amax(numpy.abs(f_x - f_x_old)))
if GRAPH:
plt.imshow(numpy.abs(f_x), vmin=0)#, vmax=1)
# Objective: get a good indicator function fx
while numpy.linalg.norm(fx - fxOld) > 1e-8: # NOTE: Using a different epsilon here (OS: could be an input perhaps)
if graphShow: plt.clf()
fxOld = fx.copy()
fk = numpy.fft.fft2(fx)
fk[kTrust] = fkNaiive[kTrust]
#fx = numpy.fft.fftshift(numpy.fft.ifft2(fk)) # importing HACK from 1D version
fx = numpy.fft.ifft2(fk)
fx = fx + l*(_PIprojector(fx) - fx)
if graphShow:
plt.suptitle(f'Iteration Number = {it}', fontsize=10)
plt.imshow(numpy.real(fx), vmin=0, vmax=1)
plt.colorbar()
# plt.plot(numpy.real(f_x),numpy.imag(f_x),'.')
plt.pause(0.00001)
count += 1
plt.pause(0.001)
it += 1
if count > maxIterations:
# f_x = numpy.zeros_like(f_x) # NOTE: THIS WILL RETURN NO MATCHES WHEN TOMOPACK DOESN'T CONVERGE
print('\nKILLED LOOP')
if it > maxIterations:
if verbose: print('\tNo convergence before max iterations"')
# fx = numpy.zeros_like(fx) # NOTE: Should we return a zero indicator function?
break
# plt.show()
# print('Took ' + str(count) + ' iterations to complete.')
if GRAPH:
#not_k_trust_ind = numpy.where(~k_trust)
plt.figure(figsize=[8,6])
if graphShow:
from matplotlib.colors import LogNorm
plt.figure(figsize=[8, 6])
plt.subplot(331)
plt.subplot(321)
plt.title(r'$p(x)$')
plt.imshow(radioMM)
# plt.plot(pPredictedIJ[:,0], pPredictedIJ[:,1], 'rx')
plt.colorbar()
plt.subplot(332)
plt.subplot(322)
plt.title(r'$\psi(x)$')
plt.imshow(psiMM, vmin=radioMM.min(), vmax=radioMM.max())
plt.colorbar()
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')
plt.subplot(334)
plt.subplot(323)
plt.title(r'$|\tilde{p}(k)|$')
plt.imshow(numpy.abs(radioMM_FFT))#, vmin=0, vmax=10)
plt.imshow(numpy.abs(numpy.fft.fftshift(radioMMFFT)), norm=LogNorm(vmin=0.1, vmax=10**5))
plt.colorbar()
plt.subplot(335)
plt.subplot(324)
plt.title(r'$|\tilde{\psi}(k)|$')
plt.imshow(numpy.abs(psiMM_FFT), vmin=0, vmax=10)
plt.imshow(numpy.abs(numpy.fft.fftshift(psiMMFFT)), norm=LogNorm(vmin=0.1, vmax=10**5))
plt.colorbar()
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)
plt.subplot(337)
plt.title(r'Estimated $|f(k)|$')
plt.imshow(numpy.abs(f_k))
plt.colorbar()
plt.subplot(338)
plt.subplot(325)
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.imshow(numpy.real(fx), vmin=0, vmax=1)
plt.colorbar()
#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()
psiMMFFTtrusted = psiMMFFT.copy()
psiMMFFTtrusted[~kTrust] = numpy.nan
plt.subplot(326)
plt.title(r'$|\tilde{\psi}(ktrust)|$')
plt.imshow(numpy.abs(numpy.fft.fftshift(psiMMFFTtrusted)), norm=LogNorm(vmin=0.1, vmax=10**5))
plt.colorbar()
plt.subplots_adjust(hspace=0.5)
plt.show()
#plt.savefig('tt3D.png',dpi=200)
#return pPredictedIJ
return numpy.real(f_x)
return numpy.real(fx)
def indicatorFunctionToDetectorPositions(f_x, scanFixedNumber=None, massThreshold=0.5, debug=False):
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment