Commit 68a8108d authored by Olga Stamati's avatar Olga Stamati
Browse files

add motion kernel function and test

parent 2273a6b8
Pipeline #94329 passed with stage
in 1 minute and 49 seconds
......@@ -25,7 +25,7 @@ class TestOneSphere(unittest.TestCase):
except OSError:
pass
def test_pointToDetectorPixelRange(self):
def _test_pointToDetectorPixelRange(self):
print("\n\ntest_pointToDetectorPixelRange(): Starting up")
sourceDetectorDistMM=100
......@@ -62,7 +62,7 @@ class TestOneSphere(unittest.TestCase):
print("test_pointToDetectorPixelRange(): Done")
def test_singleSphereToDetectorPixelRange(self):
def _test_singleSphereToDetectorPixelRange(self):
# This is a very light test for reasonable behaviour, the real test will come within test_projectorsOneSphere
print("\n\ntest_singleSphereToDetectorPixelRange(): Starting up")
......@@ -83,7 +83,7 @@ class TestOneSphere(unittest.TestCase):
print("test_singleSphereToDetectorPixelRange(): Done")
def test_projectorsOneSphere(self):
def _test_projectorsOneSphere(self):
print("\n\ntest_projectorsOneSphere(): Starting up")
radiusMM = 2
......@@ -312,6 +312,18 @@ class TestOneSphere(unittest.TestCase):
print("test_projectorsOneSphere(): Done")
def test_computeMotionKernel(self):
# horizontal kernel
kernel1 = radioSphere.projectSphere.computeMotionKernel(numpy.array([0, 2, 0]), numpy.array([0, 5, 0]))
self.assertAlmostEqual(kernel1.sum(), 1, places=4)
# vertical kernel
kernel2 = radioSphere.projectSphere.computeMotionKernel(numpy.array([0, 0, 2]), numpy.array([0, 0, 5]))
self.assertAlmostEqual(kernel2.sum(), 1, places=4)
# no kernel
kernel3 = radioSphere.projectSphere.computeMotionKernel(numpy.array([0, 0, 2]), numpy.array([0, 0, 2.2]))
self.assertTrue(kernel3 is None)
#def test_projectorsThreeSpheres(self):
#print("\n\ntest_projectorsThreeSpheres(): Starting up")
#radiusMM = 2
......
......@@ -380,3 +380,53 @@ def mm2gl(radioMM, calib=None):
if calib is not None:
print("projectSphere.mm2gl() I need to be implemented")
return numpy.exp(-radioMM)
def computeMotionKernel(posPrev, posPost, displacementThreshold=1):
"""
Helper function that calculates unique motion kernel in 2D (detector plane) based on two given positions.
Its size is given by the maximum displacement and its direction corresponds to the direction of the displacement vector
Parameters
----------
posPrev : 2D numpy array of floats
Nx3 positions of the particles in MM (step i-1)
posPost : 2D numpy array of floats
Nx3 positions of the particles in MM (step i+1)
displacementThreshold : float
Threshold displacement in MM
Only if displacement is bigger than this value, a kernel is calculated
Returns
-------
kernel : 2D numpy array of floats
Motion kernel
"""
# Calculate displacement only for y-z direction (detector plane)
dx = posPost[1:]-posPrev[1:]
kernelR = None
if any(y > displacementThreshold for y in abs(dx)):
# create an empty array with size the max of the displacement
kernelSize = int(numpy.ceil(numpy.abs(dx).max()/2)*2 + 1) #round to the nearest odd number
#kernelSize = int(numpy.ceil(numpy.abs(dxSphere).max()/2)*2 + 3) #round to the nearest odd number
kernel = numpy.zeros((kernelSize, kernelSize))
# fill in a horizontal step (or hat) function
kernel[int((kernelSize - 1)/2), :] = numpy.ones(kernelSize)
# rotate kernel to follow the direction of the displacement vector
direction = dx/numpy.linalg.norm(dx)
theta = numpy.rad2deg(numpy.arctan(direction[1]/direction[0]))
kernelR = scipy.ndimage.rotate(kernel, -theta, reshape=False, order=3, mode="reflect", prefilter=False)
# normalise and clean interpolation small values (coming from the ndimage rotation)
kernelR /= kernelSize
kernelR[kernelR < 0.01] = 0
#kernelR = kernel/kernelSize
return kernelR
Supports Markdown
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