Commit 6af2860b authored by Olga Stamati's avatar Olga Stamati
Browse files

tests for tomopack() and modify and tests for indicatorFunctionToDetectorPositions()

parent 06ce4084
Pipeline #79348 passed with stage
in 2 minutes and 42 seconds
"""
Simple tests on tomopack to make sure everything works as expected.
The layout is as per python's unittest but right it should just
be executed, and all the `self.assertEqual\True`s should pass.
If a new function is added it should be called at the bottom of the function
"""
import numpy
import unittest
import matplotlib.pyplot as plt
import radioSphere
class TestDetection(unittest.TestCase):
def tearDown(self):
try:
pass
except OSError:
pass
def test_tomopackOneSphere(self):
COR = 50
r = 1
detectorResolution = [400, 600]
# Case 1: simplest case with 1 sphere in the middle of the detector and parallel projection
p1 = radioSphere.projectSphere.projectSphereMM(numpy.array([[COR, 0, 0]]),
numpy.array([r]),
detectorResolution = detectorResolution,
sourceDetectorDistMM = numpy.inf)
psi = radioSphere.projectSphere.projectSphereMM(numpy.array([[COR, 0, 0]]),
numpy.array([r]),
detectorResolution = detectorResolution,
sourceDetectorDistMM = numpy.inf)
f_x1 = radioSphere.detectSpheres.tomopack(p1, psi)
peaksJI1 = radioSphere.detectSpheres.indicatorFunctionToDetectorPositions(f_x1, scanFixedNumber=1)
self.assertEqual(peaksJI1[0][0], detectorResolution[0]//2)
self.assertEqual(peaksJI1[0][1], detectorResolution[1]//2)
# Case 2: 1 sphere offset from the middle of the detector and parallel projection
p2 = radioSphere.projectSphere.projectSphereMM(numpy.array([[COR, -5.7, 7.4]]),
numpy.array([r]),
detectorResolution = detectorResolution,
sourceDetectorDistMM = numpy.inf)
f_x2 = radioSphere.detectSpheres.tomopack(p2, psi)
peaksJI2 = radioSphere.detectSpheres.indicatorFunctionToDetectorPositions(f_x2, massThreshold=0.1)
spherePosDetector2 = numpy.where(p2 == p2.max())
estimatedSpherePosDetector2 = numpy.where(f_x2 == f_x2.max())
# 1px max diff
self.assertTrue(abs(spherePosDetector2[0]-peaksJI2[0][0]) <= 1)
self.assertTrue(abs(spherePosDetector2[1]-peaksJI2[0][1]) <= 1)
# Case 3: divergent beam, simplest case with 1 sphere in the middle of the detector
p3 = radioSphere.projectSphere.projectSphereMM(numpy.array([[COR, 0, 0]]),
numpy.array([r]),
detectorResolution = detectorResolution,
sourceDetectorDistMM = 100)
psi3 = radioSphere.projectSphere.projectSphereMM(numpy.array([[COR, 0, 0]]),
numpy.array([r]),
detectorResolution = detectorResolution,
sourceDetectorDistMM = 100)
f_x3 = radioSphere.detectSpheres.tomopack(p3, psi3)
peaksJI3 = radioSphere.detectSpheres.indicatorFunctionToDetectorPositions(f_x3, scanFixedNumber=1)
self.assertEqual(peaksJI3[0][0], detectorResolution[0]//2)
self.assertEqual(peaksJI3[0][1], detectorResolution[1]//2)
# Case 4: divergent beam, 1 sphere offset from the middle
p4 = radioSphere.projectSphere.projectSphereMM(numpy.array([[COR, -1.2, -2.6]]),
numpy.array([r]),
detectorResolution = detectorResolution,
sourceDetectorDistMM = 100)
f_x4 = radioSphere.detectSpheres.tomopack(p4, psi3)
peaksJI4 = radioSphere.detectSpheres.indicatorFunctionToDetectorPositions(f_x4, scanFixedNumber=1)
spherePosDetector4 = numpy.where(p4 == p4.max())
# 1px max diff
self.assertTrue(abs(spherePosDetector4[0][0]-peaksJI4[0][0]) <= 1)
self.assertTrue(abs(spherePosDetector4[1][0]-peaksJI4[0][1]) <= 1)
# tomopack 10 random spheres with divergent beam
def test_tomopack10Spheres(self):
COR = 50
r = 1
detectorResolution = [100, 150]
sourceDetectorDistMM = 100
pixelSizeMM = 0.1
yMin, yMax = -1.5, 1.5
zMin, zMax = -2, 2
pos = numpy.zeros((10, 3))
pos[:, 0] += COR
pos[:, 1] = numpy.linspace(yMin, yMax, 10)
pos[:, 2] = zMin + r + (zMax - zMin - 2*r)*numpy.random.rand(10)
radii = numpy.ones_like(pos[:, 0])*r
p = radioSphere.projectSphere.projectSphereMM(pos,
radii,
detectorResolution = detectorResolution,
pixelSizeMM = pixelSizeMM,
sourceDetectorDistMM = sourceDetectorDistMM)
psi = radioSphere.projectSphere.projectSphereMM(numpy.array([[COR, 0, 0]]),
numpy.array([r]),
detectorResolution = detectorResolution,
pixelSizeMM = pixelSizeMM,
sourceDetectorDistMM = sourceDetectorDistMM)
f_x = radioSphere.detectSpheres.tomopack(p, psi)
peaksJI = radioSphere.detectSpheres.indicatorFunctionToDetectorPositions(f_x, scanFixedNumber=pos.shape[0])
posTomopackXYZ = numpy.zeros([peaksJI.shape[0], 3])
posTomopackXYZ[:, 0] += COR
# TODO: this should be a separate helper function for both a divergent and a parallel beam
for i in range(posTomopackXYZ.shape[0]):
# detector I gives real position Y in mm and scale by zoom
posTomopackXYZ[i, 1] = -1*(peaksJI[i,1] - p.shape[1]/2.0)*pixelSizeMM
posTomopackXYZ[i, 1] /= sourceDetectorDistMM/COR
# detector J gives real position Z in mm and scale by zoom
posTomopackXYZ[i, 2] = -1*(peaksJI[i,0] - p.shape[0]/2.0)*pixelSizeMM
posTomopackXYZ[i, 2] /= sourceDetectorDistMM/COR
# sort tomopack guess to calculate residual directly on positions
posTomopackXYZ = posTomopackXYZ[numpy.argsort(posTomopackXYZ[:, 1])]
self.assertTrue(numpy.linalg.norm(pos - posTomopackXYZ) < r)
if __name__ == "__main__":
unittest.main()
......@@ -12,7 +12,8 @@ import matplotlib.pyplot as plt
import tifffile
# _kernel = numpy.ones((1,3,3))/9. # square kernel
_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
_kernel = numpy.array([[1,2,1], [2,4,2], [1,2,1]])/16. # gaussian kernel 2D
# non-naiive approach
def _PI_projector(f):
......@@ -26,12 +27,12 @@ def _PI_projector(f):
# 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):
# Detect unpadded 2D image first:
if len(f.shape) == 2:
f = f[numpy.newaxis, ...]
twoD = True
else:
twoD = 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
f_abs = numpy.abs(f)
......@@ -39,49 +40,40 @@ def _filter_maxima(f, debug=False, removeNegatives=False):
# Remove negatives -- however keeping negatives can help sharpen the result of the convolution
if removeNegatives: f_abs[f<0] = 0
masses = scipy.ndimage.convolve(f_abs,_kernel)
fxConvolved = scipy.ndimage.convolve(f_abs, _kernel)
fxConvolvedMaxFiltered = scipy.ndimage.maximum_filter(fxConvolved, size=(3,3))
peaks = masses == scipy.ndimage.maximum_filter(masses,size=(1,3,3))
peaks = fxConvolved == fxConvolvedMaxFiltered
masses = peaks*fxConvolved
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.subplot(2, 3, 1)
plt.imshow(f)
plt.title('f_x')
plt.colorbar()
#plt.imshow(g)
#plt.colorbar()
#plt.title('Local filtered maxima')
#plt.show()
plt.subplot(2, 3, 2)
plt.imshow(fxConvolved)
plt.colorbar()
plt.title('fxConvolved')
#plt.imshow(peaks)
#plt.colorbar()
#plt.title('Location of local peaks')
#plt.show()
plt.subplot(2, 3, 3)
plt.imshow(fxConvolvedMaxFiltered)
plt.colorbar()
plt.title('fxConvolvedMaxFiltered')
plt.imshow(masses[0])
plt.subplot(2, 3, 4)
plt.imshow(peaks)
plt.colorbar()
plt.title('Local masses (first slice if 3D)')
plt.show()
plt.title('peaks')
plt.imshow(masses[0]>0.5)
plt.subplot(2, 3, 5)
plt.imshow(masses)
plt.colorbar()
plt.title('Local masses > 0.5 (first slice if 3D)')
plt.title('masses')
plt.show()
#plt.imshow(peaks*masses)
#plt.colorbar()
#plt.title('Final weighted estimation')
#plt.show()
if twoD: return peaks[0]*masses[0]
else: return peaks *masses
return masses
def tomopack(radioMM, psiMM, maxIterations=50, l=0.1, epsilon='iterative', kTrustRatio=0.75, GRAPH=False):
'''
......@@ -297,7 +289,7 @@ def tomopack(radioMM, psiMM, maxIterations=50, l=0.1, epsilon='iterative', kTrus
return numpy.real(f_x)
def indicatorFunctionToDetectorPositions(f_x, debug=False, particlesPresenceThreshold=0.9):
def indicatorFunctionToDetectorPositions(f_x, scanFixedNumber=None, massThreshold=0.5, debug=False):
'''
Takes an approximation of the indicator function f_x from tomopack and returns positions on the detector
......@@ -306,28 +298,47 @@ def indicatorFunctionToDetectorPositions(f_x, debug=False, particlesPresenceThre
f_x : 2D numpy array of floats
The value of the indicator function for each pixel on the detector
scanFixedNumber : int (optional)
The number of spheres scanned, If None peaks will be thresholded based on massThreshold
Default = None
massThreshold : float (optional)
Threshold for accepting the result of the filtering of the indicator function
Default = 0.5
debug : bool, optional
Show debug graphs (especially in the maximum filtering)
Default = False
particlesPresence : float, optional
Threshold for accepting a particle
Default = 0.9
Returns
-------
peaksPOSnJI : 2D numpy array
positions of spheres centres on the detector (rows, coloumns (JI))
'''
assert(len(f_x.shape) == 2), "detectSpheres.indicatorFunctionToDetectorPositions(): need 2D array"
f_x = _filter_maxima(f_x, debug=debug)
masses = _filter_maxima(f_x, debug=debug)
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]]
# print(f_x[pos[0],pos[1]])
for i in range(int(val)):
pPredictedIJ.append([pos[1], pos[0]])
pPredictedIJ = numpy.array(pPredictedIJ)
#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]]
## print(f_x[pos[0],pos[1]])
#for i in range(int(val)):
#pPredictedIJ.append([pos[1], pos[0]])
#pPredictedIJ = numpy.array(pPredictedIJ)
return pPredictedIJ
if scanFixedNumber:
# get the indices of all of the peaks, from highest to lowest
sortedPeakIndices = numpy.argsort(masses, axis=None)[::-1]
# get just the first scanFixedNumber of those and put them into a scanFixedNumber x 3 array
peaksPOSnJI = numpy.vstack(numpy.unravel_index(sortedPeakIndices[:scanFixedNumber], masses.shape)).T
else:
filteredPeaks = masses > massThreshold
peaksPOSnJI = numpy.argwhere(filteredPeaks)
return peaksPOSnJI
def psiSeriesScanTo3DPositions(radio, # not (necessarily) MM
......@@ -661,6 +672,7 @@ def tomopackDivergentScanTo3DPositions(radioMM, radiusMM,
)
allPeaks = fXconvolvedSeries == fXconvolvedMaximumFiltered
masses = allPeaks*fXconvolvedSeries
if verbose:
tifffile.imsave(cacheFile[:-4] + '_masses.tif', masses.astype('<f4'))
tifffile.imsave(cacheFile[:-4] + '_peaks.tif', allPeaks.astype('<f4'))
......
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