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

[skip-ci] spam-pixelSearch with mask

parent 0eb68345
......@@ -54,7 +54,6 @@ for key in sorted(argsDict):
# Do something about these variables later...
halfWindowSize = numpy.array(args.HWS)
minMaskCoverage = args.MASK_COVERAGE
greyThreshold = [args.GREY_LOW_THRESH, args.GREY_HIGH_THRESH]
# Fill in search range as a dictionary
......@@ -82,15 +81,18 @@ if args.LAB1 is not None:
lab1 = tifffile.imread(args.LAB1.name).astype(spam.label.labelType)
boundingBoxes = spam.label.boundingBoxes(lab1)
nodePositions = spam.label.centresOfMass(lab1, boundingBoxes=boundingBoxes)
numberOfNodes = nodePositions.shape[0]
### Otherwise we are in node spacing and half-window size mode
else:
nodePositions, nodesDim = spam.DIC.makeGrid(im1.shape, args.NS)
boundingBoxes = numpy.zeros((nodePositions.shape[0], 6), dtype=int)
boundingBoxes[0::2] = nodePositions-halfWindowSize
boundingBoxes[1::2] = nodePositions+halfWindowSize
print(boundingBoxes)
exit()
# 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()
if args.MASK1:
im1mask = tifffile.imread(args.MASK1.name) != 0
......@@ -237,23 +239,22 @@ Returns
For each node, the NCC score obtained
"""
def getImagettes(nodeNumber, nodePositions, PhiField, searchRange, halfWindowSize, im1, im2, minMaskVolume, greyThreshold):
def getImagettes(nodePosition, Phi, searchRange, boundingBox, im1, im2, im1mask, minMaskCoverage, greyThreshold):
returnStatus = 2
initialDisplacement = PhiField[nodeNumber, 0:3, 3].astype(int)
initialDisplacement = Phi[0:3, 3].astype(int)
# Add initial displacement guess to search range
searchRangeForThisNode = searchRange
# point in im2 that we are searching around
#searchCentre = nodePosition + initialDisplacement # <-- this is only true in global coordinates
searchCentre = nodePosition - boundingBox[0::2] - searchRange[0::2]
#searchCentre = halfWindowSize - searchRange[0::2]
#Add initial displacement to searchRange
searchRangeForThisNode[0::2] += initialDisplacement
searchRangeForThisNode[1::2] += initialDisplacement
# point in im2 that we are searching around
searchCentre = halfWindowSize - searchRangeForThisNode[0::2]
searchRange[0::2] += initialDisplacement
searchRange[1::2] += initialDisplacement
# 2020-09-25 OS and EA: Prepare startStop array for imagette 1 to be extracted with new slicePadded
PhiNoDisp = PhiField[nodeNumber]
PhiNoDisp = Phi
PhiNoDisp[0:3,-1] = 0.0
# If F is not the identity, create a pad to be able to apply F to imagette 1
......@@ -261,11 +262,12 @@ def getImagettes(nodeNumber, nodePositions, PhiField, searchRange, halfWindowSiz
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)))
applyPhiPad = int(0.5*numpy.ceil(max(boundingBox[1::2]-boundingBox[0::2])))
startStopIm1 = [int(boundingBox[0] - applyPhiPad), int(boundingBox[1] + applyPhiPad),
int(boundingBox[2] - applyPhiPad), int(boundingBox[3] + applyPhiPad),
int(boundingBox[4] - applyPhiPad), int(boundingBox[5] + applyPhiPad)]
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)
......@@ -278,24 +280,43 @@ def getImagettes(nodeNumber, nodePositions, PhiField, searchRange, halfWindowSiz
imagette1paddedDef = spam.DIC.applyPhi(imagette1padded, PhiNoDisp)
# undo padding
imagette1def = imagette1paddedDef[applyPhiPad:-applyPhiPad,
applyPhiPad:-applyPhiPad,
applyPhiPad:-applyPhiPad]
applyPhiPad:-applyPhiPad,
applyPhiPad:-applyPhiPad]
### Check mask
if im1mask is None:
# no mask1 --> always pas this test (e.g., labelled image)
maskVolumeCondition = True
imagette1mask = None
else:
imagette1mask = spam.helpers.slicePadded(im1mask, boundingBox) != 0
maskVolumeCondition = imagette1mask.mean() > minMaskCoverage
# Make sure imagette is not 0-dimensional in any dimension
# Check minMaskVolume
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[0]),
int(nodePositions[nodeNumber, 0] + halfWindowSize[0] + searchRangeForThisNode[1] + 1),
int(nodePositions[nodeNumber, 1] - halfWindowSize[1] + searchRangeForThisNode[2]),
int(nodePositions[nodeNumber, 1] + halfWindowSize[1] + searchRangeForThisNode[3] + 1),
int(nodePositions[nodeNumber, 2] - halfWindowSize[2] + searchRangeForThisNode[4]),
int(nodePositions[nodeNumber, 2] + halfWindowSize[2] + searchRangeForThisNode[5] + 1)]
imagette2 = spam.helpers.slicePadded(im2, startStopIm2)
# Failed minMaskVolume or greylevel condition
# ------------ Grey threshold low --------------- and -------------- Grey threshold high -----------
if numpy.nanmean(imagette1def) > greyThreshold[0] and numpy.nanmean(imagette1def) < greyThreshold[1]:
if maskVolumeCondition:
# 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(boundingBox[0] + searchRange[0] ),
int(boundingBox[1] + searchRange[1] + 1),
int(boundingBox[2] + searchRange[2] ),
int(boundingBox[3] + searchRange[3] + 1),
int(boundingBox[4] + searchRange[4] ),
int(boundingBox[5] + searchRange[5] + 1)]
imagette2 = spam.helpers.slicePadded(im2, startStopIm2)
# Failed minMaskVolume condition
else:
returnStatus = -5
imagette1def = None
imagette2 = None
# Failed greylevel condition
else:
returnStatus = -5
imagette1def = None
......@@ -308,10 +329,11 @@ def getImagettes(nodeNumber, nodePositions, PhiField, searchRange, halfWindowSiz
imagette2 = None
return {'imagette1': imagette1def,
'imagette1mask': imagette1mask,
'imagette2': imagette2,
'returnStatus': returnStatus,
'initialDisplacement': initialDisplacement,
'searchRangeForThisNode': searchRangeForThisNode,
'searchRange': searchRange,
'searchCentre': searchCentre}
mpi = False
......@@ -333,13 +355,6 @@ else:
numberOfWorkers = 1
workersActive = numpy.array([0])
# start setting up
numberOfNodes = nodePositions.shape[0]
# Check minMaskVolume
minMaskVolume = int(minMaskCoverage * (1+halfWindowSize[0]*2)*
(1+halfWindowSize[1]*2)*
(1+halfWindowSize[2]*2))
# Create pixelSearchCC vector
pixelSearchCC = numpy.zeros((numberOfNodes))
......@@ -364,17 +379,18 @@ while finishedNodes != numberOfNodes:
# Get the next node off the queue
nodeNumber = q.get()
imagetteReturns = getImagettes(nodeNumber, nodePositions, PhiField, searchRange, halfWindowSize, im1, im2, minMaskVolume, greyThreshold)
imagetteReturns = getImagettes(nodePositions[nodeNumber], PhiField[nodeNumber], searchRange, boundingBoxes[nodeNumber], im1, im2, im1mask, args.MASK_COVERAGE, 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']
m = {'nodeNumber': nodeNumber,
'im1': imagetteReturns['imagette1'],
'im2': imagetteReturns['imagette2'],
'im1mask': imagetteReturns['imagette1mask'],
'searchRange': imagetteReturns['searchRange'],
'searchCentre': imagetteReturns['searchCentre'],
'initialDisplacement': imagetteReturns['initialDisplacement']
}
# print "\tBoss: sending node {} to worker {}".format( nodeNumber, worker )
......@@ -388,8 +404,9 @@ while finishedNodes != numberOfNodes:
else: # Not MPI
returns = spam.DIC.correlate.pixelSearch(imagetteReturns['imagette1'],
imagetteReturns['imagette2'],
searchRange=imagetteReturns['searchRangeForThisNode'],
searchCentre=imagetteReturns['searchCentre'])
imagette1mask = imagetteReturns['imagette1mask'],
searchRange = imagetteReturns['searchRange'],
searchCentre = imagetteReturns['searchCentre'])
initialDisplacement = imagetteReturns['initialDisplacement']
writeReturns = True
......@@ -524,7 +541,8 @@ if args.VTK:
#if tag == 1:
#pixelSearchReturns = spam.DIC.correlate.pixelSearch(m['im1'],
#m['im2'],
#searchRange=m['searchRangeForThisNode'],
#imagette1mask=m['im1mask'],
#searchRange=m['searchRange'],
#searchCentre=m['searchCentre'])
## print "\t\tI am worker {} Sending result for node {}".format( mpiRank, node )
#mpiComm.send([mpiRank, m['nodeNumber'], pixelSearchReturns, m['initialDisplacement']], dest=boss, tag=2)
......
......@@ -762,7 +762,7 @@ def registerMultiscale(im1, im2, binStart, binStop=1, im1mask=None, PhiInit=None
return reg
def pixelSearch(imagette1, imagette2, searchCentre=None, searchRange=[0, 0, 0, 0, 0, 0]):
def pixelSearch(imagette1, imagette2, imagette1mask=None, searchCentre=None, searchRange=[0, 0, 0, 0, 0, 0]):
"""
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).
......@@ -776,10 +776,15 @@ def pixelSearch(imagette1, imagette2, searchCentre=None, searchRange=[0, 0, 0, 0
Parameters
----------
imagette1 : 3D numpy array of floats
Image 1 is the smaller reference image
Imagette 1 is the smaller reference image
imagette2 : 3D numpy array of floats
Image 2 is the bigger image inside which to search
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
(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
......@@ -787,7 +792,7 @@ def pixelSearch(imagette1, imagette2, searchCentre=None, searchRange=[0, 0, 0, 0
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
Default = 0 in all directions
Returns
-------
......@@ -814,7 +819,11 @@ def pixelSearch(imagette1, imagette2, searchCentre=None, searchRange=[0, 0, 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( imagette1.shape, imagette1mask.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)
......
......@@ -2102,6 +2102,13 @@ def pixelSearch(parser):
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('-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('-pf',
'-phiFile',
dest='PHIFILE',
......@@ -2194,13 +2201,6 @@ def pixelSearch(parser):
#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,
......@@ -2301,9 +2301,11 @@ def pixelSearch(parser):
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
print("I have been passed a labelled image and so I am disactivating node spacing and half-window size and mask and setting mask coverage to 0")
args.HWS = None
args.NS = None
args.MASK1 = None
args.MASK_COVERAGE = 0
# Output file name prefix
if args.PREFIX is None:
......
......@@ -372,14 +372,31 @@ class TestFunctionDVC(unittest.TestCase):
# CASE 1: local coordinates
ps1 = spam.DIC.pixelSearch(im[subVolSlice1].copy(), imDef,
searchCentre=p,
searchRange=[-5, 5]*3,
)
# print(ps1)
searchRange=[-5, 5]*3)
self.assertTrue(ps1["cc"] > 0.9)
for i in range(3):
self.assertAlmostEqual(tVector[i], ps1["t"][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)
for i in range(3):
self.assertAlmostEqual(tVector[i], ps1["t"][i], places=4)
# CASE 1: local coordinates with mask
imagette1mask = numpy.ones_like(im[subVolSlice1])
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)
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])
......
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