Commit bb62b7c1 authored by Edward Andò's avatar Edward Andò
Browse files

new pixel search that checks all combinations and returns a position in im2...

new pixel search that checks all combinations and returns a position in im2 (simpler but less useful)
parent ed4db083
......@@ -762,14 +762,12 @@ def registerMultiscale(im1, im2, binStart, binStop=1, im1mask=None, PhiInit=None
return reg
def pixelSearch(imagette1, imagette2, imagette1mask=None, searchCentre=None, searchRange=[0, 0, 0, 0, 0, 0]):
def pixelSearch(imagette1, imagette2, imagette1mask=None, imagette2mask=None):
"""
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).
The normalised correlation coefficient (NCC) is computed for EVERY combination of the displacements of imagette1 defined within imagette2.
What is returned is the highest NCC value obtained and the position where it was obtained with respect to the origin in imagette2.
Values of NCC > 0.99 can generally be trusted.
......@@ -782,78 +780,47 @@ def pixelSearch(imagette1, imagette2, imagette1mask=None, searchCentre=None, sea
Imagette 2 is the bigger image inside which to search
imagette1mask : 3D numpy array of bools, optional
A mask for imagette1 to define which areas to correlate
A mask for imagette1 to define which pixels to include in the correlation
(True = Correlate these pixels, False = Skip),
Default = no mask
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 : 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
imagette2mask : 3D numpy array of bools, optional
A mask for imagette2 to define which pixels to include in the correlation
(True = Correlate these pixels, False = Skip),
Default = no mask
Returns
-------
Dictionary:
't' : 3-component vector
Z, Y, X displacements in pixels from searchCentre of the best position of the centre of imagette1
p : 3-component vector
Z, Y, X position with respect to origin in imagette2 of imagette1 to get best NCC
'cc' : float
Normalised correlation coefficient, ~0.5 is random, >0.99 is very good correlation.
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 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"
#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(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')
if imagette1mask is not None:
assert imagette1.shape == imagette1mask.shape, "spam.DIC.correlate.pixelSearch: imagette1mask ({}) should have same size as imagette1 ({})".format( imagette1mask.shape, imagette1.shape)
imagette1 = imagette1.astype('<f4')
imagette1[imagette1mask == 0] = numpy.nan
# 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]
print("searchRange in:", searchRange)
# Check that the search range doesn't go outside imagette2, if it does, redact it
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("searchRange out:", searchRange)
if imagette2mask is not None:
assert imagette2.shape == imagette2mask.shape, "spam.DIC.correlate.pixelSearch: imagette2mask ({}) should have same size as imagette2 ({})".format( imagette2mask.shape, imagette2.shape)
imagette2 = imagette2.astype('<f4')
imagette2[imagette2mask == 0] = numpy.nan
# Run the actual pixel search
# print imagette1.shape, imagette2.shape, searchCentre, searchRange
returns = numpy.zeros(4, dtype='<f4')
DICToolkit.pixelSearch(imagette1.astype("<f4"),
imagette2.astype("<f4"),
searchCentre.astype("<f4"),
numpy.array(searchRange).astype("<f4"),
returns)
# Collect and pack returns
return {'t': returns[0:3],
'cc': returns[3]}
return numpy.array(returns[0:3]), returns[3]
def globalCorrelation(im1, im2, mesh, convergenceCriterion=0.01, debugFiles=False, maxIterations=20, initialDisplacements=None, boundaryConditions=None, prefix="./"):
......
......@@ -75,6 +75,4 @@ void computeGMresidualAndPhase(py::array_t<float> im1,
void pixelSearch(py::array_t<float> im1,
py::array_t<float> im2,
py::array_t<float> startPos, /* Position in image 1 of the middle of the search range, Z, Y, X */
py::array_t<float> searchRange, /* minus and plus searching directions [ [ z-, z+ ], [ y-, y+ ], [ x-, x+ ] ] */
py::array_t<float> argoutdata);
......@@ -23,8 +23,6 @@ float InvSqrt(float x)
/* Image sizes, ZYX and images*/
void pixelSearch(py::array_t<float> im1Numpy,
py::array_t<float> im2Numpy,
py::array_t<float> startPosNumpy, /* Position in image 2 of the middle of the search range, Z, Y, X */
py::array_t<float> searchRangeNumpy, /* minus and plus searching directions [ [ z-, z+ ], [ y-, y+ ], [ x-, x+ ] ] */
py::array_t<float> argoutdataNumpy )
{
......@@ -32,10 +30,6 @@ void pixelSearch(py::array_t<float> im1Numpy,
float *im1 = (float*) im1Buf.ptr;
py::buffer_info im2Buf = im2Numpy.request();
float *im2 = (float*) im2Buf.ptr;
py::buffer_info startPosBuf = startPosNumpy.request();
float *startPos = (float*) startPosBuf.ptr;
py::buffer_info searchRangeBuf = searchRangeNumpy.request();
float *searchRange = (float*) searchRangeBuf.ptr;
py::buffer_info argoutdataBuf = argoutdataNumpy.request();
float *argoutdata = (float*) argoutdataBuf.ptr;
......@@ -51,7 +45,7 @@ void pixelSearch(py::array_t<float> im1Numpy,
// size_t index1, index2;
/* loop variables for 3D search range */
long zDisp, yDisp, xDisp;
size_t zTop, yTop, xTop;
/* loop variables for 3D CC calculation */
unsigned int z, y, x;
......@@ -85,100 +79,80 @@ void pixelSearch(py::array_t<float> im1Numpy,
/* Go through search range in im2 -- z, y, x positions are offsets of the window,
Consequently the first iteration here at z=y=x=0, is comparing im1 with the top
Corner of im2 */
// std::cout << "start pos:" << startPos[0] << " " << startPos[1] << " " << startPos[2] << std::endl;
for ( zDisp=(long)searchRange[0]; zDisp <= (long)searchRange[1]; zDisp++ )
// coordinates of the top corner of im1 in the coordinates of im2
for ( zTop=0; zTop <= im2z-im1z; zTop++ )
{
for ( yDisp = (long)searchRange[2]; yDisp <= (long)searchRange[3]; yDisp++ )
for ( yTop = 0; yTop <= im2y-im1y; yTop++ )
{
for ( xDisp = (long)searchRange[4]; xDisp <= (long)searchRange[5]; xDisp++ )
for ( xTop = 0; xTop <= im2x-im1x; xTop++ )
{
// std::cout << zDisp << " " << yDisp << " " << xDisp << std::endl;
/* Calculate this bit as far out as possible to benefit from recycling */
long zTop, yTop, xTop;
long zBot, yBot, xBot;
zTop = (long)startPos[0]+zDisp-(im1z-1)/2;
yTop = (long)startPos[1]+yDisp-(im1y-1)/2;
xTop = (long)startPos[2]+xDisp-(im1x-1)/2;
zBot = zTop+im1z;
yBot = yTop+im1y;
xBot = xTop+im1x;
// printf("zTop %i yTop %i xTop %i\n", zTop, yTop, xTop);
// // printf("zDisp %i yDisp %i xDisp %i\n", zDisp, yDisp, xDisp);
std::cout << "Tops: "<< zTop << " " << yTop << " " << xTop << "\n" << std::endl;
std::cout << "Bots: "<< zBot << " " << yBot << " " << xBot << "\n" << std::endl;
/* Check we're not outside the boundaries... */
if ( zTop >= 0 && yTop >= 0 && xTop >= 0 &&
zBot <= (int)im2z && yBot <= (int)im2y && xBot <= (int)im2x )
// std::cout << zTop << " " << yTop << " " << xTop << std::endl;
/* reset calculations */
/* three components to our NCC calculation (see documentation/C-remi.odt) */
float a,b,c;
a = b = c = 0;
/* CC calculation Loop z-first (for numpy) */
for ( z=0; z<im1z; z++ )
{
std::cout << "yo" << std::endl;
/* reset calculations */
/* three components to our NCC calculation (see documentation/C-remi.odt) */
float a,b,c;
a = b = c = 0;
/* CC calculation Loop z-first (for numpy) */
for ( z=0; z<im1z; z++ )
/* More speedups, precalculate slice offset*/
size_t zOffset1 = z *im1yXim1x;
size_t zOffset2 = (z+zTop)*im2yXim2x;
for ( y=0; y<im1y; y++ )
{
/* More speedups, precalculate slice offset*/
unsigned long zOffset1 = z *im1yXim1x;
unsigned long zOffset2 = (z+zTop)*im2yXim2x;
/* More speedups, precalculate column offset*/
size_t yOffset1 = y *im1x;
size_t yOffset2 = (y+yTop)*im2x;
for ( y=0; y<im1y; y++ )
for ( x=0; x<im1x; x++ )
{
/* More speedups, precalculate column offset*/
unsigned long yOffset1 = y *im1x;
unsigned long yOffset2 = (y+yTop)*im2x;
/* build index to 1D image1 */
size_t index1 = zOffset1 + yOffset1 + x;
for ( x=0; x<im1x; x++ )
/* build index to 1D image2 */
size_t index2 = zOffset2 + yOffset2 + (x+xTop);
/* 2015-10-22: EA -- skip NaNs in the reference image
* NaNs in C are not even equal to themselves, so we'll use this simple check. */
if ( im1[ index1 ] == im1[ index1 ] and im2[ index2 ] == im2[ index2 ] )
{
/* build index to 1D image1 */
size_t index1 = zOffset1 + yOffset1 + x;
/* 2015-10-22: EA -- skip NaNs in the reference image
* NaNs in C are not even equal to themselves, so we'll use this simple check. */
if ( im1[ index1 ] == im1[ index1 ] )
{
/* build index to 1D image2 */
size_t index2 = zOffset2 + yOffset2 + (x+xTop);
// fetch 1 pixel from both images
im1px = im1[ index1 ];
im2px = im2[ index2 ];
// Our little bits of the NCC
a = a + im1px * im2px;
b = b + im1px * im1px;
c = c + im2px * im2px;
}
// fetch 1 pixel from both images
im1px = im1[ index1 ];
im2px = im2[ index2 ];
// Our little bits of the NCC
a = a + im1px * im2px;
b = b + im1px * im1px;
c = c + im2px * im2px;
}
}
}
/* End CC calculation loop */
}
/* End CC calculation loop */
/* once the sums are done, add up and sqrt
* assemble bits and calculate CC */
/* once the sums are done, add up and sqrt
* assemble bits and calculate CC */
cc = a * InvSqrt( b * c );
cc = a * InvSqrt( b * c );
// cc = a / std::sqrt( b * c );
// printf( "\tC: pixel_search: cc = %f\n", cc );
// printf( "\t-> CC@(z=%i,y=%i,x=%i) = %f\n", z, y, x, cc );
/* If this cc is higher than the previous best, update our best... */
if ( cc > ccMax )
{
xMax = xDisp;
yMax = yDisp;
zMax = zDisp;
ccMax = cc;
/* If this cc is higher than the previous best, update our best... */
if ( cc > ccMax )
{
xMax = xTop;
yMax = yTop;
zMax = zTop;
ccMax = cc;
// printf("##### CC UPDATED #####");
// printf( "\t-> New CC_max@(z=%i,y=%i,x=%i) = %f\n", zMax, yMax, xMax, cc );
}
}
}
}
......
......@@ -59,6 +59,7 @@ class TestFunctionDVC(unittest.TestCase):
except OSError:
pass
def test_applyPhi(self):
a = numpy.random.rand(10, 10, 10)
Phi = spam.deformation.computePhi({'t': [0.0, 3.0, 3.0]})
......@@ -357,61 +358,53 @@ class TestFunctionDVC(unittest.TestCase):
for i in range(3):
self.assertAlmostEqual(t7[c][i] - t[c][i], 0, places=1)
def test_pixelSearch(self):
tVector = numpy.random.randint(1, 4, 3)
# print("test_pixelSearch(): initial random displacement:", tVector)
#print("test_pixelSearch(): initial random displacement:", tVector)
p = numpy.random.randint(15, 35, 3)
imDef = spam.DIC.applyPhi(im, Phi=spam.deformation.computePhi({"t": tVector}))
halfWindowSize = [5, 5, 5]
halfWindowSize = [5, 5, 5]
# search for the translation of a random point
subVolSlice1 = (slice(int(p[0] - halfWindowSize[0]), int(p[0] + halfWindowSize[0] + 1)),
slice(int(p[1] - halfWindowSize[1]), int(p[1] + halfWindowSize[1] + 1)),
slice(int(p[2] - halfWindowSize[2]), int(p[2] + halfWindowSize[2] + 1)))
subVolSliceRef = (slice(int(p[0] - halfWindowSize[0]), int(p[0] + halfWindowSize[0] + 1)),
slice(int(p[1] - halfWindowSize[1]), int(p[1] + halfWindowSize[1] + 1)),
slice(int(p[2] - halfWindowSize[2]), int(p[2] + halfWindowSize[2] + 1)))
# CASE 1: local coordinates
ps1 = spam.DIC.pixelSearch(im[subVolSlice1].copy(), imDef,
searchCentre=p,
searchRange=[-5, 5]*3)
self.assertTrue(ps1["cc"] > 0.9)
topDisp = [-5, -5, -5]
subVolSliceDef = (slice(int(p[0] - halfWindowSize[0] + topDisp[0]), int(p[0] + halfWindowSize[0] + 1 + 5)),
slice(int(p[1] - halfWindowSize[1] + topDisp[1]), int(p[1] + halfWindowSize[1] + 1 + 5)),
slice(int(p[2] - halfWindowSize[2] + topDisp[2]), int(p[2] + halfWindowSize[2] + 1 + 5)))
# CASE 0: one image with itself:
ps0 = spam.DIC.pixelSearch(im[subVolSliceRef].copy(), im[subVolSliceRef].copy())
self.assertTrue(ps0[1] > 0.9)
for i in range(3):
self.assertAlmostEqual(tVector[i], ps1["t"][i], places=4)
self.assertAlmostEqual(0, ps0[0][i], places=4)
# CASE 1: local coordinates
ps1 = spam.DIC.pixelSearch(im[subVolSlice1].copy(), imDef,
searchCentre=p,
searchRange=[-5, 5]*3)
self.assertTrue(ps1["cc"] > 0.9)
ps1 = spam.DIC.pixelSearch(im[subVolSliceRef].copy(), imDef[subVolSliceDef].copy())
#print(ps1)
self.assertTrue(ps1[1] > 0.9)
for i in range(3):
self.assertAlmostEqual(tVector[i], ps1["t"][i], places=4)
self.assertAlmostEqual(tVector[i], ps1[0][i]+topDisp[0], places=4)
# CASE 1: local coordinates with mask
imagette1mask = numpy.ones_like(im[subVolSlice1])
# CASE 2: local coordinates with mask
imagette1mask = numpy.ones_like(im[subVolSliceRef])
imagette1mask[(0,-1), :] = 0
imagette1mask[:, (0,-1)] = 0
ps1b = spam.DIC.pixelSearch(im[subVolSlice1].copy(), imDef,
searchCentre=p,
searchRange=[-5, 5]*3,
imagette1mask=imagette1mask)
self.assertTrue(ps1b["cc"] > 0.9)
ps2 = spam.DIC.pixelSearch(im[subVolSliceRef].copy(), imDef[subVolSliceDef].copy(), imagette1mask=imagette1mask)
self.assertTrue(ps2[1] > 0.9)
for i in range(3):
self.assertAlmostEqual(tVector[i], ps1b["t"][i], places=4)
# CASE 1b: wrong searchRange and no searchCentre
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=[-5, 5]*3,
# globalCoords = True)
self.assertAlmostEqual(tVector[i], ps2[0][i]+topDisp[0], places=4)
# self.assertTrue(ps2["cc"]>0.9)
# for i in range(3):
# self.assertAlmostEqual(tVector[i], ps2["t"][i], places=4)
# CASE XX: Catch asertion error when im1 is bigger than im2
self.assertRaises(AssertionError, spam.DIC.pixelSearch, imDef, im[subVolSliceRef].copy())
# Test for imagette2.shape >= imagette1.shape
self.assertRaises(AssertionError, spam.DIC.pixelSearch, imDef, im[subVolSlice1].copy(), imDef, searchRange=[-5, 5]*3)
def test_makeGrid(self):
# 3D image
......
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