Commit 2273a6b8 authored by Olga Stamati's avatar Olga Stamati
Browse files

optimisation: add motion kernel for each sphere

parent 104f2d5d
Pipeline #94314 passed with stage
in 1 minute and 47 seconds
......@@ -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, DEMparameter=0.01, limitsX=[None, None], GRAPH=False, NDDEM_output=False, blur=False, showConvergence=False, optimiseGL=False, optimiseGLcalib=None, outDir=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.
......@@ -77,9 +77,15 @@ def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20,
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
......@@ -181,7 +187,7 @@ 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)
......@@ -208,7 +214,6 @@ def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20,
pixelSizeMM=pixelSizeMM,
detectorResolution=detectorResolution)
#if removeBackground: isBackgroundMask[limits[sphere][0,0]:limits[sphere][0,1], limits[sphere][1,0]:limits[sphere][1,1]] = False
for sphere in range(len(radiiMM)):
......@@ -221,9 +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
......@@ -259,7 +281,17 @@ 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
......@@ -285,7 +317,7 @@ 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)
......
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