Commit 104f2d5d authored by Olga Stamati's avatar Olga Stamati
Browse files

optimisation: add limits in the x-ray direction

parent 5e77c079
Pipeline #94306 passed with stage
in 1 minute and 46 seconds
......@@ -7,7 +7,7 @@ def DEM_step(posMMin, radiiMM, k=0.01):
Parameters
----------
posMM : 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
......
......@@ -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, 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=False, 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.
......@@ -70,6 +70,10 @@ def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20,
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
......@@ -151,6 +155,12 @@ def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20,
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
......@@ -286,6 +296,10 @@ def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20,
#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")
......@@ -310,10 +324,8 @@ def optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20,
# check for overlaps and if yes, apply DEM regularisation
if any(delta[~diag] < 0):
print("\nFound overlaps")
import radioSphere.DEM
xyzMMDEM, k = radioSphere.DEM.DEM_step(xyzMM, radiiMM, k=DEMparameter)
print("k", k)
if verbose > 1: print(" DEM changed positions by: ", numpy.linalg.norm(xyzMM-xyzMMDEM), end='')
xyzMM = xyzMMDEM
......
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