Commit 8d256574 authored by Edward Andò's avatar Edward Andò
Browse files

spam-pixelSearch script with pixelSearchOnGrid removed from grid.py

parent c21c64c7
Pipeline #53592 failed with stages
in 13 minutes and 14 seconds
This diff is collapsed.
......@@ -27,11 +27,9 @@ field, if there is a homogeneous background movement to measure
import numpy
import argparse
#import spam.helpers
#import spam.mesh
#import spam.label
import spam.helpers
import spam.DIC
#import spam.deformation
import spam.deformation
import tifffile
import os
......
......@@ -146,6 +146,7 @@ TTK_modules_pb = ["mesh", "label", "filters", "measurements", "kalisphera", "DIC
###############################################################
scripts = ['scripts/spam-mmr',
'scripts/spam-reg',
'scripts/spam-pixelSearch',
'scripts/spam-gdic',
"scripts/spam-ddic",
"scripts/spam-ldic",
......
......@@ -86,269 +86,6 @@ def makeGrid(imageSize, nodeSpacing):
return nodePositions, nodesDim
def pixelSearchOnGrid(im1, im2, nodePositions, halfWindowSize, searchRange, PhiField=None, minMaskCoverage=0.5, im1mask=None, greyThreshold=[-numpy.inf, numpy.inf], mpi=False):
"""
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 ]
mpi : bool, optional (default = False)
Are we being called by an MPI run?
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
"""
def getImagettes(nodeNumber, nodePositions, PhiField, searchRange, halfWindowSize, im1, im2, minMaskVolume, greyThreshold):
returnStatus = 2
initialDisplacement = PhiField[nodeNumber, 0:3, 3].astype(int)
# Add initial displacement guess to search range
searchRangeForThisNode = {'zRange': [searchRange['zRange'][0] + initialDisplacement[0], searchRange['zRange'][1] + initialDisplacement[0]],
'yRange': [searchRange['yRange'][0] + initialDisplacement[1], searchRange['yRange'][1] + initialDisplacement[1]],
'xRange': [searchRange['xRange'][0] + initialDisplacement[2], searchRange['xRange'][1] + initialDisplacement[2]]}
# point in im2 that we are searching around
searchCentre = [halfWindowSize[0] - searchRangeForThisNode['zRange'][0],
halfWindowSize[1] - searchRangeForThisNode['yRange'][0],
halfWindowSize[2] - searchRangeForThisNode['xRange'][0]]
# 2020-09-25 OS and EA: Prepare startStop array for imagette 1 to be extracted with new slicePadded
PhiNoDisp = PhiField[nodeNumber]
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)):
applyPhiPad = 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
applyPhiPad = int(0.25*numpy.ceil(max(halfWindowSize)))
startStopIm1 = [int(nodePositions[nodeNumber, 0] - halfWindowSize[0] - applyPhiPad), int(nodePositions[nodeNumber, 0] + halfWindowSize[0] + applyPhiPad + 1),
int(nodePositions[nodeNumber, 1] - halfWindowSize[1] - applyPhiPad), int(nodePositions[nodeNumber, 1] + halfWindowSize[1] + applyPhiPad + 1),
int(nodePositions[nodeNumber, 2] - halfWindowSize[2] - applyPhiPad), int(nodePositions[nodeNumber, 2] + halfWindowSize[2] + applyPhiPad + 1)]
# In either case, extract imagette1, now guaranteed to be the right size
imagette1padded = spam.helpers.slicePadded(im1, startStopIm1)
# If F is not the identity, apply F and undo crop
if numpy.allclose(PhiNoDisp, numpy.eye(4)):
# 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
imagette1paddedDef = spam.DIC.applyPhi(imagette1padded, PhiNoDisp)
# undo padding
imagette1def = imagette1paddedDef[applyPhiPad:-applyPhiPad,
applyPhiPad:-applyPhiPad,
applyPhiPad:-applyPhiPad]
# Make sure imagette is not 0-dimensional in any dimension
if numpy.all(numpy.array(imagette1def.shape) > 0):
if numpy.nanmean(imagette1def) > greyThreshold[0] and numpy.nanmean(imagette1def) < greyThreshold[1] and len(imagette1def.ravel()) > minMaskVolume:
# Slice for image 2
## 2020-09-25 OS and EA: Prepare startStop array for imagette 2 to be extracted with new slicePadded
## Extract it...
startStopIm2 = [int(nodePositions[nodeNumber, 0] - halfWindowSize[0] + searchRangeForThisNode['zRange'][0]),
int(nodePositions[nodeNumber, 0] + halfWindowSize[0] + searchRangeForThisNode['zRange'][1] + 1),
int(nodePositions[nodeNumber, 1] - halfWindowSize[1] + searchRangeForThisNode['yRange'][0]),
int(nodePositions[nodeNumber, 1] + halfWindowSize[1] + searchRangeForThisNode['yRange'][1] + 1),
int(nodePositions[nodeNumber, 2] - halfWindowSize[2] + searchRangeForThisNode['xRange'][0]),
int(nodePositions[nodeNumber, 2] + halfWindowSize[2] + searchRangeForThisNode['xRange'][1] + 1)]
imagette2 = spam.helpers.slicePadded(im2, startStopIm2)
# Failed minMaskVolume or greylevel condition
else:
returnStatus = -5
imagette1def = None
imagette2 = None
# Failed 0-dimensional imagette test
else:
returnStatus = -5
imagette1def = None
imagette2 = None
return {'imagette1': imagette1def,
'imagette2': imagette2,
'returnStatus': returnStatus,
'initialDisplacement': initialDisplacement,
'searchRangeForThisNode': searchRangeForThisNode,
'searchCentre': searchCentre}
if mpi:
import mpi4py.MPI
mpiComm = mpi4py.MPI.COMM_WORLD
mpiSize = mpiComm.Get_size()
mpiRank = mpiComm.Get_rank()
mpiStatus = mpi4py.MPI.Status()
boss = mpiSize - 1
numberOfWorkers = mpiSize - 1
workersActive = numpy.zeros(numberOfWorkers)
else:
# not mpi
numberOfWorkers = 1
workersActive = numpy.array([0])
# start setting up
numberOfNodes = nodePositions.shape[0]
# Check input sanity
if type(halfWindowSize) == int or type(halfWindowSize) == float:
halfWindowSize = [halfWindowSize] * 3
# Check minMaskVolume
minMaskVolume = int(minMaskCoverage * (1+halfWindowSize[0]*2)*
(1+halfWindowSize[1]*2)*
(1+halfWindowSize[2]*2))
# Check F field
if PhiField is None:
PhiField = numpy.zeros((numberOfNodes, 4, 4))
for nodeNumber in range(numberOfNodes):
PhiField[nodeNumber] = numpy.eye(4)
# Create pixelSearchCC vector
pixelSearchCC = numpy.zeros((numberOfNodes))
# Add nodes to a queue -- mostly useful for MPI
q = queue.Queue()
for node in range(numberOfNodes):
q.put(node)
finishedNodes = 0
writeReturns = False
print("\n\tStarting Pixel search")
widgets = [progressbar.FormatLabel(''), ' ', progressbar.Bar(), ' ', progressbar.AdaptiveETA()]
pbar = progressbar.ProgressBar(widgets=widgets, maxval=numberOfNodes)
pbar.start()
while finishedNodes != numberOfNodes:
# If there are workers not working, satify their requests...
# Note: this condition is alyas true if we are not in MPI and there are jobs to do
if workersActive.sum() < numberOfWorkers and not q.empty():
worker = numpy.where(workersActive == False)[0][0]
# Get the next node off the queue
nodeNumber = q.get()
imagetteReturns = getImagettes(nodeNumber, nodePositions, PhiField, searchRange, halfWindowSize, im1, im2, minMaskVolume, greyThreshold)
if imagetteReturns['returnStatus'] == 2:
if mpi:
# build message for pixel search worker
m = {'nodeNumber': nodeNumber,
'im1': imagetteReturns['imagette1'],
'im2': imagetteReturns['imagette2'],
'searchRangeForThisNode': imagetteReturns['searchRangeForThisNode'],
'searchCentre': imagetteReturns['searchCentre'],
'initialDisplacement': imagetteReturns['initialDisplacement']
}
# print "\tBoss: sending node {} to worker {}".format( nodeNumber, worker )
mpiComm.send(m, dest=worker, tag=1)
# Mark this worker as working
workersActive[worker] = True
# NOTE: writeReturns is defined later when receiving messages
else: # Not MPI
returns = spam.DIC.correlate.pixelSearch(imagetteReturns['imagette1'],
imagetteReturns['imagette2'],
searchRange=imagetteReturns['searchRangeForThisNode'],
searchCentre=imagetteReturns['searchCentre'])
initialDisplacement = imagetteReturns['initialDisplacement']
writeReturns = True
else: # Regardless of MPI or single proc
# Failed to extract imagettes or something
pixelSearchCC[nodeNumber] = 0.0
finishedNodes += 1
PhiField[nodeNumber, 0:3, 0:3] = numpy.eye(3)
PhiField[nodeNumber, 0:3, 3] = numpy.nan
# Otherwise spend time looking waiting for replies from workers
elif mpi:
message = mpiComm.recv(source=mpi4py.MPI.ANY_SOURCE, tag=2, status=mpiStatus)
tag = mpiStatus.Get_tag()
if tag == 2:
worker = message[0]
nodeNumber = message[1]
returns = message[2]
initialDisplacement = message[3]
# print "\tBoss: received node {} from worker {}".format( nodeNumber, worker )
workersActive[worker] = False
writeReturns = True
else:
print("\tBoss: Don't recognise tag ", tag)
# If we have new DVC returns, save them in our output matrices
if writeReturns:
finishedNodes += 1
writeReturns = False
# set translation for this node
PhiField[nodeNumber, 0:3, 3] = numpy.array(returns['transformation']['t'])
pixelSearchCC[nodeNumber] = returns['cc']
widgets[0] = progressbar.FormatLabel(" CC={:0>7.5f} ".format(pixelSearchCC[nodeNumber]))
pbar.update(finishedNodes)
pbar.finish()
print("\n")
return {'PhiField': PhiField,
'pixelSearchCC': pixelSearchCC}
def registerOnGrid(im1, im2, nodePositions, halfWindowSize, PhiField=None, margin=None, maxIterations=None, deltaPhiMin=None, updateGradient=None, interpolationOrder=None, minMaskCoverage=0.5, im1mask=None, greyThreshold=[-numpy.inf, numpy.inf], mpi=False, killWorkersWhenDone=True):
"""
This function handles grid-based local correlation, performing a "register" subpixel refinement.
......
......@@ -2075,3 +2075,239 @@ def register(parser):
return args
def pixelSearch(parser):
parser.add_argument('im1',
metavar='im1',
type=argparse.FileType('r'),
help="Greyscale image of reference state for correlation")
parser.add_argument('im2',
metavar='im2',
type=argparse.FileType('r'),
help="Greyscale image of deformed state for correlation")
parser.add_argument('-lab1',
'--labelledFile1',
dest='LAB1',
default=None,
type=argparse.FileType('r'),
help="Path to tiff file containing a labelled image 1 that defines zones to correlate. Disactivates -hws and -ns options")
parser.add_argument('-mf1',
'--maskFile1',
dest='MASK1',
default=None,
type=argparse.FileType('r'),
help="Path to tiff file containing the mask of image 1 -- masks zones not to correlate, which should be == 0")
parser.add_argument('-pf',
'-phiFile',
dest='PHIFILE',
default=None,
type=argparse.FileType('r'),
help="Path to TSV file containing the deformation function field (required)")
parser.add_argument('-pfb',
'--phiFile-bin-ratio',
type=int,
default=1,
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('-pfni',
'--neighbours-for-phi-field-interpolation',
type=int,
default=6,
dest='NEIGHBOURS',
help="Number of neighbours for field interpolation. Default = 6")
#parser.add_argument('-regs',
#'--subtract-registration',
#action="store_true",
#dest='REGSUB',
#help='Subtract rigid part of input registration from output displacements? Only works if you load a registration TSV. Default = False')
parser.add_argument('-sr',
'--search-range',
nargs=6,
type=int,
default=[-3, 3, -3, 3, -3, 3],
dest='SEARCH_RANGE',
help='Z- Z+ Y- Y+ X- X+ ranges (in pixels) for the pxiel search. Requires pixel search to be activated. Default = +-3px')
#parser.add_argument('-psf',
#'--pixel-search-filter',
#type=int,
#default=0,
#dest='PS_FILTER',
#help='Median filter pixel search results. Default = 0')
# Default: node spacing equal in all three directions
parser.add_argument('-ns',
'--node-spacing',
nargs=1,
type=int,
default=None,
dest='NS',
help="Node spacing in pixels (assumed equal in all 3 directions -- see -ns3 for different setting). Default = 10px")
# Possible: node spacing different in all three directions
parser.add_argument('-ns3',
'--node-spacing-3',
nargs=3,
type=int,
default=None,
dest='NS',
help="Node spacing in pixels (different in 3 directions). Default = 10, 10, 10px")
# Default: window size equal in all three directions
parser.add_argument('-hws',
'--half-window-size',
nargs=1,
type=int,
default=None,
dest='HWS',
help="Half correlation window size, measured each side of the node pixel (assumed equal in all 3 directions -- see -hws3 for different setting). Default = 10 px")
# Possible: node spacing different in all three directions
parser.add_argument('-hws3',
'--half-window-size-3',
nargs=3,
type=int,
default=None,
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('-bb',
#'--binning-begin',
#type=int,
#default=4,
#dest='BIN_BEGIN',
#help='Initial binning to apply to input images for initial registration. Default = 4')
#parser.add_argument('-be',
#'--binning-end',
#type=int,
#default=1,
#dest='BIN_END',
#help='Binning level to stop at for initial registration. Default = 1')
parser.add_argument('-mc',
'--mask-coverage',
type=float,
default=0.5,
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('-glt',
'--grey-low-threshold',
type=float,
default=-numpy.inf,
dest='GREY_LOW_THRESH',
help="Grey threshold on mean of reference imagette BELOW which the correlation is not performed. Default = -infinity")
parser.add_argument('-ght',
'--grey-high-threshold',
type=float,
default=numpy.inf,
dest='GREY_HIGH_THRESH',
help="Grey threshold on mean of reference imagette ABOVE which the correlation is not performed. Default = infinity")
parser.add_argument('-od',
'--out-dir',
type=str,
default=None,
dest='OUT_DIR',
help="Output directory, default is the dirname of gmsh file")
parser.add_argument('-pre',
'--prefix',
type=str,
default=None,
dest='PREFIX',
help='Prefix for output files (without extension). Default is basename of mesh file')
parser.add_argument('-def',
'--save-deformed-image1',
action="store_true",
default=False,
dest='DEF',
help="Activate the saving of a deformed image 1 (as <im1>-reg-def.tif)")
parser.add_argument('-vtk',
'--VTKout',
action="store_true",
dest='VTK',
help='Activate VTK output format. Default = False')
parser.add_argument('-notsv',
'--noTSVout',
action="store_false",
dest='TSV',
help='Disactivate TSV output format. Default = False')
parser.add_argument('-tif',
'-tiff',
'--TIFFout',
'--TIFout',
action="store_true",
dest='TIFF',
help='Activate TIFFoutput format. Default = False')
args = parser.parse_args()
# If we have no out dir specified, deliver on our default promise -- this can't be done inline before since parser.parse_args() has not been run at that stage.
if args.OUT_DIR is None:
args.OUT_DIR = os.path.dirname(args.im1.name)
# However if we have no dir, notice this and make it the current directory.
if args.OUT_DIR == "":
args.OUT_DIR = "./"
else:
# Check existence of output directory
try:
if args.OUT_DIR:
os.makedirs(args.OUT_DIR)
else:
args.DIR_out = os.path.dirname(args.im1.name)
except OSError:
if not os.path.isdir(args.OUT_DIR):
raise
# Catch interdependent node spacing and correlation window sizes
if args.NS is None:
print("\nUsing default node spacing: "),
if args.HWS is None:
print("2x default half window size"),
args.HWS = [10]
print("({}) which is".format(args.HWS[0])),
args.NS = [args.HWS[0] * 2]
else:
print("2x user-set half window size"),
if len(args.HWS) == 1:
print("({}) which is".format(args.HWS[0])),
args.NS = [int(args.HWS[0] * 2)]
elif len(args.HWS) == 3:
print("({} -- selecting smallest) which is".format(args.HWS)),
args.NS = [int(min(args.HWS) * 2)]
print(args.NS)
# Catch 3D options
if len(args.NS) == 1:
args.NS = [args.NS[0], args.NS[0], args.NS[0]]
if len(args.HWS) == 1:
args.HWS = [args.HWS[0], args.HWS[0], args.HWS[0]]
if args.LAB1 is not None:
print("I have been passed a labelled image and so I am disactivating node spacing and half-window size")
args.HWS = None
args.NS = None
# Output file name prefix
if args.PREFIX is None:
args.PREFIX = os.path.splitext(os.path.basename(args.im1.name))[0] + "-" + os.path.splitext(os.path.basename(args.im2.name))[0] + "-pixelSearch"
return args
Markdown is supported
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