Commit 88d3c2f5 authored by Edward Andò's avatar Edward Andò
Browse files

[skip-ci] spam-ldic now working, with no pixelSearch and -pfd option;...

[skip-ci] spam-ldic now working, with no pixelSearch and -pfd option; spam-pixelSearch only missing -skp
parent 445b31c5
Pipeline #58427 skipped
This diff is collapsed.
......@@ -22,18 +22,20 @@ You should have received a copy of the GNU General Public License along with
this program. If not, see <http://www.gnu.org/licenses/>.
"""
import numpy
import os
import argparse
import multiprocessing
import progressbar
import tifffile
import numpy
import spam.helpers
#import spam.mesh
import spam.label
import spam.DIC
import spam.deformation
import multiprocessing
import progressbar
import tifffile
import os
# Define argument parser object
......@@ -49,12 +51,6 @@ argsDict = vars(args)
for key in sorted(argsDict):
print("\t{}: {}".format(key, argsDict[key]))
# Do something about these variables later...
halfWindowSize = numpy.array(args.HWS)
greyThreshold = [args.GREY_LOW_THRESH, args.GREY_HIGH_THRESH]
# Fill in search range
searchRange = numpy.array([args.SEARCH_RANGE[0], args.SEARCH_RANGE[1], args.SEARCH_RANGE[2], args.SEARCH_RANGE[3], args.SEARCH_RANGE[4], args.SEARCH_RANGE[5]])
......@@ -90,10 +86,8 @@ else:
# start setting up
numberOfNodes = nodePositions.shape[0]
boundingBoxes = numpy.zeros((numberOfNodes, 6), dtype=int)
boundingBoxes[:, 0::2] = nodePositions - halfWindowSize
boundingBoxes[:, 1::2] = nodePositions + halfWindowSize
#print(boundingBoxes)
#exit()
boundingBoxes[:, 0::2] = nodePositions - args.HWS
boundingBoxes[:, 1::2] = nodePositions + args.HWS
if args.MASK1:
im1mask = tifffile.imread(args.MASK1.name) != 0
......@@ -255,7 +249,7 @@ def pixelSearchOneNode(nodeNumber):
imagetteReturns['returnStatus'] = 0
else:
imagetteReturns = spam.DIC.getImagettes(nodePositions[nodeNumber], PhiField[nodeNumber].copy(), searchRange.copy(), boundingBoxes[nodeNumber], im1, im2, im1mask, im2mask, args.MASK_COVERAGE, greyThreshold, twoD)
imagetteReturns = spam.DIC.getImagettes(im1, boundingBoxes[nodeNumber], PhiField[nodeNumber].copy(), im2, searchRange.copy(), im1mask=im1mask, im2mask=im2mask, minMaskCoverage=args.MASK_COVERAGE, greyThreshold=[args.GREY_LOW_THRESH, args.GREY_HIGH_THRESH], applyF=args.APPLY_F)
# If getImagettes was successful (size check and mask coverage check)
if imagetteReturns['returnStatus'] == 1:
......
......@@ -85,80 +85,113 @@ def makeGrid(imageSize, nodeSpacing):
return nodePositions, nodesDim
"""
This function handles grid-based local correlation, offering an initial rough dispalcement-only guess.
At the moment matching of windows is done with a Normalised-Correlation-Coefficient approach.
Parameters
----------
im1 : 3D numpy array
A 3D image of greylevels defining a reference configuration for the pixel search
im2 : 3D numpy array
A deformed 3D image of greylevels
nodePositions : nPoints*3 numpy array
Array containing Z, Y, X positions of each point in the grid, as returned by ``makeGrid`` for example
defined at im1
halfWindowSize : 3-item list or int
Size of subvolumes to perform the image correlation on, as a data range taken either side of the voxel on which the node is placed.
The subvolume will be 2*halfWindowSize + 1 pixels on each side.
A general recommendation is to make this half the node spacing
searchRange : dictionary
Search range as a dictionary containing 3 keys: 'zRange', 'yRange', and 'xRange',
Each of which contains a list with two items
PhiField : nPoints*4*4 numpy array, optional
Optional field of ``F`` transformation operators defined for each node.
Currently, only the translational components of F will be taken into account.
Default = No displacement
minMaskCoverage : float, optional
Minimum number of pixels in a subvolume for it to be correlated (only considered in the case of im1mask).
Default = 0.5
im1mask : 3D boolean numpy array, optional
A mask for im1 which is true in the zones to correlate.
Default = None
greyThreshold : list of two floats, optional
Threshold for the mean greylevel in each im1 subvolume.
If the mean is below the first value or above the second value, the grid point is not correlated.
Default = [ -inf, inf ]
twoD : bool, optional
Are passed imagette1 and imagette2 3D images?
Returns
-------
Dictionary containing:
Keys
PhiField : nNodes*4*4 numpy array of floats
For each node, the measured transformatio operator (displacement only)
pixelSearchCC : nNodes numpy array of floats
For each node, the NCC score obtained
"""
# WARNING ------------------------- VVVVVVVVVVV--- gets easily overwritten, pass a .copy()!
def getImagettes(nodePosition, Phi, searchRange, boundingBox, im1, im2, im1mask, im2mask, minMaskCoverage, greyThreshold, twoD=False):
def getImagettes(im1, boundingBox, Phi, im2, searchRange, im1mask=None, im2mask=None, minMaskCoverage=0.0, greyThreshold=[-numpy.inf, numpy.inf], applyF='all'):
"""
This function is responsible for extracting correlation windows ("imagettes") from two larger images (im1 and im2).
Both spam.correlate.pixelSearch and spam.correlate.register[Multiscale] want a fixed, smaller imagette1
and a larger imagette 2 in which to search/interpolate.
Parameters
----------
im1 : 3D numpy array
This is the large input reference image
boundingBox : 6-component numpy array of ints
This defines the windows to extract from im1.
It is defined as [ Zbot Ztop Ybot Ytop Xbot Xtop ]
Note: for 2D Zbot and Ztop should both be zero
Phi : 4x4 numpy array of floats
Phi matrix representing the movement of imagette1,
if not equal to `I`, imagette1 is deformed by the non-translation parts of Phi (F)
and the displacement is added to the search range (see below)
im2 : 3D numpy array
This is the large input deformed image
searchRange : 6-component numpy array of ints
This defines where imagette2 should be extracted with respect to imagette1's position in im1.
The 6 components correspond to [ Zbot Ztop Ybot Ytop Xbot Xtop ].
If Z, Y and X values are the same, then imagette2 will be displaced and the same size as imagette1.
If 'bot' is lower than 'top', imagette2 will be larger in that dimension
im1mask : 3D numpy array, optional
This needs to be same size as im1, but can be `None` if no mask is wanted.
This defines a mask for zones to correlate in im1, 0 means zone not to correlate
Default = None
im2mask : 3D numpy array, optional
This needs to be same size as im2, but can be `None` if no mask is wanted.
This defines a mask for zones to correlate in im2, 0 means zone not to correlate
Default = None
minMaskCoverage : float, optional
Threshold for imagette1 non-mask coverage, i.e. how much of imagette1 can be full of mask
before it is rejected with returnStatus = -5?
Default = 0
greyThreshold : two-component list of floats, optional
Bottom and top threshold values for mean value of imagette1 to reject it with returnStatus = -5
Default = no threshold
applyF : string, optional
If a non-identity Phi is passed, should the F be applied to the returned imagette1?
Options are: 'all', 'rigid', 'no'
Default = 'all'
Note: as of January 2021, it seems to make more sense to have this as 'all' for pixelSearch, and 'no' for local DIC
Returns
-------
Dictionary :
'imagette1' : 3D numpy array,
'imagette1mask': 3D numpy array of same size as imagette1 or None,
'imagette2': 3D numpy array, bigger or equal size to imagette1
'imagette2mask': 3D numpy array of same size as imagette2 or None,
'returnStatus': int,
Describes success in extracting imagette1 and imagette2.
If == 1 success, otherwise negative means failure.
'pixelSearchOffset': 3-component list of ints
Coordinates of the top of the pixelSearch range in im1, i.e., the displacement that needs to be
added to the raw pixelSearch output to make it a im1 -> im2 displacement
"""
returnStatus = 1
imagette1mask = None
imagette2mask = None
initialDisplacement = numpy.round(Phi[0:3, 3]).astype(int)
assert len(im1.shape) == len(im2.shape) == 3, "3D images needed for im1 and im2, if you have 2D images please pad them with im[numpy.newaxis, ...]"
if im1mask is not None: assert len(im1mask.shape) == 3, "3D image needed for im1mask, if you have 2D images please pad them with im[numpy.newaxis, ...]"
if im2mask is not None: assert len(im2mask.shape) == 3, "3D image needed for im2mask, if you have 2D images please pad them with im[numpy.newaxis, ...]"
# Detect 2D images
if im1.shape[0] == 1:
twoD = True
# Impose no funny business in z if in twoD
boundingBox[0:2] = 0
searchRange[0:2] = 0
else:
twoD = False
# Catch bad bounding boxes:
if numpy.all((boundingBox[1::2] - boundingBox[0::2]) > 0) or ( numpy.all(boundingBox[3::2] - boundingBox[2::2] > 0) and twoD):
# 2020-09-25 OS and EA: Prepare startStop array for imagette 1 to be extracted with new slicePadded
PhiNoDisp = Phi.copy()
PhiNoDisp[0:3,-1] -= initialDisplacement
#PhiNoDisp[0:3,-1] = 0.0
if applyF == 'rigid':
PhiNoDisp = spam.deformation.computeRigidPhi(PhiNoDisp)
#PhiNoDisp[0:3,-1] = 0.0
# If F is not the identity, create a pad to be able to apply F to imagette 1
if numpy.allclose(PhiNoDisp, numpy.eye(4)):
if numpy.allclose(PhiNoDisp, numpy.eye(4)) or applyF == 'no':
applyPhiPad = (0, 0, 0)
else:
# 2020-10-06 OS and EA: Add a pad to each dimension of 25% of max(halfWindowSize) to allow space to apply F (no displacement) to imagette1
......@@ -175,13 +208,14 @@ def getImagettes(nodePosition, Phi, searchRange, boundingBox, im1, im2, im1mask,
imagette1padded = spam.helpers.slicePadded(im1, startStopIm1)
# If F is not the identity, apply F and undo crop
if numpy.allclose(PhiNoDisp, numpy.eye(4)):
if numpy.allclose(PhiNoDisp, numpy.eye(4)) or applyF == 'no':
# In this case there is is no padding (despite the name) and we can just keep going
imagette1def = imagette1padded
else:
# apply F to imagette 1 padded
if twoD: imagette1paddedDef = spam.DIC.applyPhiPython(imagette1padded, PhiNoDisp)
else: imagette1paddedDef = spam.DIC.applyPhi(imagette1padded, PhiNoDisp)
else: imagette1paddedDef = spam.DIC.applyPhi( imagette1padded, PhiNoDisp)
# undo padding
if twoD:
......@@ -200,7 +234,7 @@ def getImagettes(nodePosition, Phi, searchRange, boundingBox, im1, im2, im1mask,
imagette1mask = None
else:
imagette1mask = spam.helpers.slicePadded(im1mask, boundingBox + numpy.array([0, 1, 0, 1, 0, 1])) != 0
maskVolumeCondition = imagette1mask.mean() >= minMaskCoverage
maskVolumeCondition = (imagette1mask != 0).mean() >= minMaskCoverage
# Make sure imagette is not 0-dimensional in any dimension
......@@ -220,8 +254,9 @@ def getImagettes(nodePosition, Phi, searchRange, boundingBox, im1, im2, im1mask,
int(boundingBox[5] + initialDisplacement[2] + searchRange[5] + 1)]
imagette2 = spam.helpers.slicePadded(im2, startStopIm2)
if im2mask is not None: imagette2mask = spam.helpers.slicePadded(im2mask, startStopIm2)
if im2mask is not None:
imagette2mask = spam.helpers.slicePadded(im2mask, startStopIm2)
# Failed minMaskVolume condition
else:
returnStatus = -5
......
......@@ -42,6 +42,13 @@ def ldicParser(parser):
type=argparse.FileType('r'),
help="A space-separated list of two or more 3D greyscale tiff files to correlate, in order")
parser.add_argument('-np',
'--number-of-processes',
default=None,
type=int,
dest='PROCESSES',
help="Number of parallel processes to use. Default = multiprocessing.cpu_count()")
parser.add_argument('-mf1',
'--maskFile1',
dest='MASK1',
......@@ -70,6 +77,20 @@ def ldicParser(parser):
dest='PHIFILE_BIN_RATIO',
help="Ratio of binning level between loaded Phi file and current calculation. If the input Phi file has been obtained on a 500x500x500 image and now the calculation is on 1000x1000x1000, this should be 2. Default = 1")
parser.add_argument('-pfd',
'--phiFile-direct',
action="store_true",
#default=1,
dest='PHIFILE_DIRECT',
help="Trust the Phi file completely? This option ignores and overrides -pfni and requires same nodes in same positions. Default = False")
parser.add_argument('-pfni',
'--neighbours-for-phi-field-interpolation',
type=int,
default=6,
dest='NEIGHBOURS',
help="Number of neighbours for field interpolation. Default = 6")
parser.add_argument('-glt',
'--grey-low-threshold',
type=float,
......@@ -120,67 +141,63 @@ def ldicParser(parser):
dest='HWS',
help="Half correlation window size, measured each side of the node pixel (different in 3 directions). Default = 10, 10, 10px")
parser.add_argument('-nogp',
'--no-grid-point',
action="store_false",
dest='GRID_POINT',
help='Disactivates grid-point registration')
parser.add_argument('-gpm',
'--grid-point-margin',
parser.add_argument('-m',
'-mar',
'-margin',
nargs=1,
type=int,
default=[3],
dest='GRID_POINT_MARGIN',
help="Margin in pixels for grid-point registration. Default = 3 px")
dest='MARGIN',
help="Margin in pixels for the correlation of each local subvolume. Default = 3 px")
parser.add_argument('-gpm3',
'--grid-point-margin3',
parser.add_argument('-m3',
'-mar3',
'-margin3',
nargs=3,
type=int,
default=None,
dest='GRID_POINT_MARGIN',
help="Subpixel margin for grid-point registration. Default = [3, 3, 3] px")
dest='MARGIN',
help="Margin in pixels for the correlation of each local subvolume. Default = 3 px")
parser.add_argument('-gpi',
'--grid-point-max-iterations',
parser.add_argument('-it',
'--max-iterations',
type=int,
default=50,
dest='GRID_POINT_MAX_ITERATIONS',
help="Maximum iterations for grid-point registration. Default = 50")
dest='MAX_ITERATIONS',
help="Maximum iterations for local correlation. Default = 50")
parser.add_argument('-gpp',
'--grid-point-min-delta-phi',
parser.add_argument('-dp',
'--min-delta-phi',
type=numpy.float,
default=0.001,
dest='GRID_POINT_MIN_PHI_CHANGE',
help="Minimum change in Phi to consider grid-point registration as converged. Default = 0.001")
dest='MIN_DELTA_PHI',
help="Minimum change in Phi for local convergence. Default = 0.001")
parser.add_argument('-gpo',
'--grid-point-interpolation-order',
parser.add_argument('-o',
'--interpolation-order',
type=int,
default=1,
dest='GRID_POINT_INTERPOLATION_ORDER',
help="Interpolation order for grid-point registration. Default = 1")
dest='INTERPOLATION_ORDER',
help="Image interpolation order for local correlation. Default = 1, i.e., linear interpolation")
parser.add_argument('-gpmc',
'--grid-point-mask-coverage',
parser.add_argument('-mc',
'--mask-coverage',
type=float,
default=0.5,
dest='GRID_POINT_MASK_COVERAGE',
dest='MASK_COVERAGE',
help="In case a mask is defined, tolerance for a subvolume's pixels to be masked before it is skipped with RS=-5. Default = 0.5")
parser.add_argument('-gpug',
'--grid-point-update-gradient',
parser.add_argument('-ug',
'--update-gradient',
action="store_true",
dest='GRID_POINT_UPDATE_GRADIENT',
help='Update gradient in grid-point registration? More computation time but more robust and possibly fewer iterations.')
dest='UPDATE_GRADIENT',
help='Update gradient in local correlation? More computation time but sometimes more robust and possibly fewer iterations.')
parser.add_argument('-sef',
'--series-Ffile',
action="store_true",
dest='SERIES_PHIFILE',
help='During a total analysis, activate use of previous Ffield for next correlation')
#parser.add_argument('-sef',
#'--series-Ffile',
#action="store_true",
#dest='SERIES_PHIFILE',
#help='During a total analysis, activate use of previous Ffield for next correlation')
parser.add_argument('-sei',
'--series-incremental',
......@@ -210,9 +227,9 @@ def ldicParser(parser):
parser.add_argument('-tsv',
'--TSVout',
action="store_true",
action="store_false",
dest='TSV',
help='Activate TSV output format. Default = False')
help='Activate TSV output format. Default = True')
parser.add_argument('-tif',
'-tiff',
......@@ -286,19 +303,18 @@ def ldicParser(parser):
if len(args.HWS) == 1:
args.HWS = [args.HWS[0], args.HWS[0], args.HWS[0]]
if len(args.GRID_POINT_MARGIN) == 1:
args.GRID_POINT_MARGIN = [args.GRID_POINT_MARGIN[0], args.GRID_POINT_MARGIN[0], args.GRID_POINT_MARGIN[0]]
if len(args.MARGIN) == 1:
args.MARGIN = [args.MARGIN[0], args.MARGIN[0], args.MARGIN[0]]
if type(args.GRID_POINT_MAX_ITERATIONS) == list:
args.GRID_POINT_MAX_ITERATIONS = args.GRID_POINT_MAX_ITERATIONS[0]
if type(args.MAX_ITERATIONS) == list:
args.MAX_ITERATIONS = args.MAX_ITERATIONS[0]
# Catch and overwrite 2D options
if twoD:
args.NS[0] = 1
args.HWS[0] = 0
args.GRID_POINT_MARGIN[0] = 0
args.PSR[0] = 0
args.PSR[1] = 0
args.MARGIN[0] = 0
# Behaviour undefined for series run and im1 mask since im1 will change, complain and continue
if args.MASK1 is not None and args.SERIES_INCREMENTAL:
......@@ -318,12 +334,8 @@ def ldicParser(parser):
print("#############################################################")
print("#############################################################")
if args.SERIES_PHIFILE:
args.TSV = True
# just keep the name for this one
if args.PHIFILE is not None:
args.PHIFILE = args.PHIFILE.name
#if args.SERIES_PHIFILE:
#args.TSV = True
return args
......@@ -1815,7 +1827,7 @@ def register(parser):
'--update-gradient',
action="store_true",
dest='UPDATE_GRADIENT',
help='Update gradient during newton iterations? More computation time but more robust and possibly fewer iterations. Default = False')
help='Update gradient during newton iterations? More computation time but sometimes more robust and possibly fewer iterations. Default = False')
parser.add_argument('-it',
'--max-iterations',
......@@ -1836,7 +1848,7 @@ def register(parser):
type=int,
default=1,
dest='INTERPOLATION_ORDER',
help="Interpolation order for image interpolation. Default = 1")
help="Image interpolation order. Default = 1, i.e., linear interpolation")
parser.add_argument('-g',
'--graph',
......@@ -1976,6 +1988,14 @@ def pixelSearch(parser):
dest='NEIGHBOURS',
help="Number of neighbours for field interpolation. Default = 6")
parser.add_argument('-F',
'--apply-F',
type=str,
default='all',
dest='APPLY_F',
help="Apply the F part of Phi guess? Accepted values are:\n\t\"all\": apply all of F" +
"\n\t\"rigid\": apply rigid part (mostly rotation) \n\t\"no\": don't apply it \"all\" is default")
#parser.add_argument('-regs',
#'--subtract-registration',
#action="store_true",
......
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