Commit 239c4152 authored by Benjy Marks's avatar Benjy Marks
Browse files

Merge branch 'master' of gricad-gitlab.univ-grenoble-alpes.fr:ttk/radioSphere

parents 0348adb5 68a8108d
......@@ -16,6 +16,7 @@ import radioSphere.DEM
import numpy
import unittest
from scipy.spatial.distance import cdist
class TestOptimiser(unittest.TestCase):
def tearDown(self):
......@@ -28,20 +29,28 @@ class TestOptimiser(unittest.TestCase):
# Just touching should not move
pos = numpy.array([[0., 0., 5.], [0., 0., -5.]])
rad = numpy.array([5., 5.])
posNew1 = radioSphere.DEM.DEM_step(pos, rad)
posNew1, _ = radioSphere.DEM.DEM_step(pos, rad)
self.assertEqual(numpy.sum(numpy.abs(pos-posNew1)), 0)
# Not touching should not move
pos = numpy.array([[0., 0., 7.], [0., 0., -7.]])
rad = numpy.array([5., 5.])
posNew2 = radioSphere.DEM.DEM_step(pos, rad)
posNew2, _ = radioSphere.DEM.DEM_step(pos, rad)
self.assertEqual(numpy.sum(numpy.abs(pos-posNew2)), 0)
# Overlapping should move
pos = numpy.array([[0., 0., 3.], [0., 0., -3.]])
rad = numpy.array([5., 5.])
posNew3 = radioSphere.DEM.DEM_step(pos, rad)
posNew3, _ = radioSphere.DEM.DEM_step(pos, rad, k=0.001)
self.assertEqual(numpy.sum(numpy.abs(pos-posNew3)) > 0, True)
delta = cdist(posNew3, posNew3) - 2*rad[0] # calculate new distances
self.assertFalse(any(delta[~numpy.eye(len(rad)).astype('bool')] < 0)) # no overlaps
# Check k influence
self.assertAlmostEqual(numpy.mean(numpy.abs(posNew3), axis=0)[2], rad[0], places=2)
posNew4, _ = radioSphere.DEM.DEM_step(pos, rad, k=0.1)
self.assertAlmostEqual(int(numpy.mean(numpy.abs(posNew4), axis=0)[2]), int(rad[0]), places=1)
if __name__ == '__main__':
unittest.main()
......@@ -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
......
import numpy as np
import numpy
from scipy.spatial.distance import cdist
def DEM_step(locMMin, radiiMM, k=0.5):
def DEM_step(posMMin, radiiMM, k=0.01):
"""
Lightweight DEM with assumption of same size spheres as mechanical regularisation
Parameters
----------
locMMin : Nx3 2D numpy array of floats
posMMin : Nx3 2D numpy array of floats
xyz positions of spheres in mm, with the origin being the middle of the detector
radiiMM : 1D numpy array of floats
......@@ -15,29 +15,36 @@ def DEM_step(locMMin, radiiMM, k=0.5):
k : float, optional
Stiffness and timestep wrapped into one
Default = 0.1
Default = 0.01
Returns
-------
locMMin : output positions
posMM : output positions
"""
locMM = locMMin.copy()
posMM = posMMin.copy()
if radiiMM.min() != radiiMM.max():
print("DEM.DEM_step(): WARINING I assume all radii are the same, taking first one")
#k = 0.1 # stiffness and timestep wrapped into one
np = len(locMM)
delta = cdist(locMM,locMM) - 2*radiiMM[0] # assuming all radii the same
# print(delta)
for i in range(0,np):
for j in range(i+1,np):
if delta[i,j] < 0:
# print(i,j,i+j)
branch_vector = locMM[i] - locMM[j]
F = -k*delta[i,j]*branch_vector
locMM[i] += F
locMM[j] -= F
return locMM
nSpheres = len(posMM)
delta = cdist(posMM, posMM) - 2*radiiMM[0] # assuming all radii the same
diag = numpy.eye(nSpheres).astype('bool')
# detect overlaps and apply DEM regularisation only for these
while any(delta[~diag] < 0):
for i in range(0, nSpheres):
for j in range(i+1, nSpheres):
if delta[i, j] < 0:
# print(i,j,i+j)
branchVector = posMM[i] - posMM[j]
F = -k*delta[i, j] * branchVector
posMM[i] += F
posMM[j] -= F
delta = cdist(posMM,posMM) - 2*radiiMM[0]
k += 0.01
return posMM, k
if __name__ == '__main__':
......
......@@ -10,7 +10,7 @@ import time
numpy.set_printoptions(suppress=True)
def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20, minDeltaMM=0.01, perturbationMM=None, sourceDetectorDistMM=100, pixelSizeMM=0.1, detectorResolution=[512,512], verbose=0, DEMcorr=False, GRAPH=False, NDDEM_output=False, blur=None, showConvergence=False, optimiseGL=False, optimiseGLcalib=None):
def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20, minDeltaMM=0.01, perturbationMM=None, sourceDetectorDistMM=100, pixelSizeMM=0.1, detectorResolution=[512,512], verbose=0, DEMcorr=False, DEMparameter=0.01, limitsX=[None, None], GRAPH=False, NDDEM_output=False, blur=None, motionKernel=None, showConvergence=False, optimiseGL=False, optimiseGLcalib=None, outDir=None):
"""
This function takes in a reference projection (radioMM) in mm and sphere positions
and moves them around to minimise the residual as far as possible.
......@@ -65,12 +65,27 @@ def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20,
Perform a soft DEM correction between all iterations?
Default = False
DEMparameter : float, optional
If DEMcorr is activated
Stiffness and timestep wrapped into one
Default = 0.01
limitsX : two-component list of floats, optional
Minimum and maximum zoom position in the direction of the X-ray beam
Default = no limits
NDDEM_output : bool, optional
Save every iteration so that it can be viewed using NDDEM - sorry Eddy
blur : float, optional
sigma of blur to pass to scipy.ndimage.gaussian_filter to
blur the radiograph at the end of everything
blur : dictionary of floats, optional
Dictionary containing a unique sigma for each sphere
to pass to `scipy.ndimage.gaussian_filter()` to blur the synthetic radiographs
Default = None
motionKernel : dictionary of floats,, optional
Dictionary containing a unique motion kernel (2d array of floats) for each sphere
to convolve the synthetic radiographs
Default = None
showConvergence : bool, optional
show a graph of the system converging with iterations
......@@ -145,6 +160,14 @@ def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20,
elif isinstance(perturbationMM, float):
perturbationMM = (perturbXscaling*perturbationMM, perturbationMM, perturbationMM)
print(f"optimiseSensitivityFields(): using a perturbation of:\n\tX: {perturbationMM[0]:0.3f}mm Y: {perturbationMM[1]:0.3f}mm Z: {perturbationMM[2]:0.3f}mm")
if not all(x is None for x in limitsX):
minZoom = max(limitsX)
maxZoom = min(limitsX)
limitsX = True
else:
limitsX = False
# initialise variables
iterations = 0
xyzMM = xyzGuessesMM.copy().astype('<f4')
......@@ -155,8 +178,8 @@ def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20,
dX = numpy.inf
# outputForFigure = []
while iterations < iterationsMax and dX > minDeltaMM:
if verbose > 0: print("\tIteration Number", iterations, end='')
if verbose > 1: print("\tperturbationMM:\t", perturbationMM)
if verbose > 0: print("\tIt:", iterations, end='')
#if verbose > 1: print("\tperturbationMM:\t", perturbationMM)
# Generate radio with current guess of position (xyzMM)
guessedRadioMM = radioSphere.projectSphere.projectSphereMM( xyzMM,
......@@ -164,11 +187,14 @@ def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20,
sourceDetectorDistMM=sourceDetectorDistMM,
pixelSizeMM=pixelSizeMM,
detectorResolution=detectorResolution,
blur=blur)
#blur=blur
)
if optimiseGL: guessedRadio = radioSphere.mm2gl(guessedRadioMM, calib=optimiseGLcalib)
else: guessedRadio = guessedRadioMM
#print(numpy.square(residualMM).sum(), end="")
# update residual
residual = radio - guessedRadio
## create an empty mask that stores locations with particles in the ROI
#if removeBackground: isBackgroundMask = numpy.ones_like(radioMM,dtype=bool)
......@@ -190,11 +216,6 @@ def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20,
#if removeBackground: isBackgroundMask[limits[sphere][0,0]:limits[sphere][0,1], limits[sphere][1,0]:limits[sphere][1,1]] = False
# update residual
residual = radio - guessedRadio
if verbose > 1: tifffile.imsave("./residuals/residual-iteration-{:04d}.tif".format(iterations), residual.astype('f4'))
for sphere in range(len(radiiMM)):
if verbose > 2: print("\tSphere {} of {}".format(sphere+1, len(radiiMM)))
# Compute local reference projection for this sphere with ROI activated around the sphere
......@@ -205,7 +226,26 @@ def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20,
detectorResolution=detectorResolution,
ROIcentreMM=xyzMM[sphere].copy(),
ROIradiusMM=radiiMM[sphere]*1.2,
blur=blur)
#blur=blur
)
# crop ROI for each sphere
radioCropSphere = radio[limits[sphere][0,0]:limits[sphere][0,1], limits[sphere][1,0]:limits[sphere][1,1]]
guessedRadioCropSphere = guessedRadio[limits[sphere][0,0]:limits[sphere][0,1], limits[sphere][1,0]:limits[sphere][1,1]]
# apply kernel to each ROI
if motionKernel:
if motionKernel[sphere] is not None:
guessedRadioCropSphere = scipy.ndimage.convolve(guessedRadioCropSphere, motionKernel[sphere])
sphereRefProjectionMM = scipy.ndimage.convolve(sphereRefProjectionMM, motionKernel[sphere])
if blur:
if blur[sphere] is not None:
guessedRadioCropSphere = scipy.ndimage.gaussian_filter(guessedRadioCropSphere, sigma=blur[sphere])
sphereRefProjectionMM = scipy.ndimage.gaussian_filter(sphereRefProjectionMM, sigma=blur[sphere])
residualCropSphere = radioCropSphere - guessedRadioCropSphere
if optimiseGL: sphereRefProjection = radioSphere.mm2gl(sphereRefProjectionMM, calib=optimiseGLcalib)
else: sphereRefProjection = sphereRefProjectionMM
......@@ -222,7 +262,7 @@ def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20,
plt.imshow(residual, cmap='coolwarm', vmin=radioSphere.mm2gl(-0.1), vmax=radioSphere.mm2gl(0.1))
else:
plt.title(f"Current Residual (mm) LUT: [-0.1, 0.1] iteration={iterations}")
plt.imshow(residual, cmap='coolwarm', vmin=-0.1, vmax=0.1)
plt.imshow(residual, cmap='coolwarm', vmin=-radiiMM[sphere], vmax=radiiMM[sphere])
plt.pause(0.01)
......@@ -241,12 +281,23 @@ def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20,
detectorResolution=detectorResolution,
ROIcentreMM=xyzMM[sphere].copy(),
ROIradiusMM=radiiMM[sphere]*1.2,
blur=blur)
#blur=blur)
)
# apply kernel to each ROI
if motionKernel:
if motionKernel[sphere] is not None:
perturbedProjectionMM = scipy.ndimage.convolve(perturbedProjectionMM, motionKernel[sphere])
if blur:
if blur[sphere] is not None:
perturbedProjectionMM = scipy.ndimage.gaussian_filter(perturbedProjectionMM, sigma=blur[sphere])
if optimiseGL: perturbedProjection = radioSphere.mm2gl(perturbedProjectionMM, calib=optimiseGLcalib)
else: perturbedProjection = perturbedProjectionMM
sensXYZ[i] = sphereRefProjection - perturbedProjection
if GRAPH:
plt.subplot(1,4,2+i)
if optimiseGL:
......@@ -266,16 +317,21 @@ def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20,
LSQret = scipy.optimize.least_squares(optimiseMe,
[1.0, 1.0, 1.0],
args=[residual[limits[sphere][0,0]:limits[sphere][0,1], limits[sphere][1,0]:limits[sphere][1,1]], sensXYZ],
args=[residualCropSphere, sensXYZ],
verbose=False,
method='lm',
diff_step=1.0)
diff_step=1)
if LSQret['success'] == True:
#print('LSQ Success!')
#print(LSQret['x'])
#step = LSQret['x'] * perturbationMM[i]
#for i in range(3):
# xyzMM[sphere][i] -= step[i] * perturbationMM[i]
xyzMM[sphere] -= LSQret['x'] * perturbationMM
# check if current position falls inside the limits, and if not bring it back
if limitsX:
if xyzMM[sphere, 0] > minZoom: xyzMM[sphere, 0] = minZoom
if xyzMM[sphere, 0] < maxZoom: xyzMM[sphere, 0] = maxZoom
else:
print("LSQ failed to converge")
......@@ -285,17 +341,30 @@ def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20,
### End optimisation iterations
### Now a soothing DEM step can be applied to the current updated XYZ positions
#### Now a soothing DEM step can be applied to the current updated XYZ positions
#if DEMcorr:
#import radioSphere.DEM
#xyzMMnew, k = radioSphere.DEM.DEM_step(xyzMM, radiiMM, k=0.1)
##if verbose > 0: print(" DEM changed positions by: ", numpy.linalg.norm(xyzMM-xyzMMnew), end='')
#xyzMM = xyzMMnew
if DEMcorr:
import radioSphere.DEM
xyzMMnew = radioSphere.DEM.DEM_step(xyzMM,radiiMM)
if verbose > 0: print(" DEM changed positions by: ", numpy.linalg.norm(xyzMM-xyzMMnew), end='')
xyzMM = xyzMMnew
from scipy.spatial.distance import cdist
nSpheres = len(xyzMM)
delta = cdist(xyzMM, xyzMM) - 2*radiiMM[0]
diag = numpy.eye(nSpheres).astype('bool')
# check for overlaps and if yes, apply DEM regularisation
if any(delta[~diag] < 0):
import radioSphere.DEM
xyzMMDEM, k = radioSphere.DEM.DEM_step(xyzMM, radiiMM, k=DEMparameter)
if verbose > 1: print(" DEM changed positions by: ", numpy.linalg.norm(xyzMM-xyzMMDEM), end='')
xyzMM = xyzMMDEM
dX = numpy.linalg.norm(xyzMM-xyzMMprev)
# outputForFigure.append([numpy.linalg.norm(step),numpy.sqrt(numpy.sum(residualMM.flatten()**2))])
# if verbose > 0: print(" |deltaMM|: ", numpy.linalg.norm(xyzMM-xyzMMprev), end='')
print(f" |deltaMM|: {dX:0.5f}", end='\r')
print(f" LSQ: {LSQret['x']} | |dMM|: {dX:0.5f}", end='\r')
xyzMMprev = xyzMM.copy()
......@@ -312,7 +381,7 @@ def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20,
plt.ylabel('Sum of squared residuals (GL^2)')
else:
plt.ylabel('Sum of squared residuals (mm^2)')
plt.pause(1e-2)
plt.pause(2)
plt.show()
iterations += 1
......@@ -328,9 +397,8 @@ def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20,
print("\n\toptimiseSensitivityFields(): hit max iterations ({})".format(iterations))
plt.ioff()
return xyzMM
return xyzMM
def _optimiseSensitivityFieldsMultiProj(radioMM, xyzGuessesMM, radiiMM, transformationCentresMM, transformationMatrices, iterationsMax=20, minDeltaMM=0.01, perturbationMM=numpy.array([0.5, 0.5, 0.5]), sourceDetectorDistMM=100, pixelSizeMM=0.1, detectorResolution=[512,512], solver='leastsquares', verbose=0, DEMcorr=False, GRAPH=False):
"""
......
......@@ -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