Skip to content

GitLab

Commit 8932b939 authored by Olga Stamati's avatar Olga Stamati
Browse files

big update in detectSpheres module and tests

parent b8f38577
Pipeline #82642 passed with stage
in 1 minute and 48 seconds
Expand all Hide whitespace changes
Inline Side-by-side
tests/test_detection.py
View file @ 8932b939
......@@ -11,14 +11,61 @@ import radioSphere
import numpy
import unittest
import matplotlib.pyplot as plt
import shutil
class TestDetection(unittest.TestCase):
def tearDown(self):
try:
shutil.rmtree("./output")
pass
except OSError:
pass
def test_fxToPeaks(self):
# 2D fx
fx = numpy.zeros((200, 100))
fx[40, 70] = 0.5
peaks = radioSphere.detectSpheres.indicatorFunctionToDetectorPositions(fx, scanFixedNumber=1)
self.assertTrue(numpy.allclose([0, 40, 70], peaks))
# 3D fx
fx = numpy.zeros((30, 200, 100))
fx[5, 40, 70] = 0.5
peaksa = radioSphere.detectSpheres.indicatorFunctionToDetectorPositions(fx, scanFixedNumber=1)
peaksb = radioSphere.detectSpheres.indicatorFunctionToDetectorPositions(fx, massThreshold=0.4)
self.assertTrue(numpy.allclose([5, 40, 70], peaksa))
self.assertTrue(numpy.allclose([5, 40, 70], peaksb))
def test_peaksTo3DPos(self):
# peaks at middle of detector so real pos should be at 0 of coord system
peaks = numpy.array([[2, 100, 50]])
CORxPos = numpy.arange(3)
posXYZ = radioSphere.detectSpheres.detectorPeaksTo3DCoordinates(peaks, CORxPos, detectorResolution=[200, 100])
self.assertTrue(numpy.allclose(posXYZ, [2, 0, 0]))
# backwards compatibility with projectSphere
detectorResolution=[200, 100]
sourceDetectorDistMM = 50
pixelSizeMM = 0.5
peaks = [2, 20, 80]
CORxPos = numpy.arange(3)
posXYZ = radioSphere.detectSpheres.detectorPeaksTo3DCoordinates(numpy.array([peaks]), CORxPos,
pixelSizeMM=pixelSizeMM,
detectorResolution=detectorResolution,
sourceDetectorDistMM=sourceDetectorDistMM)
p = radioSphere.projectSphere.projectSphereMM(posXYZ,
numpy.ones(posXYZ.shape[0]),
detectorResolution=[200, 100],
sourceDetectorDistMM=sourceDetectorDistMM,
pixelSizeMM=pixelSizeMM)
sphereCentreDetector = numpy.argwhere(p == p.max())[0]
# 1px tol
self.assertTrue(numpy.allclose(peaks[1:], sphereCentreDetector, atol=1))
def test_tomopackOneSphere(self):
COR = 50
r = 1
......@@ -35,8 +82,8 @@ class TestDetection(unittest.TestCase):
detectorResolution = detectorResolution,
sourceDetectorDistMM = numpy.inf)
f_x1 = radioSphere.detectSpheres.tomopack(p1, psi)
peaksJI1 = radioSphere.detectSpheres.indicatorFunctionToDetectorPositions(f_x1, scanFixedNumber=1)
fx1 = radioSphere.detectSpheres.tomopack(p1, psi)
peaksJI1 = radioSphere.detectSpheres.indicatorFunctionToDetectorPositions(fx1, scanFixedNumber=1)[:, 1:]
self.assertEqual(peaksJI1[0][0], detectorResolution[0]//2)
self.assertEqual(peaksJI1[0][1], detectorResolution[1]//2)
......@@ -47,17 +94,17 @@ class TestDetection(unittest.TestCase):
detectorResolution = detectorResolution,
sourceDetectorDistMM = numpy.inf)
f_x2 = radioSphere.detectSpheres.tomopack(p2, psi)
peaksJI2 = radioSphere.detectSpheres.indicatorFunctionToDetectorPositions(f_x2, massThreshold=0.1)
fx2 = radioSphere.detectSpheres.tomopack(p2, psi)
peaksJI2 = radioSphere.detectSpheres.indicatorFunctionToDetectorPositions(fx2, massThreshold=0.1)[:, 1:]
spherePosDetector2 = numpy.where(p2 == p2.max())
estimatedSpherePosDetector2 = numpy.where(f_x2 == f_x2.max())
estimatedSpherePosDetector2 = numpy.where(fx2 == fx2.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
# Case 3: A bit of zoom, 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,
......@@ -68,27 +115,27 @@ class TestDetection(unittest.TestCase):
detectorResolution = detectorResolution,
sourceDetectorDistMM = 100)
f_x3 = radioSphere.detectSpheres.tomopack(p3, psi3)
peaksJI3 = radioSphere.detectSpheres.indicatorFunctionToDetectorPositions(f_x3, scanFixedNumber=1)
fx3 = radioSphere.detectSpheres.tomopack(p3, psi3)
peaksJI3 = radioSphere.detectSpheres.indicatorFunctionToDetectorPositions(fx3, scanFixedNumber=1)[:, 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
# Case 4: A bit of zoom, 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)
fx4 = radioSphere.detectSpheres.tomopack(p4, psi3)
peaksJI4 = radioSphere.detectSpheres.indicatorFunctionToDetectorPositions(fx4, scanFixedNumber=1)[:, 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
# Tomopack 10 random spheres
def test_tomopack10Spheres(self):
COR = 50
r = 1
......@@ -118,25 +165,94 @@ class TestDetection(unittest.TestCase):
sourceDetectorDistMM = sourceDetectorDistMM)
f_x = radioSphere.detectSpheres.tomopack(p, psi)
peaksJI = radioSphere.detectSpheres.indicatorFunctionToDetectorPositions(f_x, scanFixedNumber=pos.shape[0])
fx = radioSphere.detectSpheres.tomopack(p, psi)
posTomopackXYZ = numpy.zeros([peaksJI.shape[0], 3])
posTomopackXYZ[:, 0] += COR
# call + test helper which transforms fx guess to detector peaks
peaksJI = radioSphere.detectSpheres.indicatorFunctionToDetectorPositions(fx, scanFixedNumber=pos.shape[0])
# 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
# call + test helper which transforms det peaks to 3d coord
posTomopackXYZ = radioSphere.detectSpheres.detectorPeaksTo3DCoordinates(peaksJI,
numpy.ones(peaksJI.shape[0])*COR,
pixelSizeMM=pixelSizeMM,
detectorResolution=p.shape,
sourceDetectorDistMM=sourceDetectorDistMM
)
# sort tomopack guess to calculate residual directly on positions
posTomopackXYZ = posTomopackXYZ[numpy.argsort(posTomopackXYZ[:, 1])]
self.assertTrue(numpy.linalg.norm(pos - posTomopackXYZ) < r)
self.assertTrue(numpy.mean(numpy.linalg.norm(pos - posTomopackXYZ, axis=1) < 0.5*r))
def test_tomopackDivergentBeamOneSphere(self):
r = 1
COR = 50
detectorResolution = [400, 600]
# Case 1: simplest case with 1 sphere in the middle of the detector
pos1 = numpy.array([[COR, 0, 0]])
p1 = radioSphere.projectSphere.projectSphereMM(pos1,
numpy.array([r]),
detectorResolution = detectorResolution,
sourceDetectorDistMM = 100)
posXYZ1a = radioSphere.detectSpheres.tomopackDivergentScanTo3DPositions(p1, r, CORxRef=COR, CORxNumber=5, scanFixedNumber=1, saveSeries=False, verbose=False)
posXYZ1b = radioSphere.detectSpheres.tomopackDivergentScanTo3DPositions(p1, r, CORxRef=COR, CORxNumber=5, saveSeries=False, massThreshold=0.2)
self.assertTrue(numpy.allclose(pos1, posXYZ1a[0], atol=0.1))
self.assertTrue(numpy.allclose(pos1, posXYZ1b[0], atol=0.1))
# Case 2: 1 sphere shifted from COR and centre of detector
pos2 = numpy.array([[numpy.random.uniform(46, 48), numpy.random.uniform(5, -5), numpy.random.uniform(-3, 3)]])
#pos2 = numpy.array([[47.20396842, 0, 0]])
#pos2 = numpy.array([[46.22945869, -1.63507697, 0.96052506]])
p2 = radioSphere.projectSphere.projectSphereMM(pos2,
numpy.array([r]),
detectorResolution= detectorResolution,
sourceDetectorDistMM=100)
posXYZ2a = radioSphere.detectSpheres.tomopackDivergentScanTo3DPositions(p2, r, CORxMin=COR-4*r, CORxMax=COR+4*r, CORxNumber=20, scanFixedNumber=1, saveSeries=False, verbose=False)
posXYZ2b = radioSphere.detectSpheres.tomopackDivergentScanTo3DPositions(p2, r, CORxMin=COR-4*r, CORxMax=COR+4*r, CORxNumber=20, saveSeries=False, verbose=False, massThreshold=0.1)
self.assertTrue(numpy.allclose(pos2, posXYZ2a[0], atol=0.5))
self.assertTrue(numpy.allclose(pos2, posXYZ2b[0], atol=0.5))
def test_tomopackDivergentBeam10Spheres(self):
r = 1
COR = 50
detectorResolution = [400, 600]
pixelSizeMM = 0.1
xMin, xMax = 40, 60
yMin, yMax = -5, 5
zMin, zMax = -7, 7
pos = numpy.zeros((10, 3))
pos[:, 0] = numpy.arange(xMin, xMax, 2*r)
pos[:, 1] = yMin + r + (yMax - yMin - 2*r)*numpy.random.rand(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,
numpy.ones_like(pos[:, 0])*r,
pixelSizeMM=pixelSizeMM,
detectorResolution= detectorResolution,
sourceDetectorDistMM=100)
posXYZa = radioSphere.detectSpheres.tomopackDivergentScanTo3DPositions(p, r,
CORxMin=COR-11*r, CORxMax=COR+11*r, CORxNumber=50,
scanFixedNumber=pos.shape[0],
saveSeries=True, saveSeriesDirectory="./output")
posXYZb = radioSphere.detectSpheres.tomopackDivergentScanTo3DPositions(p, r, CORxMin=COR-11*r, CORxMax=COR+11*r, CORxNumber=50,
fXseries = "./output/fXseries.tif", psiXseries = "./output/psiXseries.tif",
massThreshold=0.1,
saveSeries=False)
# sort guess based on X axis to calculate residual directly on position
posXYZa = posXYZa[numpy.argsort(posXYZa[:, 0])]
posXYZb = posXYZb[numpy.argsort(posXYZb[:, 0])]
self.assertTrue(numpy.mean(numpy.linalg.norm(pos - posXYZa, axis=1)) < 1.5*r)
self.assertTrue(numpy.mean(numpy.linalg.norm(pos - posXYZb, axis=1)) < 1.5*r)
if __name__ == "__main__":
unittest.main()
tools/detectSpheres.py
View file @ 8932b939
This diff is collapsed.
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