Commit 1c238ccd authored by Edward Andò's avatar Edward Andò
Browse files

[skip-ci] documented and safened pxielSearch()

parent 4848890d
......@@ -38,7 +38,7 @@ import os
# Define argument parser object
parser = argparse.ArgumentParser(description="spam-register "+spam.helpers.optionsParser.GLPv3descriptionHeader +\
parser = argparse.ArgumentParser(description="spam-pixelSearch "+spam.helpers.optionsParser.GLPv3descriptionHeader +\
"This script performs a pixel search from im1 to im2\n",
formatter_class=argparse.RawTextHelpFormatter)
......@@ -57,9 +57,7 @@ minMaskCoverage = args.MASK_COVERAGE
greyThreshold = [args.GREY_LOW_THRESH, args.GREY_HIGH_THRESH]
# Fill in search range as a dictionary
searchRange = {'zRange': [args.SEARCH_RANGE[0], args.SEARCH_RANGE[1]],
'yRange': [args.SEARCH_RANGE[2], args.SEARCH_RANGE[3]],
'xRange': [args.SEARCH_RANGE[4], args.SEARCH_RANGE[5]]}
searchRange = [args.SEARCH_RANGE[0], args.SEARCH_RANGE[1], args.SEARCH_RANGE[2], args.SEARCH_RANGE[3], args.SEARCH_RANGE[4], args.SEARCH_RANGE[5]]
# Load reference image
im1 = tifffile.imread(args.im1.name)
......
......@@ -548,7 +548,6 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiRigid=False, PhiInitBinRat
'deltaPhiNorm': deltaPhiNorm}
def registerMultiscale(im1, im2, binStart, binStop=1, im1mask=None, PhiInit=None, PhiRigid=False, PhiInitBinRatio=1.0, margin=None, maxIterations=100, deltaPhiMin=0.0001, updateGradient=False, interpolationOrder=1, interpolator='C', verbose=False, imShowProgress=False, forceChangeScale=False):
"""
Perform multiscale subpixel image correlation between im1 and im2.
......@@ -763,9 +762,16 @@ def registerMultiscale(im1, im2, binStart, binStop=1, im1mask=None, PhiInit=None
return reg
def pixelSearch(imagette1, imagette2, searchCentre=None, searchRange={'zRange': [-2, 2], 'yRange': [-2, 2], 'xRange': [-2, 2]}):
def pixelSearch(imagette1, imagette2, searchCentre=None, searchRange=[0, 0, 0, 0, 0, 0]):
"""
Interface to Pixel Search C code.
This function performs a pixel-by-pixel search in 3D for a small reference imagette1 within a larger imagette2.
An optional searchCentre is defined with respect to imagette2 (otherwise imagette2's centre is taken).
For every combination of the displacements of imagette1 defined in the searchRange,
the Normalised Correlation Coefficient is computed for imagette1 in that position
in imagette2.
The result is the displacement, with respect to searchCentre, of the highest NCC measured (which is also returned).
Values of NCC > 0.99 can generally be trusted.
Parameters
----------
......@@ -773,56 +779,57 @@ def pixelSearch(imagette1, imagette2, searchCentre=None, searchRange={'zRange':
Image 1 is the smaller reference image
imagette2 : 3D numpy array of floats
Image 2 is the bigger image
Image 2 is the bigger image inside which to search
searchCentre : 3-component list or numpy array, optional
The point in im2 around which the search range is defined
Default is middle of imagette2
searchRange : dictionary, optional
This dictionary contains 3 keys: 'zRange', 'yRange', and 'xRange',
Each of which contains a list with two items.
Default +-2 pixels in every direction
searchRange : 6-component list or numpy array of ints, optional but highly reccommended
Contains -z, +z, -y, +y, -x, +x search ranges in pixels
Default: 0 in all directions
Returns
-------
Dictionary:
'transformation': dictionary
Dictionary containing:
't' : 3-component vector
't' : 3-component vector
Z, Y, X displacements in pixels from searchCentre of the best position of the centre of imagette1
Z, Y, X displacements in pixels
'cc':float
'cc' : float
Normalised correlation coefficient, ~0.5 is random, >0.99 is very good correlation.
Note
----
It important to remember that the C code runs MUCH faster in its current incarnation when it has
a cut-out im2 to deal with (this is related to processor optimistaions).
When the globalCoords flag is true we should search with all of im2 in memory, but (default) we will cut it out.
"""
if searchCentre is None:
searchCentre = (numpy.array(imagette2.shape, dtype='<f4') - 1) / 2.0
else:
searchCentre = numpy.array(searchCentre, dtype='<f4')
# If dictionary overwrite with numpy array
if type(searchRange) is dict:
searchRange = numpy.array([searchRange['zRange'],
searchRange['yRange'],
searchRange['xRange']], dtype='<f4')
else:
print("correlate.pixelSearch(): searchRange should be a dict with 'zRange', 'yRange', and 'xRange' entries. Exiting")
return
#assert(numpy.all(imagette2.shape >= imagette1.shape)), "spam.DIC.correlate.pixelSearch(): imagette2 should be bigger or equal to imagette1 in all dimensions"
# Not enough space to move im1 around
if numpy.any(imagette2.shape < imagette1.shape):
return {'transformation': {'t': numpy.array([numpy.nan, numpy.nan, numpy.nan])},
'cc': 0.0}
#Note
#----
#It important to remember that the C code runs A BIT faster in its current incarnation when it has
#a cut-out im2 to deal with (this is related to processor optimistaions).
#Cutting out imagette2 to just fit around the search range might save a bit of time
assert(len(searchRange) == 6), "spam.DIC.correlate.pixelSearch: searchRange should have 6 items"
assert(numpy.all(imagette2.shape >= imagette1.shape)), "spam.DIC.correlate.pixelSearch(): imagette2 should be bigger or equal to imagette1 in all dimensions"
## Not enough space to move im1 around
#if numpy.any(imagette2.shape < imagette1.shape):
#return {'t': numpy.array([numpy.nan, numpy.nan, numpy.nan]),
#'cc': 0.0}
if searchCentre is None: searchCentre = numpy.round((numpy.array(imagette2.shape, dtype='<f4') - 1) / 2).astype(int)
else: searchCentre = numpy.array(searchCentre, dtype='<f4')
# Compute HWS and offsets
imagette1halfWindowSize = numpy.round((numpy.array(imagette1.shape, dtype='<f4') - 1) / 2).astype(int)
topDiff = searchCentre.astype(int) - imagette1halfWindowSize + searchRange[0::2]
botDiff = numpy.array(imagette2.shape) - searchCentre.astype(int) - imagette1halfWindowSize - searchRange[1::2]
# Check that the search range doesn't go outside imagette2, if it does, redact it
#print("SR in", searchRange)
if topDiff[0] < 0: searchRange[0] -= topDiff[0]
if topDiff[1] < 0: searchRange[2] -= topDiff[1]
if topDiff[2] < 0: searchRange[4] -= topDiff[2]
if botDiff[0] < 0: searchRange[1] += botDiff[0]
if botDiff[1] < 0: searchRange[3] += botDiff[1]
if botDiff[2] < 0: searchRange[5] += botDiff[2]
#print("SR out", searchRange)
# Run the actual pixel search
# print imagette1.shape, imagette2.shape, searchCentre, searchRange
......@@ -830,11 +837,11 @@ def pixelSearch(imagette1, imagette2, searchCentre=None, searchRange={'zRange':
DICToolkit.pixelSearch(imagette1.astype("<f4"),
imagette2.astype("<f4"),
searchCentre.astype("<f4"),
searchRange.astype("<f4"),
numpy.array(searchRange).astype("<f4"),
returns)
# Collect and pack returns
return {'transformation': {'t': returns[0:3]},
return {'t': returns[0:3],
'cc': returns[3]}
......
......@@ -372,33 +372,29 @@ class TestFunctionDVC(unittest.TestCase):
# CASE 1: local coordinates
ps1 = spam.DIC.pixelSearch(im[subVolSlice1].copy(), imDef,
searchCentre=p,
searchRange={"zRange": [-5, 5], "yRange": [-5, 5], "xRange": [-5, 5]},
searchRange=[-5, 5]*3,
)
# print(ps1)
self.assertTrue(ps1["cc"] > 0.9)
for i in range(3):
self.assertAlmostEqual(tVector[i], ps1["transformation"]["t"][i], places=4)
self.assertAlmostEqual(tVector[i], ps1["t"][i], places=4)
# CASE 1b: wrong searchRange and no searchCentre
ps1b = spam.DIC.pixelSearch(im[subVolSlice1].copy(), imDef, searchRange=[-5, 5])
self.assertRaises(AssertionError, spam.DIC.pixelSearch, im[subVolSlice1].copy(), imDef, searchRange=[-5, 5])
# CASE 2: global coordinates
# ps2 = spam.DICpixelSearch(im[subVolSlice1].copy(), imDef,
searchCentre = p,
searchRange={"zRange": [-5, 5], "yRange": [-5, 5], "xRange": [-5, 5]},
searchRange=[-5, 5]*3,
# globalCoords = True)
# self.assertTrue(ps2["cc"]>0.9)
# for i in range(3):
# self.assertAlmostEqual(tVector[i], ps2["transformation"]["t"][i], places=4)
# self.assertAlmostEqual(tVector[i], ps2["t"][i], places=4)
# Test for imagette2.shape >= imagette1.shape
ps2 = spam.DIC.pixelSearch(imDef, im[subVolSlice1].copy(),
searchCentre=p,
searchRange={"zRange": [-5, 5], "yRange": [-5, 5], "xRange": [-5, 5]},
)
self.assertTrue(numpy.isnan(ps2["transformation"]["t"][0]))
self.assertRaises(AssertionError, spam.DIC.pixelSearch, imDef, im[subVolSlice1].copy(), imDef, searchRange=[-5, 5]*3)
def test_makeGrid(self):
# 3D image
......@@ -409,64 +405,6 @@ class TestFunctionDVC(unittest.TestCase):
nodePositions, nodesDim2D = spam.DIC.grid.makeGrid((1, 50, 50), 10)
self.assertEqual((1, 4, 4), nodesDim2D)
def test_pixelSearchOnGrid(self):
# don't pass optional parameters
nodePositions = spam.DIC.grid.makeGrid(im.shape, 10)[0]
Phi = spam.deformation.computePhi({'t': [0, 2, 1]})
imDef = spam.DIC.applyPhi(im, Phi=Phi)
ps = spam.DIC.pixelSearchOnGrid(im, imDef,
nodePositions, 5,
searchRange={"zRange": [0, 0], "yRange": [-3, 3], "xRange": [-3, 3]})
self.assertEqual(ps["PhiField"][:, 0, -1].sum(), 0)
self.assertEqual(ps["PhiField"][:, 1, -1].sum(), nodePositions.shape[0] * 2)
self.assertEqual(ps["PhiField"][:, 2, -1].sum(), nodePositions.shape[0])
# force to fail for greyThreshold condition
psb = spam.DIC.pixelSearchOnGrid(im, imDef,
nodePositions, 5,
searchRange={"zRange": [0, 0], "yRange": [-3, 3], "xRange": [-3, 3]},
greyThreshold=[numpy.inf, numpy.inf])
for node in range(nodePositions.shape[0]):
self.assertTrue(numpy.allclose(psb["PhiField"][node, 0:3, 0:3], numpy.eye(3)))
self.assertEqual(numpy.isfinite(psb["PhiField"][:, 0:3, -1]).sum(), 0)
# pass default optional parameters
ps2 = spam.DIC.pixelSearchOnGrid(im, imDef,
nodePositions, 5,
searchRange={"zRange": [0, 0], "yRange": [-3, 3], "xRange": [-3, 3]},
im1mask=numpy.ones_like(im),
minMaskCoverage=0.5,
greyThreshold=[-numpy.inf, numpy.inf])
self.assertEqual(ps2["PhiField"][:, 0, -1].sum(), 0)
self.assertEqual(ps2["PhiField"][:, 1, -1].sum(), nodePositions.shape[0] * 2)
self.assertEqual(ps2["PhiField"][:, 2, -1].sum(), nodePositions.shape[0])
# small rotation to test correct creation of imagette1def
nodePositions = spam.DIC.grid.makeGrid(im.shape, 10)[0]
Phi = spam.deformation.computePhi({'r': [2, 0, 0]})
imDef = spam.DIC.applyPhi(im, Phi=Phi)
# Compute PhiField
PhiField = numpy.zeros((nodePositions.shape[0], 4, 4))
for nodeNumber, nodePosition in enumerate(nodePositions):
# I
PhiField[nodeNumber] = numpy.eye(4)
# F
PhiField[nodeNumber, 0:3, 0:3] = Phi[0:3, 0:3]
# t
PhiField[nodeNumber, 0:3, -1 ] = spam.deformation.decomposePhi(Phi.copy(), PhiCentre=imCentre, PhiPoint=nodePosition)['t']
ps3 = spam.DIC.pixelSearchOnGrid(im, imDef,
nodePositions, 5,
searchRange={"zRange": [-3, 3], "yRange": [-3, 3], "xRange": [-3, 3]},
PhiField=PhiField)
for nodeNumber, nodePosition in enumerate(nodePositions):
self.assertAlmostEqual(numpy.sum(PhiField[nodeNumber] - ps3["PhiField"][nodeNumber]), 0, places=1)
def test_registerOnGrid(self):
im = spam.datasets.loadSnow()
nodePositions = spam.DIC.grid.makeGrid(im.shape, 20)[0]
......
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