grid example now complete in tests (spam-pixelSearch + spam-ldic +...

grid example now complete in tests (spam-pixelSearch + spam-ldic + spam-regularStrain). Registration checked in 2D

scripts/spam-ldic

@@ -18,8 +18,9 @@ You should have received a copy of the GNU General Public License along with
 this program.  If not, see <http://www.gnu.org/licenses/>.
 """
 
-from __future__ import print_function
 import os
+os.environ['OPENBLAS_NUM_THREADS'] = '1'
+
 import argparse
 import tifffile
 import multiprocessing
@@ -116,7 +117,7 @@ def correlateOneNode(nodeNumber):
         - deltaPhiNorm
     """
     PhiInit = PhiField[nodeNumber]
-    imagetteReturns = spam.DIC.getImagettes(im1, boundingBoxes[nodeNumber], PhiInit.copy(), im2, margin, im1mask=im1mask, minMaskCoverage=args.MASK_COVERAGE, greyThreshold=[args.GREY_LOW_THRESH, args.GREY_HIGH_THRESH], applyF='no')
+    imagetteReturns = spam.DIC.getImagettes(im1, nodePositions[nodeNumber], args.HWS, PhiInit.copy(), im2, margin, im1mask=im1mask, minMaskCoverage=args.MASK_COVERAGE, greyThreshold=[args.GREY_LOW_THRESH, args.GREY_HIGH_THRESH], applyF='no', twoD=twoD)
 
     if imagetteReturns['returnStatus'] == 1:
         # compute displacement that will be taken by the getImagettes

scripts/spam-pixelSearch

@@ -23,6 +23,8 @@ this program.  If not, see <http://www.gnu.org/licenses/>.
 """
 
 import os
+os.environ['OPENBLAS_NUM_THREADS'] = '1'
+
 import argparse
 import multiprocessing
 import progressbar
@@ -85,9 +87,6 @@ else:
     nodePositions, nodesDim = spam.DIC.makeGrid(im1.shape, args.NS)
     # start setting up
     numberOfNodes = nodePositions.shape[0]
-    boundingBoxes = numpy.zeros((numberOfNodes, 6), dtype=int)
-    boundingBoxes[:, 0::2] = nodePositions - args.HWS
-    boundingBoxes[:, 1::2] = nodePositions + args.HWS
 
     if args.MASK1:
         im1mask = tifffile.imread(args.MASK1.name) != 0
@@ -249,7 +248,7 @@ def pixelSearchOneNode(nodeNumber):
             imagetteReturns['returnStatus']        = 0
 
     else:
-        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)
+        imagetteReturns = spam.DIC.getImagettes(im1, nodePositions[nodeNumber], args.HWS, 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, twoD=twoD)
 
     # If getImagettes was successful (size check and mask coverage check)
     if imagetteReturns['returnStatus'] == 1:

tests/test_scripts.py

@@ -89,7 +89,7 @@ class testAll(unittest.TestCase):
 
         pass
 
-    def _test_gridDICtools(self):
+    def test_gridDICtools(self):
         # Make test directory
         if not os.path.isdir(testFolder):
             os.mkdir(testFolder)
@@ -115,7 +115,8 @@ class testAll(unittest.TestCase):
         tifffile.imsave(testFolder + "snow-def.tif", snowDef)
         os.chmod(testFolder + "snow-def.tif", 0o666)
         # Plus 2D version
-        tifffile.imsave(testFolder + "snow-def-2D.tif", spam.DIC.applyPhiPython(snowRef2D[numpy.newaxis, ...], Phi=Phi))
+        #tifffile.imsave(testFolder + "snow-def-2D.tif", spam.DIC.applyPhiPython(snowRef2D[numpy.newaxis, ...], Phi=Phi))
+        tifffile.imsave(testFolder + "snow-def-2D.tif", snowDef[50])
         os.chmod(testFolder + "snow-def-2D.tif", 0o666)
 
         snowDefOnlyDisp = spam.DIC.applyPhi(snowRef, Phi=spam.deformation.computePhi({'t': refTranslation}))
@@ -155,6 +156,14 @@ class testAll(unittest.TestCase):
                                            'returnStatus': 2,
                                            'deltaPhiNorm': 1})
 
+        # generate a fake "eye reg" initial guess which is close in 2D
+        spam.helpers.writeRegistrationTSV(testFolder + "snow-grid-2D-ereg.tsv", [0, snowRef.shape[1]-1/2, snowRef.shape[2]-1/2],
+                                          {'Phi': spam.deformation.computePhi({'r': 0.8*refRotation, 't': 0.8*refTranslation}),
+                                           'error': 0,
+                                           'iterations': 0,
+                                           'returnStatus': 2,
+                                           'deltaPhiNorm': 1})
+
         # generate a fake "eye reg" initial guess which is close in 2D
         spam.helpers.writeRegistrationTSV(testFolder + "snow-grid-2D-ereg-onlyRot.tsv", [0, snowRef.shape[1]-1/2, snowRef.shape[2]-1/2],
                                           {'Phi': spam.deformation.computePhi({'r': refRotation}),
@@ -166,7 +175,7 @@ class testAll(unittest.TestCase):
         #######################################################
         ### Step 2 check spam-reg functionality
         #######################################################
-        ## Just run a simple registration
+        ## Just run a simple registration, with a guess from "ereg"
         exitCode = subprocess.call(["spam-reg", "-bb", "2",
                                     testFolder + "snow-ref.tif", testFolder + "snow-def.tif",
                                     "-od", testFolder + "", "-pre", "snow-grid",
@@ -179,6 +188,24 @@ class testAll(unittest.TestCase):
         self.assertTrue(numpy.allclose(refTranslation, transformation['t'], atol=0.01))
         self.assertTrue(numpy.allclose(refRotation,    transformation['r'], atol=0.01))
 
+
+
+        ## Just run a simple registration in 2D!
+        exitCode = subprocess.call(["spam-reg", "-bb", "2",
+                                    testFolder + "snow-ref-2D.tif", testFolder + "snow-def-2D.tif",
+                                    "-od", testFolder + "", "-pre", "snow-grid-2D",
+                                    "-mar", "20",
+                                    "-pf", testFolder + "snow-grid-2D-ereg.tsv",
+                                    "-mf1", testFolder + "snow-mask-ref-2D.tif"
+                                    ])
+        self.assertEqual(exitCode, 0)
+        reg2dResult = spam.helpers.readCorrelationTSV(testFolder + "snow-grid-2D-registration.tsv")
+        transformation2d = spam.deformation.decomposePhi(reg2dResult['PhiField'][0])
+        #print(transformation['t'])
+        #print(transformation['r'])
+        self.assertTrue(numpy.allclose(refTranslation, transformation2d['t'], atol=0.01))
+        self.assertTrue(numpy.allclose(refRotation,    transformation2d['r'], atol=0.01))
+
         #######################################################
         ### Step 3 check spam-pixelSearch functionality
         #######################################################
@@ -280,89 +307,138 @@ class testAll(unittest.TestCase):
         self.assertTrue(numpy.isclose(refTranslation[2], numpy.median(PSresult3e['PhiField'][PSresult3e['returnStatus']==1,2,-1]), atol=0.1))
 
 
+        #######################################################
+        ### Step 4 check spam-pixelSearchPropagate
+        #######################################################
+        # ???
 
+        #######################################################
+        ### Step 5 check spam-ldic loading previous results
+        #######################################################
 
-        exit()
+        ### Step 5a: load registration
+        exitCode = subprocess.call(["spam-ldic",
+                                    "-pf", testFolder + "snow-grid-registration.tsv",
+                                    "-glt", "5000", "-hws", "10", "-ns", "10",
+                                    testFolder + "snow-ref.tif", testFolder + "snow-def.tif",
+                                    "-od", testFolder + "",
+                                    "-pre", "snow-grid-5a",
+                                    "-mf1", testFolder + "snow-mask-ref.tif",
+                                    "-mc", "1.0", "-tif"
+                                    ])
+        self.assertEqual(exitCode, 0)
+        LDICresult5a = spam.helpers.readCorrelationTSV(testFolder + "snow-grid-5a-ldic.tsv")
+        dims = LDICresult5a['fieldDims']
 
-        ### TODO: ldic here starting for our favourite pixel search result
+        # Did more than 80% of points converge?
+        returnStatusCropped = LDICresult5a['returnStatus'].reshape(dims)[2:-2, 2:-2, 2:-2]
+        self.assertTrue(0.8 < numpy.sum(returnStatusCropped==2)/returnStatusCropped.size)
 
-        # Check TSV for local results
-        #######################################################
-        TSV = spam.helpers.readCorrelationTSV(testFolder + "snow-ref-snow-def.tsv")
-        converged = TSV['returnStatus'] == 2
-        # The "majority" should converge
-        self.assertEqual(numpy.mean(converged.astype(float)) > 0.4, True)
+        # median Rotation within 0.5 deg?
+        PhiFieldCropped = LDICresult5a['PhiField'].reshape(dims[0], dims[1], dims[2], 4, 4)[2:-2, 2:-2, 2:-2, :, :]
+        self.assertTrue(numpy.isclose(refRotation[0], numpy.median(spam.deformation.decomposePhiField(PhiFieldCropped[returnStatusCropped==2], components='r')['r'][:,0]), atol=0.5))
 
-        Phis = TSV['PhiField']
-        decomposedPhisConverged = spam.deformation.decomposeFfield(Phis[:,0:3,0:3], components=['r'])
+        # And the z-displacement is low
+        self.assertTrue(numpy.isclose(0, numpy.median(PhiFieldCropped[returnStatusCropped==2,0,-1]), atol=1.0))
 
-        # There's Z-rotation so translations need to be decomposed properly... just check Z
-        self.assertAlmostEqual(refTranslation[0], numpy.mean(Phis[converged,0,-1]), places=1)
 
-        # Check only rotations, with a rotation there will be a non-trivial displacement field
-        for i in range(3):
-            self.assertAlmostEqual(refRotation[i],    numpy.mean(decomposedPhisConverged['r'][converged, i]), places=1)
+        ### Step 5b: load pixelSearch results with -pfd (same grid) + ug
+        exitCode = subprocess.call(["spam-ldic",
+                                    "-pf", testFolder + "snow-grid-3a-pixelSearch.tsv", "-pfd",
+                                    "-glt", "5000", "-hws", "10",
+                                    testFolder + "snow-ref.tif", testFolder + "snow-def.tif",
+                                    "-od", testFolder + "",
+                                    "-pre", "snow-grid-5b",
+                                    "-mf1", testFolder + "snow-mask-ref.tif",
+                                    "-mc", "1.0", "-ug"
+                                    ])
+        self.assertEqual(exitCode, 0)
+        LDICresult5b = spam.helpers.readCorrelationTSV(testFolder + "snow-grid-5b-ldic.tsv")
 
-        #######################################################
-        ## OK onto the regularStrain
-        #######################################################
-        # make sure it runs the help without error
-        exitCode = subprocess.call(["spam-regularStrain", "--help"], stdout=FNULL, stderr=FNULL)
+        # Did more than 80% of points converge?
+        self.assertTrue(0.8 < numpy.sum(LDICresult5b['returnStatus']==2)/len(LDICresult5b['returnStatus']))
+
+        # median Rotation within 0.5 deg?
+        self.assertTrue(numpy.isclose(refRotation[0], numpy.median(spam.deformation.decomposePhiField(LDICresult5a['PhiField'][LDICresult5a['returnStatus']>0], components='r')['r'][:,0]), atol=0.5))
+
+        # And the z-displacement is low
+        self.assertTrue(numpy.isclose(0, numpy.median(LDICresult5b['PhiField'][LDICresult5b['returnStatus']==2,0,-1]), atol=1.0))
+
+
+        ### Step "5c" check 2D mode
+        exitCode = subprocess.call(["spam-ldic",
+                                    "-pf", testFolder + "snow-grid-2D-registration.tsv",
+                                    "-glt", "5000", "-hws", "10", "-ns", "5",
+                                    testFolder + "snow-ref-2D.tif", testFolder + "snow-def-2D.tif",
+                                    "-od", testFolder + "", "-tif",
+                                    "-pre", "snow-grid-2d",
+                                    "-mf1", testFolder + "snow-mask-ref-2D.tif", "-mc", "1.0"
+                                    ])
         self.assertEqual(exitCode, 0)
+        LDICresult5c = spam.helpers.readCorrelationTSV(testFolder + "snow-grid-2d-ldic.tsv")
+
+        # Did more than 50% of points converge?
+        self.assertTrue(0.5 < numpy.sum(LDICresult5c['returnStatus']==2)/len(LDICresult5c['returnStatus']))
+
+        # median Rotation within 0.5 deg?
+        self.assertTrue(numpy.isclose(refRotation[0], numpy.median(spam.deformation.decomposePhiField(LDICresult5c['PhiField'][LDICresult5c['returnStatus']>0], components='r')['r'][:,0]), atol=0.5))
 
-        # Now let's calculate rotation vector with Geers elements (radius=1) and TSV output
-        exitCode = subprocess.call(["spam-regularStrain", "snow-ref-snow-def.tsv", "-comp", "r", "-tsv", "-mask", "-r", "1"],
+        # And the z-displacement is low
+        self.assertTrue(numpy.isclose(0, numpy.median(LDICresult5c['PhiField'][LDICresult5c['returnStatus']==2,0,-1]), atol=1.0))
+
+
+
+
+        #######################################################
+        ## Step 6 onto the regularStrain
+        #######################################################
+        ### Step 6a Now let's calculate rotation vector with Geers elements (radius=1) and TSV output
+        exitCode = subprocess.call(["spam-regularStrain",
+                                    "snow-grid-5a-ldic.tsv",
+                                    "-comp", "r",
+                                    "-tsv", "-mask",
+                                    "-r", "1",
+                                    "-od", testFolder, "-pre", "snow-grid-6a",
+                                    "-tif"],
                                    stdout=FNULL, stderr=FNULL)
         self.assertEqual(exitCode, 0)
 
         # Now read output
-        strain = spam.helpers.readStrainTSV("snow-ref-snow-def-strain-Geers.tsv")
-        dims = [strain['fieldDims'][0], strain['fieldDims'][1], strain['fieldDims'][2], 3]
+        strain6a = spam.helpers.readStrainTSV("snow-grid-6a-strain-Geers.tsv")
+        dims6a = [strain6a['fieldDims'][0], strain6a['fieldDims'][1], strain6a['fieldDims'][2], 3]
 
         # Check each rotation component while having a 2-pixel crop of the boundaries
         for i in range(3):
-            self.assertAlmostEqual(numpy.mean(strain['r'].reshape(dims)[2:-2, 2:-2, 2:-2, i]),
-                                   refRotation[i],
-                                   places=1)
-
-        # Now let's calculate small and large strains with Q8 elements and TSV output
-        exitCode = subprocess.call(["spam-regularStrain", "--Q8", "snow-ref-snow-def.tsv",
+            #print(numpy.nanmean(strain6a['r'].reshape(dims6a)[2:-2, 2:-2, 2:-2, i]))
+            self.assertTrue(numpy.isclose(numpy.nanmean(strain6a['r'].reshape(dims6a)[2:-2, 2:-2, 2:-2, i]),
+                                          refRotation[i],
+                                          atol=0.5))
+
+        ### Step 6b Now let's calculate small and large strains with Q8 elements and TSV output
+        exitCode = subprocess.call(["spam-regularStrain", "--Q8",
+                                    "snow-grid-5a-ldic.tsv",
                                     "-comp", "U", "e", "vol", "volss", "dev", "devss",
-                                    "-tsv", "-vtk", "-mask"],
+                                    "-tsv", "-vtk", "-mask",
+                                    "-od", testFolder, "-pre", "snow-grid-6b"],
                                    stdout=FNULL, stderr=FNULL)
         self.assertEqual(exitCode, 0)
-        strain = spam.helpers.readStrainTSV("snow-ref-snow-def-strain-Q8.tsv")
+        strain6b = spam.helpers.readStrainTSV("snow-grid-6b-strain-Q8.tsv")
+        dims6b = [strain6b['fieldDims'][0], strain6b['fieldDims'][1], strain6b['fieldDims'][2], 3, 3]
         # Check small strains, should be ~0
-        self.assertAlmostEqual(numpy.nanmean(strain['e']), 0, places=1)
+        self.assertTrue(numpy.allclose(numpy.nanmean(strain6b['e'].reshape(dims6b)[2:-2, 2:-2, 2:-2], axis=(0,1,2)), numpy.zeros((3,3)), atol=0.02))
 
-        Umean = numpy.nanmean(strain['U'], axis=0)
-        self.assertAlmostEqual(numpy.nanmean(Umean-numpy.eye(3)), 0, places=1)
+        Umean6b = numpy.nanmean(strain6b['U'].reshape(dims6b)[2:-2, 2:-2, 2:-2], axis=(0, 1, 2))
+        self.assertTrue(numpy.allclose(Umean6b,numpy.eye(3), atol=0.02))
 
         # dev and vol both ~0
-        self.assertAlmostEqual(numpy.nanmean(strain['vol']),   0, places=1)
-        self.assertAlmostEqual(numpy.nanmean(strain['volss']), 0, places=1)
-        self.assertAlmostEqual(numpy.nanmean(strain['dev']),   0, places=1)
-        self.assertAlmostEqual(numpy.nanmean(strain['devss']), 0, places=1)
+        self.assertTrue(numpy.isclose(numpy.nanmean(strain6b['vol'].reshape(  strain6b['fieldDims'])[2:-2, 2:-2, 2:-2]), 0, atol=0.02))
+        self.assertTrue(numpy.isclose(numpy.nanmean(strain6b['volss'].reshape(strain6b['fieldDims'])[2:-2, 2:-2, 2:-2]), 0, atol=0.02))
+        self.assertTrue(numpy.isclose(numpy.nanmean(strain6b['dev'].reshape(  strain6b['fieldDims'])[2:-2, 2:-2, 2:-2]), 0, atol=0.02))
+        self.assertTrue(numpy.isclose(numpy.nanmean(strain6b['devss'].reshape(strain6b['fieldDims'])[2:-2, 2:-2, 2:-2]), 0, atol=0.02))
 
-        VTK = spam.helpers.readStructuredVTK("snow-ref-snow-def-strain-Q8.vtk")
-        for component in ["U", "e", "vol", "volss", "dev", "devss"]:
-            self.assertAlmostEqual(numpy.nanmean(strain[component]), numpy.nanmean(VTK[component]), places=2)
-
-        #######################################################
-        ### Back to ldic for 2D test
-        #######################################################
-        ### overwrite with 2D slices
-        tifffile.imsave(testFolder + "snow-ref.tif", snowRef[50])
-        # Only translation
-        Phi = spam.deformation.computePhi({'t': refTranslation})
-        snowDef = spam.DIC.applyPhi(snowRef, Phi=Phi)
-
-        tifffile.imsave(testFolder + "snow-def.tif", snowDef[50])
-        exitCode = subprocess.call(["spam-ldic", "-ps", "on", "-reg", "-regbb", "4", "-regbe", "1", "-regm", "0.1", "-hws", "5", '-ns', '5', testFolder + "snow-ref.tif", testFolder + "snow-def.tif", "-od", testFolder + "", "-tsv"])
-        self.assertEqual(exitCode, 0)
-        TSV = spam.helpers.readCorrelationTSV(testFolder + "snow-ref-snow-def.tsv")
-        self.assertAlmostEqual(numpy.mean(TSV['PhiField'][:,1,-1][TSV['returnStatus']==2]), refTranslation[1], places=2)
-        self.assertAlmostEqual(numpy.mean(TSV['PhiField'][:,2,-1][TSV['returnStatus']==2]), refTranslation[2], places=2)
+        #VTK6b = spam.helpers.readStructuredVTK("snow-grid-6b-strain-Q8.vtk")
+        #for component in ["U", "e", "vol", "volss", "dev", "devss"]:
+            #self.assertTrue(numpy.nanmean(strain6b[component]), numpy.nanmean(VTK6b[component]), places=2)
 
     def _test_ITKwatershed(self):
         # make sure it runs the help without error
@@ -516,7 +592,7 @@ class testAll(unittest.TestCase):
 
         self.assertEqual(numpy.mean(imDefScript[10:-10, 10:-10, 10:-10]) > numpy.mean(imDefScriptCGS[10:-10, 10:-10, 10:-10]), True)
 
-    def test_discreteDICtools(self):
+    def _test_discreteDICtools(self):
         #######################################################
         ### 1. Generate Data
         #######################################################
@@ -991,7 +1067,8 @@ class testAll(unittest.TestCase):
         exitCode = subprocess.call(["spam-ldic", "--help"], stdout=FNULL, stderr=FNULL)
         self.assertEqual(exitCode, 0)
 
-
+        exitCode = subprocess.call(["spam-regularStrain", "--help"], stdout=FNULL, stderr=FNULL)
+        self.assertEqual(exitCode, 0)
 
 
 if __name__ == '__main__':

tools/DIC/grid.py

@@ -86,8 +86,7 @@ def makeGrid(imageSize, nodeSpacing):
     return nodePositions, nodesDim
 
 
-# WARNING ------------------------- VVVVVVVVVVV--- gets easily overwritten, pass a .copy()!
-def getImagettes(im1, boundingBox, Phi, im2, searchRange, im1mask=None, im2mask=None, minMaskCoverage=0.0, greyThreshold=[-numpy.inf, numpy.inf], applyF='all'):
+def getImagettes(im1, nodePosition, halfWindowSize, Phi, im2, searchRange, im1mask=None, im2mask=None, minMaskCoverage=0.0, greyThreshold=[-numpy.inf, numpy.inf], applyF='all', twoD=False):
     """
     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
@@ -98,10 +97,14 @@ def getImagettes(im1, boundingBox, Phi, im2, searchRange, im1mask=None, im2mask=
         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
+        nodePosition : 3-component numpy array of ints
+            This defines the centre of the window to extract from im1.
+            Note: for 2D Z = 0
+
+        halfWindowSize : 3-component numpy array of ints
+            This defines the half-size of the correlation window,
+            i.e., how many pixels to extract in Z, Y, X either side of the centre.
+            Note: for 2D Z = 0
 
         Phi : 4x4 numpy array of floats
             Phi matrix representing the movement of imagette1,
@@ -142,6 +145,9 @@ def getImagettes(im1, boundingBox, Phi, im2, searchRange, im1mask=None, im2mask=
             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
 
+        twoD : bool, optional
+            Are the images two-dimensional?
+
     Returns
     -------
         Dictionary :
@@ -157,7 +163,7 @@ def getImagettes(im1, boundingBox, Phi, im2, searchRange, im1mask=None, im2mask=
             '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
@@ -171,116 +177,116 @@ def getImagettes(im1, boundingBox, Phi, im2, searchRange, im1mask=None, im2mask=
     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
+    ## Detect 2D images
+    #if im1.shape[0] == 1:
+        #twoD = True
+        ## Impose no funny business in z if in twoD
+        #halfWindowSize[0] = 0
+        #searchRange[0:2] = 0
+    #else:
+        #twoD = False
+
+    PhiNoDisp = Phi.copy()
+    PhiNoDisp[0:3,-1] -= initialDisplacement
+    if applyF == 'rigid':
+        PhiNoDisp = spam.deformation.computeRigidPhi(PhiNoDisp)
+
+    # 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)) or applyF == 'no':
+        applyPhiPad = (0, 0, 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
-        if applyF == 'rigid':
-            PhiNoDisp = spam.deformation.computeRigidPhi(PhiNoDisp)
-            #PhiNoDisp[0:3,-1] = 0.0
+        # 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.5*numpy.ceil(max(halfWindowSize)))
 
-        # 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)) 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
-            applyPhiPad = int(0.5*numpy.ceil(max(boundingBox[1::2]-boundingBox[0::2])))
-
-            if twoD: applyPhiPad = (          0, applyPhiPad, applyPhiPad)
-            else:    applyPhiPad = (applyPhiPad, applyPhiPad, applyPhiPad)
+        if twoD: applyPhiPad = (     0     , applyPhiPad, applyPhiPad)
+        else:    applyPhiPad = (applyPhiPad, applyPhiPad, applyPhiPad)
 
-        startStopIm1 = [int(boundingBox[0] - applyPhiPad[0]), int(boundingBox[1] + applyPhiPad[0] + 1),
-                        int(boundingBox[2] - applyPhiPad[1]), int(boundingBox[3] + applyPhiPad[1] + 1),
-                        int(boundingBox[4] - applyPhiPad[2]), int(boundingBox[5] + applyPhiPad[2] + 1)]
+    # 2020-09-25 OS and EA: Prepare startStop array for imagette 1 to be extracted with new slicePadded
+    startStopIm1 = [int(nodePosition[0] - halfWindowSize[0] - applyPhiPad[0]    ),
+                    int(nodePosition[0] + halfWindowSize[0] + applyPhiPad[0] + 1),
+                    int(nodePosition[1] - halfWindowSize[1] - applyPhiPad[1]    ),
+                    int(nodePosition[1] + halfWindowSize[1] + applyPhiPad[1] + 1),
+                    int(nodePosition[2] - halfWindowSize[2] - applyPhiPad[2]    ),
+                    int(nodePosition[2] + halfWindowSize[2] + applyPhiPad[2] + 1)]
 
-        # In either case, extract imagette1, now guaranteed to be the right size
-        imagette1padded = spam.helpers.slicePadded(im1, startStopIm1)
+    # 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)) or applyF == 'no':
-            # In this case there is is no padding (despite the name) and we can just keep going
-            imagette1def = imagette1padded
+    # If F is not the identity, apply F and undo crop
+    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)
+
+        # undo padding
+        if twoD:
+            imagette1def = imagette1paddedDef[              :               ,
+                                              applyPhiPad[1]:-applyPhiPad[1],
+                                              applyPhiPad[2]:-applyPhiPad[2]]
         else:
-            # apply F to imagette 1 padded
-            if twoD:    imagette1paddedDef = spam.DIC.applyPhiPython(imagette1padded, PhiNoDisp)
-            else:       imagette1paddedDef = spam.DIC.applyPhi(      imagette1padded, PhiNoDisp)
-
-            # undo padding
-            if twoD:
-                imagette1def = imagette1paddedDef[              :               ,
-                                                  applyPhiPad[1]:-applyPhiPad[1],
-                                                  applyPhiPad[2]:-applyPhiPad[2]]
-            else:
-                imagette1def = imagette1paddedDef[applyPhiPad[0]:-applyPhiPad[0],
-                                                  applyPhiPad[1]:-applyPhiPad[1],
-                                                  applyPhiPad[2]:-applyPhiPad[2]]
-
-        ### 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 + numpy.array([0, 1, 0, 1, 0, 1])) != 0
-            maskVolumeCondition = (imagette1mask != 0).mean() >= minMaskCoverage
-
-
-        # Make sure imagette is not 0-dimensional in any dimension
-        # Check minMaskVolume
-        if numpy.all(numpy.array(imagette1def.shape) > 0) or (twoD and numpy.all(numpy.array(imagette1def.shape[1:3]) > 0)):
-            # ------------ 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] + initialDisplacement[0] + searchRange[0]    ),
-                                    int(boundingBox[1] + initialDisplacement[0] + searchRange[1] + 1),
-                                    int(boundingBox[2] + initialDisplacement[1] + searchRange[2]    ),
-                                    int(boundingBox[3] + initialDisplacement[1] + searchRange[3] + 1),
-                                    int(boundingBox[4] + initialDisplacement[2] + searchRange[4]    ),
-                                    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)
-
-                # Failed minMaskVolume condition
-                else:
-                    returnStatus = -5
-                    imagette1def = None
-                    imagette2 = None
-
-            # Failed greylevel condition
+            imagette1def = imagette1paddedDef[applyPhiPad[0]:-applyPhiPad[0],
+                                              applyPhiPad[1]:-applyPhiPad[1],
+                                              applyPhiPad[2]:-applyPhiPad[2]]
+
+    ### Check mask
+    if im1mask is None:
+        # no mask1 --> always pas this test (e.g., labelled image)
+        maskVolumeCondition = True
+        imagette1mask = None
+    else:
+        startStopIm1 = [int(nodePosition[0] - halfWindowSize[0]    ),
+                        int(nodePosition[0] + halfWindowSize[0] + 1),
+                        int(nodePosition[1] - halfWindowSize[1]    ),
+                        int(nodePosition[1] + halfWindowSize[1] + 1),
+                        int(nodePosition[2] - halfWindowSize[2]    ),
+                        int(nodePosition[2] + halfWindowSize[2] + 1)]
+        imagette1mask = spam.helpers.slicePadded(im1mask, startStopIm1) != 0
+        maskVolumeCondition = (imagette1mask != 0).mean() >= minMaskCoverage
+
+
+    # Make sure imagette is not 0-dimensional in any dimension
+    # Check minMaskVolume
+    if numpy.all(numpy.array(imagette1def.shape) > 0) or (twoD and numpy.all(numpy.array(imagette1def.shape[1:3]) > 0)):
+        # ------------ 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(nodePosition[0] - halfWindowSize[0] + initialDisplacement[0] + searchRange[0]    ),
+                                int(nodePosition[0] + halfWindowSize[0] + initialDisplacement[0] + searchRange[1] + 1),
+                                int(nodePosition[1] - halfWindowSize[1] + initialDisplacement[1] + searchRange[2]    ),
+                                int(nodePosition[1] + halfWindowSize[1] + initialDisplacement[1] + searchRange[3] + 1),
+                                int(nodePosition[2] - halfWindowSize[2] + initialDisplacement[2] + searchRange[4]    ),
+                                int(nodePosition[2] + halfWindowSize[2] + initialDisplacement[2] + searchRange[5] + 1)]
+                imagette2 = spam.helpers.slicePadded(im2, startStopIm2)
+
+                if im2mask is not None:
+                    imagette2mask = spam.helpers.slicePadded(im2mask, startStopIm2)
+
+            # Failed minMaskVolume condition
             else:
                 returnStatus = -5
                 imagette1def = None
                 imagette2 = None
 
-        # Failed 0-dimensional imagette test
+        # Failed greylevel condition
         else:
             returnStatus = -5
             imagette1def = None
             imagette2 = None
 
-    # Failed bounding box test
+    # Failed 0-dimensional imagette test
     else:
-        returnStatus = -7
+        returnStatus = -5