Welcome to radioSphere’s documentation!

Sphere detection

Library of RadioSphere spheres detection functions.

detectSpheres.computeFxAndPsiSeries(radioMM, radiusMM, CORxPositions, pixelSizeMM=0.1, sourceDetectorDistMM=100, blur=0.0, maxIterations=50, l=0.5, epsilon='iterative', kTrustRatio=0.75, nProcesses=8)[source]

This function takes in a single divergent projection and positions along the ray direction, and returns the corresponding indicator function (fx) and structuring element (psi) series by running (in parallel) a series of tomopack

Parameters:

  • radioMM (2D numpy array of floats) – Radiography containing “distance” projections in mm This is typically the result of mu*log(I/I_0)

  • radiusMM (float) – Particle radius in mm

  • CORxPositions (1D numpy array of floats) – Positions in the x-ray direction on which a series of tomopack will be run

  • pixelSizeMM (float, optional) – Pixel size on detector in mm Default = 0.1

  • sourceDetectorDistMM (float, optional) – Distance between x-ray source and middle of detector. Set as numpy.inf for parallel projection Default = 100

  • blur (float, optional) – Sigma of blur to pass to scipy.ndimage.gaussian_filter to blur the synthetic psi Default = 0

  • maxIterations (int, optional) – Number of iterations to run tomopack Default = 50

  • l (float, optional) – “Lambda” parameter which controls relaxation in tomopack iterations Default = 0.5

  • epsilon (float, optional, or 'iterative') – Trust cutoff wavelengths in psi If the radiograph is in mm so is this parameter Default = ‘iterative’ – updating until kTrustRatio (below) is reached

  • kTrustRatio (float, optional) – Ratio of cutoff psi wavelengths Default = 0.75

  • nProcesses (integer, optional) – Number of processes for multiprocessing Default = number of CPUs in the system

Returns:

  • fXseries (3D numpy array of floats) – Indicator function series along the x-ray direction Its shape is (len(CORxPositions), radioMM.shape[0], radioMM.shape[1])

  • psiXseries (3D numpy array of floats) – Structuring element series along the x-ray direction Its shape is (len(CORxPositions), radioMM.shape[0], radioMM.shape[1])

detectSpheres.detectorPeaksTo3DCoordinates(peaksJI, CORxPos, detectorResolution=[512, 512], pixelSizeMM=0.1, sourceDetectorDistMM=100)[source]

This function takes positions on the detector (identified peaks of the indicator function) and a set of positions in the x-ray direction (where a series of tomopack was run) and returns the corresponding XYZ coordinates

Parameters:

  • peaksJI (2D numpy array of int) – Positions of spheres centres on the detector (rows, coloumns (JI))

  • CORxPos (1D numpy array of floats) – Positions in the x-ray direction on which a series of tomopack was run If not a divergent beam, this should be an array full of the same number – COR

  • pixelSizeMM (float, optional) – Pixel size on detector in mm Default = 0.1

  • detectorResolution (2-component list of ints, optional) – Number of pixels rows, columns of detector Default = [512,512]

  • sourceDetectorDistMM (float, optional) – Distance between x-ray source and middle of detector. Set as numpy.inf for parallel projection Default = 100

Returns:

posXYZmm – Positions of spheres centres (N) in 3D (XYZ) in mm If not a divergent beam, posX is the input COR

Return type:

(N x 3) 2D numpy array

detectSpheres.indicatorFunctionToDetectorPositions(fx, scanFixedNumber=None, massThreshold=0.5, debug=False)[source]

This function takes an indicator function fx and returns peaks on the detector. fx can be a 2D array, which is directly the raw output from tomopack, or a 3D array, if a series of tomopack was run.

In both cases, the peaks returned as a (n x 3) 2D array, where the first column corresponds to the slice of the fx series that peaks were found. If a 2D input fx is passed, the first column is always 0 (zero slice in 3D).

Parameters:

  • fx (2D or 3D numpy array of floats) – Approximation of indicator function For a divergent beam this should be a 3D array

  • scanFixedNumber (int, optional) – The number of spheres scanned If None, peaks will be thresholded based on massThreshold (below) Default = None

  • massThreshold (float, optional) – Threshold for accepting the result of the filtered indicator function Activated if scanFixedNumber (above) is None Default = 0.5

  • debug (bool, optional) – Show debug graphs (especially in the maximum filtering) Default = False

Returns:

peaksPOSnJI – Positions of spheres centres (N) on the detector (slice, rows, coloumns (JI)) Where slice is the CORx slice the resonance peak falls in (only makes sense for a 3D fx)

Return type:

(N x 3) 2D numpy array

detectSpheres.tomopack(radioMM, psiMM, maxIterations=50, l=0.5, epsilon='iterative', kTrustRatio=0.75, verbose=False, graphShow=False)[source]

‘tomopack’ FFT-based algorithm for sphere identification by Stéphane Roux. the FFT of psiMM as a structuring element to pick up projected spheres in radioMM using a parallel projection hypothesis.

Parameters:

  • radioMM (2D numpy array of floats) – Radiography containing “distance” projections in mm This is typically the result of mu*log(I/I_0)

  • psiMM (2D numpy array of floats) – A (synthetic) projection in mm of a single sphere in the centre of the detector This is the structuring element used for FFT recognition. Should be same size as radioMM

  • maxIterations (int, optional) – Number of iterations to run the detection for Default = 50

  • l (float, ptional) – “Lambda” parameter which controls relaxation in the iterations Default = 0.5

  • epsilon (float, optional, or 'iterative') – Trust cutoff wavelengths in psi. If the radiograph is in mm so is this parameter Default = ‘iterative’ – updating until kTrustRatio (below) is reached

  • kTrustRatio (float, optional) – Ratio of cutoff psi wavelengths Default = 0.75

  • verbose (bool, optional) – Get to know what the function is really thinking Default = False

  • graphShow (bool, optional) – Show evolution of fx during iterations, and final one together with radio and psi Default = False

Returns:

fx – The approximation of the indicator function Same size as radioMM, with high values where we think spheres centres are Consider using indicatorFunctionToDetectorPositions() to get the IJ detector coordinates of the centres

Return type:

2D numpy array

detectSpheres.tomopackDivergentScanTo3DPositions(radioMM, radiusMM, CORxRef=None, CORxMin=None, CORxMax=None, CORxNumber=100, fXseries=None, psiXseries=None, pixelSizeMM=0.1, sourceDetectorDistMM=100, blur=0.0, maxIterations=50, l=0.5, epsilon='iterative', kTrustRatio=0.75, scanFixedNumber=None, massThreshold=0.5, saveSeries=True, saveSeriesDirectory=None, nProcesses=8, verbose=True)[source]

This function takes in a single divergent projection, and will run tomopack generating different psis by varying their position as Centre Of Rotation (COR) in the x-direction, from CORxMin to CORxMax in CORxNumber steps.

The resulting series of indicator functions is analysed and a 3D position guess for all identified spheres is returned.

Parameters:

  • radioMM (2D numpy array of floats) – Radiography containing “distance” projections in mm This is typically the result of mu*log(I/I_0)

  • radiusMM (float) – Particle radius in mm

  • fXseries (3D numpy array of floats, optional) – Indicator function series along the x-ray direction Default = None

  • psiXseries (3D numpy array of floats, optional) – Structuring element series along the x-ray direction Default = None

  • pixelSizeMM (float, optional) – Pixel size on detector in mm Default = 0.1

  • sourceDetectorDistMM (float, optional) – Distance between x-ray source and middle of detector Default = 100

  • blur (float, optional) – Sigma of blur to pass to scipy.ndimage.gaussian_filter to blur the synthetic psi Default = 0

  • maxIterations (int, optional) – Number of iterations to run tomopack Default = 50

  • l (float, optional) – “Lambda” parameter which controls relaxation in tomopack iterations Default = 0.5

  • epsilon (float, optional, or 'iterative') – Trust cutoff wavelengths in psi for tomopack If the radiograph is in mm so is this parameter Default = ‘iterative’ – updating until kTrustRatio (below) is reached

  • kTrustRatio (float, optional) – Ratio of cutoff psi wavelengths for tomopack Default = 0.75

  • scanFixedNumber (int, optional) – The number of spheres scanned If None peaks will be thresholded based on massThreshold (below) Default = None

  • massThreshold (float, optional) – Threshold for accepting the result of the filtered indicator function Activated if scanFixedNumber (above) is None Default = 0.5

  • saveSeries (boolean, optional) – Save indicator function, psi and filtered series as tif files Default = True

  • saveSeriesDirectory (str, optional) – Directory to save the computed series Default = None

  • nProcesses (integer, optional) – Number of processes for multiprocessing Default = number of CPUs in the system

  • verbose (boolean, optional) – Print what the function is doing? Default = True

Returns:

positionsXYZmm – Positions of spheres centres in 3D (in mm)

Return type:

2D numpy array

Sphere projection

projectSphere.computeLinearBackground(radioMM, mask=None)[source]

This function computes a plane-fit for background greylevels to help with the correction of the background

Parameters:

  • radioMM (2D numpy array of floats) – The image

  • mask (2D numpy array of bools, optional) – Masked zone to fit?

Returns:

background – Same size as radioMM

Return type:

2D numpy array of floats

projectSphere.computeMotionKernel(posPrev, posPost, displacementThreshold=1)[source]

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

posPost2D numpy array of floats

Nx3 positions of the particles in MM (step i+1)

displacementThresholdfloat

Threshold displacement in MM Only if displacement is bigger than this value, a kernel is calculated

Returns:

kernel – Motion kernel

Return type:

2D numpy array of floats

projectSphere.gl2mm(radio, calib=None)[source]

This function takes a greylevel radiograph (I/I0) and returns a radiograph in mm (L), representing the path length encountered.

Parameters:

  • radio (a 2D numpy array of floats) –

  • calib (dictionary (optional)) – This contains a calibration of greylevels to mm If not passed mu is assumed to be 1.

projectSphere.mm2gl(radioMM, calib=None)[source]

This function takes a radiograph in mm (L) and returns a radiograph in greylevels (I/I0)

Parameters:

  • radioMM (a 2D numpy array of floats) –

  • calib (dictionary (optional)) – This contains a calibration of greylevels to mm If not passed mu is assumed to be 1.

projectSphere.pointToDetectorPixelRange(posMM, sourceDetectorDistMM=100, pixelSizeMM=0.1, detectorResolution=[512, 512])[source]

This function gives the detector pixel that a single point will affect. The main idea is that it will be used with ROIaroundSphere option for projectSphereMM() and in turn singleSphereToDetectorPixelRange() in order to only project the needed pixels.

Parameters:

  • posMM (1x3 2D numpy array of floats) – xyz position of sphere in mm, with the origin being the middle of the detector

  • sourceDetectorDistMM (float, optional) – Distance between x-ray source and middle of detector Set as numpy.inf for parallel projection Default = 100

  • pixelSizeMM (float, optional) – Pixel size on detector in mm Default = 0.1

  • detectorResolution (2-component list of ints, optional) – Number of pixels rows, columns of detector Default = [512,512]

Returns:

detectorPixel – row, column (j,i) coordinate on detector as per figures/projectedCoords_v2.pdf

Return type:

tuple

projectSphere.projectSphereMM(spheresPositionMM, radiiMM, ROIcentreMM=None, ROIradiusMM=None, sourceDetectorDistMM=100, pixelSizeMM=0.1, detectorResolution=[512, 512], projector='C', transformationCentreMM=None, transformationMatrix=None, blur=None, displacementsMM=None)[source]

This is the python wrapping function for the C++ projector, it gets projection geometry, list of particle positions and radii, projected in the X direction. The output is the crossed distance for each sphere in mm.

Please refer to the figures/projectedCoords_v2 for geometry

In order to allow projections from diffferent angles, an XYZ centre and a transformation matrix can be provided, which will be applied to the particle positions.

Parameters:

  • spheresPositionMM (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) – Particle radii for projection

  • ROIcentreMM (3-component vector of floats, optional) – Particle position for ROI Default = Disactivated (None)

  • ROIradiusMM (float, optional) – Particle radius for ROI Default = Disactivated (None)

  • sourceDetectorDistMM (float, optional) – Distance between x-ray source and middle of detector. Set as numpy.inf for parallel projection Default = 100

  • pixelSizeMM (float, optional) – Pixel size on detector in mm Default = 0.1

  • detectorResolution (2-component list of ints, optional) – Number of pixels rows, columns of detector Default = [512,512]

  • projector (string, optional) – Algorithm for the projector (leave this alone for now) Default = ‘C’

  • transformationCentreMM (3-component vector) – XYZ centre for a transformation Default = None

  • transformationMatrix (3x3 numpy array) – XYZ transformation matrix to apply to coordinates Default = None

  • blur (float, optional) – sigma of blur to pass to scipy.ndimage.gaussian_filter to blur the radiograph at the end of everything

Returns:

projectionMM – Radiography containing the total crossed distance through the spheres distance in mm for each beam path. To turn this into a radiography, the distances should be put into a calibrated Beer-Lambert law

Return type:

2D numpy array of floats

projectSphere.singleSphereToDetectorPixelRange(spherePositionMM, radiusMM, radiusMargin=0.1, sourceDetectorDistMM=100, pixelSizeMM=0.1, detectorResolution=[512, 512], transformationCentreMM=None, transformationMatrix=None)[source]

This function gives the detector pixel range that a single sphere will affect. The main idea is that it will be used with ROIaroundSphere option for projectSphereMM() in order to only project the needed pixels.

Parameters:

  • spherePositionMM (1x3 2D numpy array of floats) – xyz position of sphere in mm, with the origin being the middle of the detector

  • radiusMM (float) – Particle radius for projection

  • radiusMargin (float) – Multiplicative margin on radius Default = 0.1

  • ROIaroundSphere (bool, optional) – If there is only one sphere, only compute a region-of-interest radiography? Default = False

  • sourceDetectorDistMM (float, optional) – Distance between x-ray source and middle of detector Set as numpy.inf for parallel projection Default = 100

  • pixelSizeMM (float, optional) – Pixel size on detector in mm Default = 0.1

  • detectorResolution (2-component list of ints, optional) – Number of pixels rows, columns of detector Default = [512,512]

  • transformationCentreMM (3-component vector) – XYZ centre for a transformation Default = None

  • transformationMatrix (3x3 numpy array) – XYZ transformation matrix to apply to coordinates Default = None

Returns:

JIrange

Return type:

range in rows, colums of detector concerned by this grain

Position optimisation

optimisePositions.optimiseSensitivityFields(radioMM, xyzGuessesMM, radiiMM, iterationsMax=20, minDeltaMM=0.01, perturbationMM=None, displacementsMM=None, sourceDetectorDistMM=100, pixelSizeMM=0.1, detectorResolution=[512, 512], verbose=0, DEMcorr=False, DEMparameter=0.01, limitsX=None, regularisation=False, regularisationThreshold=inf, weightX=1.0, GRAPH=False, NDDEM_output=False, blur=None, motionKernel=None, showConvergence=False, optimiseGL=False, optimiseGLcalib=None, outDir=None)[source]

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. The technique used is the calculation of a sensitivity field in x, y, z with a perturbation of the synthetic radiograph. This is done locally in a ROI around each particle, which avoids doing the projection for all pixels.

Parameters:

  • radioMM (2D numpy array of floats) – Projection of distance crossed in spheres in MM

  • xyzGuessesMM (2D numpy array of floats) – Nx3 guesses for the positions of the particles in MM

  • radiiMM (1D numpy array of floats) – Particle radii for projection

  • iterationsMax (int, optional) – Number of iterations to perform Default = 20

  • minDeltaMM (float, optional) – If step size gets below this value, stop Default = 0.01

  • perturbationMM (either a single float or a 3-component list of floats, optional) – Size of the perturbation step to compute the sensitivity field Default = mean particle radius * numpy.array([SOD^2/(2*SDD), 0.5, 0.5])

  • displacementsMM (2D numpy array of floats, optional) – Nx3 displacement vectors to account for blur Default = None

  • sourceDetectorDistMM (float, optional) – Distance between x-ray source and middle of detector Default = 100

  • pixelSizeMM (float, optional) – Pixel size on detector in mm Default = 0.1

  • detectorResolution (2-component list of ints, optional) – Number of pixels height, width of detector Default = [512,512]

  • verbose (int, optional) – Print noisly? 0 = quiet except warning, 1=iterations, 2=everything Default = 1

  • DEMcorr (bool, optional) – 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 = None, no limits

  • regularisation (bool, optional) – Perform a harsh “regularisation” in time in the direction of the X-ray beam Default = False

  • regularisationThreshold (float, optional) – Threshold (MM) in distance between current X-position and input guess If this is exceeded, the X-position is forced back to the input one Default = numpy.inf

  • weightX (float, optional) – Extra weighting factor applied in the direction of the X-ray beam It is applied inside the least_squares minimisation damping down X oscillations Default = 1.

  • NDDEM_output (bool, optional) – Save every iteration so that it can be viewed using NDDEM - sorry Eddy

  • 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 Default = False

  • optimiseGL (bool, optional) – Perform optimisation on image in greylevels rather than mm? Default = False

  • optimiseGLcalib (dictionary, optional) – Calibration to apply for converting mm to greylevels and back Default = None

Returns:

xzyMM

Return type:

Corrected positions array

Attenuation calibration

calibrateAttenuation.generateFitParameters(calibSphereLog, pixelSizeMM, sourceDetectorDistMM, sourceObjectDistanceMM, radiusMM, centreYXpx, fitFunction=<function linearFit>, outputPath=False, verbose=False)[source]

Fit an attenuation law to a logged normalised radiograph of a single particle.

Parameters:

  • calibSphereLog (2D numpy array of floats) – A radiograph of a single particle divided by the background intensity and then logged.

  • pixelSizeMM (float) – The size of a single pixel in mm on the detector panel.

  • sourceDetectorDistMM (float) – The source-detector distance in mm.

  • sourceObjectDistanceMM (float) – The source-object distance in mm.

  • radiusMM (float) – The particle radius in mm.

  • centreYXpx (1D numpy array of floats) – The y and x location of the centre of the particle in pixels.

  • fitFunction (function handle (optional)) – The fitting function to use. By default fit a linear fitting law. Options are cubicFit, quadraticFit, linearFit and linearFitNoIntercept

  • outputPath (string (optional)) – A path to save the fit as an npy file.

  • verbose (bool (optional)) – Show the fit on a graph. Default is False.

Returns:

L – The path length through the sphere in mm. If the path lies outside the sphere, returns 0.

Return type:

float

calibrateAttenuation.getPathLengthThroughSphere(p, r, sdd, sod)[source]

Get the path length of a ray passing through a sphere. Schematic diagram is provided in the supplementary methods section of the first (only?) radiosphere paper.

Parameters:

  • p (float) – The distance in mm along the detector panel from the centre.

  • r (float) – The particle radius in mm.

  • sdd (float) – The source-detector distance in mm.

  • sod (float) – The source-object distance in mm.

Returns:

L – The path length through the sphere in mm. If the path lies outside the sphere, returns 0.

Return type:

float

DEM optimisation correction

DEM.DEM.DEM_step(posMMin, radiiMM, k=0.01)[source]

Lightweight DEM with assumption of same size spheres as mechanical regularisation

Parameters:

  • 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) – Particle radii for projection

  • k (float, optional) – Stiffness and timestep wrapped into one Default = 0.01

Returns:

posMM

Return type:

output positions

MercuryDPM loaders

NDDEM loaders

Indices and tables