Sketch of merge grid and discrete onto grid mode for spam-passPhiField

docs/source/kalisphera.rst

-kalisphera package
-==================
-
-Submodules
-----------
-
-kalisphera.kalisphera module
-----------------------------
-
-.. automodule:: kalisphera.kalisphera
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
-kalisphera.kalispheraToolkit module
------------------------------------
-
-.. automodule:: kalisphera.kalispheraToolkit
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
-kalisphera.kalispheraToolkit module
------------------------------------
-
-.. automodule:: kalisphera.kalispheraToolkit
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
-kalisphera.spheroid module
---------------------------
-
-.. automodule:: kalisphera.spheroid
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
-Module contents
----------------
-
-.. automodule:: kalisphera
-   :members:
-   :undoc-members:
-   :show-inheritance:

scripts/spam-filterPhiField

@@ -40,8 +40,8 @@ You should have received a copy of the GNU General Public License along with
 this program.  If not, see <http://www.gnu.org/licenses/>.
 """
 
-#import os
-#os.environ['OPENBLAS_NUM_THREADS'] = '1'
+import os
+os.environ['OPENBLAS_NUM_THREADS'] = '1'
 
 import spam.DIC
 import spam.deformation
@@ -52,7 +52,7 @@ import spam.label
 import numpy
 import multiprocessing
 import scipy.spatial
-#import progressbar
+import progressbar
 import argparse
 import tifffile
 
@@ -83,7 +83,7 @@ if PhiFromFile is None:
 # If the read Phi-file has only one line -- it's a single point registration!
 # We can either apply it to a grid or to labels
 if PhiFromFile['fieldCoords'].shape[0] == 1:
-    print(f"\tYour TSV passed with -pf {args.PHIFILE.name} is single line file. A field is required")
+    print(f"\tYour TSV passed with -pf {args.PHIFILE.name} is single line file (a registration). A field is required")
     exit()
 
 # Check if it is a discrete or gridded field
@@ -108,10 +108,10 @@ inputError         = PhiFromFile["error"]
 
 # output arrays
 outputPhiField      = numpy.zeros((inputNodePositions.shape[0], 4, 4))
-outputReturnStatus  = numpy.ones((inputNodePositions.shape[0]), dtype=float)
-outputDeltaPhiNorm  = numpy.ones((inputNodePositions.shape[0]), dtype=float)*100
+outputReturnStatus  = numpy.ones( (inputNodePositions.shape[0]), dtype=float)
+outputDeltaPhiNorm  = numpy.ones( (inputNodePositions.shape[0]), dtype=float)*100
 outputIterations    = numpy.zeros((inputNodePositions.shape[0]), dtype=float)
-outputError         = numpy.ones((inputNodePositions.shape[0]), dtype=float)*100
+outputError         = numpy.ones( (inputNodePositions.shape[0]), dtype=float)*100
 outputPixelSearchCC = numpy.zeros((inputNodePositions.shape[0]), dtype=float)
 
 
@@ -119,6 +119,7 @@ outputPixelSearchCC = numpy.zeros((inputNodePositions.shape[0]), dtype=float)
 if args.NEIGHBOUR_RADIUS is not None and args.NUMBER_OF_NEIGHBOURS is not None:
     print("Both number of neighbours and neighbour radius are set, I'm taking the radius and ignoring the number of neighbours")
     args.NUMBER_OF_NEIGHBOURS = None
+
 if args.NEIGHBOUR_RADIUS is None and args.NUMBER_OF_NEIGHBOURS is None:
     if grid:
         # Gridded input field

scripts/spam-reg

@@ -89,6 +89,7 @@ if args.INTERPOLATION_ORDER == 1:
 else:
     interpolator = 'python'
 
+
 regReturns = spam.DIC.correlate.registerMultiscale(im1,
                                                    im2,
                                                    args.BIN_BEGIN, binStop=args.BIN_END,

scripts/spam-regularStrain

@@ -84,14 +84,14 @@ else:
 # Check if the correction of the input field is asked
 if args.CORRECT_FIELD or args.CORRECT_MEDIAN_FILTER:
     print("\nspam-regularStrain: Correcting/filtering field...")
-    fieldValues = spam.deformation.correctPhiField(fileName           = args.inFile.name,
-                                                   correctBadPoints   = args.CORRECT_FIELD,
-                                                   ignoreBackGround   = ignoreBackGround,
-                                                   deltaPhiNormMin    = args.CORRECT_DELTA_PHI_NORM,
-                                                   pixelSearchCCmin   = args.CORRECT_PIXEL_SEARCH_CC,
-                                                   nNeighbours        = args.CORRECT_NEIGHBOURS,
-                                                   filterPoints       = args.CORRECT_MEDIAN_FILTER,
-                                                   filterPointsRadius = args.CORRECT_MEDIAN_FILTER_RADIUS)
+    fieldValues = spam.DIC.correctPhiField(fileName           = args.inFile.name,
+                                           correctBadPoints   = args.CORRECT_FIELD,
+                                           ignoreBackGround   = ignoreBackGround,
+                                           deltaPhiNormMin    = args.CORRECT_DELTA_PHI_NORM,
+                                           pixelSearchCCmin   = args.CORRECT_PIXEL_SEARCH_CC,
+                                           nNeighbours        = args.CORRECT_NEIGHBOURS,
+                                           filterPoints       = args.CORRECT_MEDIAN_FILTER,
+                                           filterPointsRadius = args.CORRECT_MEDIAN_FILTER_RADIUS)
     # Extract the corrected displacement field and reshape it
     disp = fieldValues[:,:3,-1].reshape(dims[0], dims[1], dims[2], 3)
 else:

tools/helpers/optionsParser.py

@@ -68,7 +68,7 @@ def ldicParser(parser):
                         dest='PHIFILE',
                         default=None,
                         type=argparse.FileType('r'),
-                        help="Path to TSV file containing initial F guess, can be single-point registration or multiple point correlation. Default = None")
+                        help="Path to TSV file containing initial Phi guess, can be single-point registration or multiple point correlation. Default = None")
 
     parser.add_argument('-pfb',
                         '--phiFile-bin-ratio',
@@ -408,7 +408,7 @@ def ddicParser(parser):
                         dest='PHIFILE',
                         default=None,
                         type=argparse.FileType('r'),
-                        help="Path to TSV file containing initial F guess, can be single-point registration or multiple point correlation. Default = None")
+                        help="Path to TSV file containing initial Phi guess, can be single-point registration or multiple point correlation. Default = None")
 
     parser.add_argument('-pfb',
                         '--phiFile-bin-ratio',
@@ -837,7 +837,7 @@ def gdicParser(parser):
                         dest='PHIFILE',
                         default=None,
                         type=argparse.FileType('r'),
-                        help="Path to TSV file containing initial F guess, can be single-point registration or multiple point correlation. Default = None")
+                        help="Path to TSV file containing initial Phi guess, can be single-point registration or multiple point correlation. Default = None")
 
     parser.add_argument('-pfb',
                         '--phiFile-bin-ratio',
@@ -1840,9 +1840,9 @@ def register(parser):
                         '-mar',
                         '--margin',
                         type=float,
-                        default=0.1,
+                        default=None,
                         dest='MARGIN',
-                        help='Interpolation margin in proportions of image size. Default = 0.1, which means 0.1 * image size from both sides in all directions')
+                        help='Interpolation margin in pixels. Default is the default for spam.DIC.registerMultiscale')
 
     parser.add_argument('-m3',
                         '-mar3',
@@ -1851,7 +1851,7 @@ def register(parser):
                         type=int,
                         default=None,
                         dest='MARGIN',
-                        help="Interpolation margin in pixels. Default = 0.1 * image size from both sides in all directions")
+                        help='ZYX interpolation margin in pixels. Default is the default for spam.DIC.registerMultiscale')
 
     parser.add_argument('-ug',
                         '--update-gradient',
@@ -1933,7 +1933,6 @@ def register(parser):
     else:
         args.PREFIX += "-registration"
 
-
     return args
 
 
@@ -2498,7 +2497,7 @@ def passPhiField(parser):
                         dest='PHIFILE',
                         default=None,
                         type=argparse.FileType('r'),
-                        help="Path to TSV file containing initial F guess, can be single-point registration or multiple point correlation. Default = None")
+                        help="Path to TSV file containing initial Phi guess, can be single-point registration or multiple point correlation. Default = None")
 
     parser.add_argument('-pfb',
                         '--phiFile-bin-ratio',
@@ -2524,7 +2523,8 @@ def passPhiField(parser):
     parser.add_argument('-lab1',
                         '--labelledFile1',
                         dest='LAB1',
-                        default=None,
+                        nargs='+',
+                        default=[],
                         type=argparse.FileType('r'),
                         help="Path to tiff file containing a labelled image 1 that defines zones to correlate. Disactivates -hws and -ns options")
 
@@ -2560,6 +2560,28 @@ def passPhiField(parser):
                         dest='im1shape',
                         help="Size of im1 in pixels Z Y X")
 
+    parser.add_argument('-pf2',
+                        '-phiFile2',
+                        dest='PHIFILE2',
+                        nargs='+',
+                        default=[],
+                        type=argparse.FileType('r'),
+                        help="Path to second spam-ddic TSV file(s). Default = None")
+
+    parser.add_argument('-pf2b',
+                        '--phiFile2-bin-ratio',
+                        type=int,
+                        default=1,
+                        dest='PHIFILE2_BIN_RATIO',
+                        help="Ratio of binning level between second 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('-mpl',
+                        '--merge-prefer-label',
+                        action="store_true",
+                        dest='MERGE_PREFER_LABEL',
+                        help='When merging grid and discrete correlation results, automatically prefer points inside labels? Default = False')
+
+
     parser.add_argument('-nr',
                         '--neighbourhood-radius-px',
                         type=float,
@@ -2641,49 +2663,56 @@ def passPhiField(parser):
         args.PREFIX += "-passed"
 
 
-    if args.LAB1 is not None: 
-        # We have a labelled image and so no nodeSpacing or halfWindowSize
-        print("I have been passed a labelled image and so I am disactivating:")
-        print("\t- node spacing")
-        args.NS            = None
-        args.im1           = None
-        args.im1shape      = None
-        # Output file name prefix
-        args.PREFIX += "-labelled"
+    if len(args.PHIFILE2) > 0:
+        print("\n\nMerge mode")
+        args.PREFIX += "-merged"
+        print()
 
     else:
-        print("No labelled image so I'm in grid mode")
-        # We are in grid, with a nodeSpacing and halfWindowSize
-        # Catch interdependent node spacing and correlation window sizes
-        if args.NS is None:
-            print("...actually no node spacing either so, output basis not defined!")
-            exit()
+        if len(args.LAB1) == 1: 
+            # We have a labelled image and so no nodeSpacing or halfWindowSize
+            print("\n\nI have been passed a labelled image and so I am disactivating:")
+            print("\t- node spacing")
+            args.NS            = None
+            args.im1           = None
+            args.im1shape      = None
+            # Output file name prefix
+            args.PREFIX += "-labelled"
+            print()
+
         else:
-            # Catch 3D options
-            if len(args.NS) == 1:
-                args.PREFIX += f"-ns{args.NS[0]}"
-                args.NS = [args.NS[0], args.NS[0], args.NS[0]]
+            print("\n\nNo labelled image so I'm in grid mode")
+            # We are in grid, with a nodeSpacing and halfWindowSize
+            # Catch interdependent node spacing and correlation window sizes
+            if args.NS is None:
+                print("...actually no node spacing either so, output basis not defined!")
+                exit()
             else:
-                # 3 NSs are passed
-                args.PREFIX += f"-ns{args.NS[0]}-{args.NS[1]}-{args.NS[2]}"
-
-            if args.im1 is None and args.im1shape is None:
-                print("In grid mode, I need to know the image size, please pass either -im1 or -im1shape")
+                # Catch 3D options
+                if len(args.NS) == 1:
+                    args.PREFIX += f"-ns{args.NS[0]}"
+                    args.NS = [args.NS[0], args.NS[0], args.NS[0]]
+                else:
+                    # 3 NSs are passed
+                    args.PREFIX += f"-ns{args.NS[0]}-{args.NS[1]}-{args.NS[2]}"
+
+                if args.im1 is None and args.im1shape is None:
+                    print("In grid mode, I need to know the image size, please pass either -im1 or -im1shape")
+                    exit()
+            # We need some way to define the image size for output
+            if args.im1 is not None:
+                print("Getting im1 dimensions by looking in the file (this ignores -im1shape)")
+                tiff = tifffile.TiffFile(args.im1.name)
+                args.im1shape = tiff.series[0].shape
+            elif args.im1shape is not None:
+                print("Trusting -im1shape dimensions as passed")
+            else:
+                print("You asked for a node spacing, but I don't know the size of the image you want me to define the grid on! Pass -im1 im.tif or -im1shape Z Y X")
                 exit()
-        # We need some way to define the image size for output
-        if args.im1 is not None:
-            print("Getting im1 dimensions by looking in the file (this ignores -im1shape)")
-            tiff = tifffile.TiffFile(args.im1.name)
-            args.im1shape = tiff.series[0].shape
-        elif args.im1shape is not None:
-            print("Trusting -im1shape dimensions as passed")
-        else:
-            print("You asked for a node spacing, but I don't know the size of the image you want me to define the grid on! Pass -im1 im.tif or -im1shape Z Y X")
-            exit()
 
-    if args.MODE not in ["all", "rigid", "disp",]:
-        print("-m should be 'all' 'rigid' or 'disp'")
-        exit()
+        if args.MODE not in ["all", "rigid", "disp",]:
+            print("-m should be 'all' 'rigid' or 'disp'")
+            exit()
 
     return args
 
@@ -2694,7 +2723,7 @@ def filterPhiField(parser):
                         dest='PHIFILE',
                         default=None,
                         type=argparse.FileType('r'),
-                        help="Path to TSV file containing initial F guess, can be single-point registration or multiple point correlation. Default = None")
+                        help="Path to TSV file containing initial Phi guess, can be single-point registration or multiple point correlation. Default = None")
 
     parser.add_argument('-pfb',
                         '--phiFile-bin-ratio',

tools/tests/test_DIC.py

@@ -6,6 +6,14 @@ import numpy
 import spam.deformation
 import spam.DIC
 import spam.datasets
+import spam.kalisphera
+import spam.label
+import scipy.ndimage
+import tifffile
+import subprocess
+
+testFolder = './'
+
 
 numpy.set_printoptions(precision=3)
 
@@ -56,6 +64,14 @@ class TestFunctionDVC(unittest.TestCase):
             os.remove("spamPhiFieldCF-corrected-N12-filteredRad3.tsv")
             os.remove("spamPhiFieldCF-ignoreBadPoints.tsv")
             os.remove("spamPhiFieldCFDel.tsv")
+            os.remove(testFolder+'Step0.tif')
+            os.remove(testFolder+'Lab0.tif')
+            os.remove(testFolder+'Step1.tif')
+            os.remove(testFolder+'Step0-Step1-registration.tsv')
+            os.remove(testFolder+'Step0-Step1.tsv')
+            os.remove(testFolder+'Step0-Step1-ddic.tsv')
+            os.remove(testFolder+'Step0-Step1-ddic.vtk')
+            #os.remove(testFolder+'merged.tsv')
         except OSError:
             pass
 
@@ -90,20 +106,20 @@ class TestFunctionDVC(unittest.TestCase):
 
     def test_correctPhiField(self):
         # case 0: wrong tsv name
-        f0 = spam.deformation.correctPhiField(fileName="./sapmPhiFieldCF.tsv")
+        f0 = spam.DIC.correctPhiField(fileName="./sapmPhiFieldCF.tsv")
 
         # case 0b: don't pass all separate arrays
-        f0b = spam.deformation.correctPhiField(fileName=None,
+        f0b = spam.DIC.correctPhiField(fileName=None,
                                      fieldCoords=fieldCoords,
                                      fieldRS=RS.copy(), fieldDPhi=DPhi.copy(),
                                      ignoreBadPoints=True)
 
         # case 1: read a tsv file and ignore the bad points
-        f1 = spam.deformation.correctPhiField(fileName="./spamPhiFieldCF.tsv", ignoreBadPoints=True, saveFile=True)
+        f1 = spam.DIC.correctPhiField(fileName="./spamPhiFieldCF.tsv", ignoreBadPoints=True, saveFile=True)
         self.assertEqual(PhiField[:, 0:3, -1].sum() - 3, numpy.nansum(f1[:, 0:3, -1]))
 
         # case 1b: read separate arrays and ignore the bad points based on subPixel
-        f1b = spam.deformation.correctPhiField(fileName=None,
+        f1b = spam.DIC.correctPhiField(fileName=None,
                                      fieldCoords=fieldCoords, fieldValues=PhiField.copy(),
                                      fieldRS=RS.copy(), fieldDPhi=DPhi.copy(),
                                      fieldPixelSearchCC=numpy.zeros_like(RS), fieldIT=IT,
@@ -111,7 +127,7 @@ class TestFunctionDVC(unittest.TestCase):
         self.assertEqual(PhiField[:, 0:3, -1].sum() - 3, numpy.nansum(f1b[:, 0:3, -1]))
 
         # case 1c: read separate arrays and ignore the bad points based on pixelSearch
-        f1c = spam.deformation.correctPhiField(fileName=None,
+        f1c = spam.DIC.correctPhiField(fileName=None,
                                      fieldCoords=fieldCoords, fieldValues=PhiField,
                                      fieldRS=RS, fieldDPhi=numpy.full((fieldCoords.shape[0]), 100),
                                      fieldPixelSearchCC=numpy.ones_like(RS), fieldIT=IT,
@@ -119,22 +135,22 @@ class TestFunctionDVC(unittest.TestCase):
         self.assertEqual(PhiField[:, 0:3, -1].sum(), f1c[:, 0:3, -1].sum())
 
         # case 2: read a tsv file correct, filter adn save it
-        f2 = spam.deformation.correctPhiField(fileName="./spamPhiFieldCF.tsv", correctBadPoints=True, filterPoints=True, saveFile=True)
+        f2 = spam.DIC.correctPhiField(fileName="./spamPhiFieldCF.tsv", correctBadPoints=True, filterPoints=True, saveFile=True)
         self.assertEqual(f2[:, 0:3, -1].sum(), f2.shape[0] * 3)
 
         # case 2b: read a tsv file and correct it with the closest neighbour
-        f2b = spam.deformation.correctPhiField(fileName="./spamPhiFieldCF.tsv", correctBadPoints=True, nNeighbours=1)
+        f2b = spam.DIC.correctPhiField(fileName="./spamPhiFieldCF.tsv", correctBadPoints=True, nNeighbours=1)
         self.assertEqual(f2b[:, 0:3, -1].sum(), fieldCoords.shape[0] * 3)
 
         # case 3: read a tsv file, correct it ignoring the background and save it
-        f3 = spam.deformation.correctPhiField(fileName="./spamPhiFieldCF.tsv", correctBadPoints=True, ignoreBackGround=True, saveFileName="./spamPhiFieldCFDel.tsv")
+        f3 = spam.DIC.correctPhiField(fileName="./spamPhiFieldCF.tsv", correctBadPoints=True, ignoreBackGround=True, saveFileName="./spamPhiFieldCFDel.tsv")
         self.assertEqual(numpy.nansum(f3[:, 0:3, -1]), (fieldCoords.shape[0] - 2) * 3)
 
         # case 3b: read separate arrays, correct the field ignoring the background and save it
         PSCC = numpy.ones_like(DPhi)
         for badPoint in range(1, PSCC.shape[0], 3):
             PSCC[badPoint] = 0.1
-        f3b = spam.deformation.correctPhiField(fieldCoords=fieldCoords, fieldValues=PhiField.copy(),
+        f3b = spam.DIC.correctPhiField(fieldCoords=fieldCoords, fieldValues=PhiField.copy(),
                                      fieldRS=RS, fieldDPhi=numpy.full((fieldCoords.shape[0]), 100),
                                      fieldPixelSearchCC=PSCC, fieldIT=IT,
                                      correctBadPoints=True, ignoreBackGround=True,
@@ -488,6 +504,137 @@ class TestFunctionDVC(unittest.TestCase):
         self.assertEqual(imagetteReturns7["returnStatus"], -5)
         self.assertTrue(imagetteReturns7["imagette1"] is None)
 
+    def test_merge(self):
+        #######################################################
+        ### We're using the DDIC test from scripts here, lightly modified
+        #######################################################
+        #First we need to create some data using DEM dataset
+        pixelSize = 0.0001
+        blurSTD = 0.8
+        noiseSTD = 0.01
+        boxSizeDEM, centres, radii = spam.datasets.loadUniformDEMboxsizeCentreRadius()
+
+        # put 0 in the middle
+        centres -= numpy.mean(centres, axis=0)
+        rMax = numpy.amax(radii)
+
+        # pad box size
+        boxSizeDEM = boxSizeDEM + 5 * rMax
+
+        # add half box size to centres
+        centres += numpy.array(boxSizeDEM)/2.0
+        boxSizePx = (boxSizeDEM / pixelSize).astype(int)
+        centresPx = centres / pixelSize
+        radiiPx = radii / pixelSize
+        box = numpy.zeros(boxSizePx, dtype="<f8")
+        spam.kalisphera.makeSphere(box, centresPx, radiiPx)
+        box[numpy.where(box > 1.0)] = 1.0
+        box[numpy.where(box < 0.0)] = 0.0
+        box = box * 0.5
+        box = box + 0.25
+        box = scipy.ndimage.filters.gaussian_filter(box, sigma=blurSTD)
+        box = numpy.random.normal(box, scale=noiseSTD)
+        binIm0 = box >= 0.5
+        #Run watershed
+        labIm0 = spam.label.ITKwatershed.watershed(binIm0)
+        #Save images 
+        tifffile.imsave(testFolder + "Step0.tif", box.astype('<f4'))
+        tifffile.imsave(testFolder + "Lab0.tif", labIm0.astype(spam.label.labelType))
+
+        #test of rigid translation and rotation
+        #Create Phi and Apply (25 px displacement on Y-axis, and 5 degree rotation along Z axis)
+        translationStep1 = [5, 2, 0]
+        rotationStep1 = [0, 0, 0]
+        transformation = {'t': translationStep1,
+                          'r': rotationStep1}
+        Phi = spam.deformation.computePhi(transformation)
+
+        # transform centres around the centres of the box
+        centresPxDeformed = numpy.zeros_like(centresPx)
+        for i, centrePx in enumerate(centresPx):
+            centresPxDeformed[i] = centrePx + spam.deformation.decomposePhi(Phi, PhiPoint=centrePx, PhiCentre=numpy.array(boxSizePx)/2.0)['t']
+        boxDeformed = numpy.zeros(boxSizePx, dtype="<f8")
+        spam.kalisphera.makeSphere(boxDeformed, centresPxDeformed, radiiPx)
+        boxDeformed[numpy.where(boxDeformed > 1.0)] = 1.0
+        boxDeformed[numpy.where(boxDeformed < 0.0)] = 0.0
+        boxDeformed = boxDeformed * 0.5
+        boxDeformed = boxDeformed + 0.25
+        boxDeformed = scipy.ndimage.filters.gaussian_filter(boxDeformed, sigma=blurSTD)
+        boxDeformed = numpy.random.normal(boxDeformed, scale=noiseSTD)
+        #Save images
+        tifffile.imsave(testFolder + "Step1.tif", boxDeformed.astype('<f4'))
+
+
+        #######################################################
+        ### Now use the ddic and ldic 
+        #######################################################
+        exitCode = subprocess.call(["spam-reg",
+                                    testFolder + "Step0.tif",
+                                    testFolder + "Step1.tif",
+                                    "-bb", "4", "-be", "2",
+                                    "-od", testFolder+""])
+        self.assertEqual(exitCode, 0)
+
+        exitCode = subprocess.call(["spam-ldic",
+                                    "-pf", testFolder + "Step0-Step1-registration.tsv",
+                                    "-glt", "0.5",
+                                    "-hws", "10",
+                                    "-ns", "10",
+                                    testFolder + "Step0.tif", testFolder + "Step1.tif",
+                                    "-od", testFolder+""])
+        self.assertEqual(exitCode, 0)
+
+        exitCode = subprocess.call(["spam-ddic",
+                                    "-pf", testFolder + "Step0-Step1-registration.tsv",
+                                    "-ld", "2",
+                                    testFolder + "Step0.tif", testFolder + "Lab0.tif", testFolder + "Step1.tif",
+                                    "-od", testFolder+""])
+        self.assertEqual(exitCode, 0)
+
+
+        # With lists
+        spam.DIC.mergeRegularGridAndDiscrete(regularGrid=spam.helpers.readCorrelationTSV(testFolder+'Step0-Step1-ldic.tsv'),
+                                                              discrete=[spam.helpers.readCorrelationTSV(testFolder+'Step0-Step1-ddic.tsv')],
+                                                              labelledImage=[tifffile.imread(testFolder+'Lab0.tif')],
+                                                              alwaysLabel=True)
+
+        # Non-working options, either one list and not the other, or different sizes, catch None returns
+        r = spam.DIC.mergeRegularGridAndDiscrete(regularGrid=spam.helpers.readCorrelationTSV(testFolder+'Step0-Step1-ldic.tsv'),
+                                                              discrete=[spam.helpers.readCorrelationTSV(testFolder+'Step0-Step1-ddic.tsv')],
+                                                              labelledImage=[tifffile.imread(testFolder+'Lab0.tif'), None],
+                                                              alwaysLabel=False)
+        self.assertEqual(r is None, True)
+
+        r = spam.DIC.mergeRegularGridAndDiscrete(regularGrid=spam.helpers.readCorrelationTSV(testFolder+'Step0-Step1-ldic.tsv'),
+                                                              discrete=[spam.helpers.readCorrelationTSV(testFolder+'Step0-Step1-ddic.tsv')],
+                                                              labelledImage=tifffile.imread(testFolder+'Lab0.tif'),
+                                                              alwaysLabel=False)
+        self.assertEqual(r is None, True)
+
+        r = spam.DIC.mergeRegularGridAndDiscrete(regularGrid=spam.helpers.readCorrelationTSV(testFolder+'Step0-Step1-ldic.tsv'),
+                                                              discrete=[spam.helpers.readCorrelationTSV(testFolder+'Step0-Step1-ddic.tsv')],
+                                                              labelledImage=None,
+                                                              alwaysLabel=False)
+        self.assertEqual(r is None, True)
+
+        # This works again, no lists
+        output = spam.DIC.mergeRegularGridAndDiscrete(regularGrid=spam.helpers.readCorrelationTSV(testFolder+'Step0-Step1-ldic.tsv'),
+                                                              discrete=spam.helpers.readCorrelationTSV(testFolder+'Step0-Step1-ddic.tsv'),
+                                                              labelledImage=tifffile.imread(testFolder+'Lab0.tif'))
+
+        for i in range(3):
+            self.assertAlmostEqual(numpy.mean(output['PhiField'][:,i,-1][output['returnStatus']==2]) - translationStep1[i], 0, places=1)
+
+        # Check that there are more RS==2 points in the merged than in the ldic
+        rsTwoLdic = numpy.sum(spam.helpers.readCorrelationTSV(testFolder+'Step0-Step1-ldic.tsv')['returnStatus']==2)
+        self.assertEqual(numpy.sum(output['returnStatus']==2) > rsTwoLdic, True)
+
+        ## Check that there some mergeSource = 1
+        #mergeSource = numpy.unique(numpy.genfromtxt(testFolder+'merged.tsv', names=True)['mergeSource'])
+        #self.assertEqual(mergeSource[0], 0)
+        #self.assertEqual(mergeSource[1], 1)
+
+
 
     def test_LQC(self):
         points, _ = spam.DIC.grid.makeGrid((10, 10, 10), 2)

tools/tests/test_deformation.py

@@ -21,14 +21,14 @@ class testAll(unittest.TestCase):
 
     def tearDown(self):
         try:
-            os.remove(testFolder+'Step0.tif')
-            os.remove(testFolder+'Lab0.tif')
-            os.remove(testFolder+'Step1.tif')
-            os.remove(testFolder+'Step0-Step1-registration.tsv')
-            os.remove(testFolder+'Step0-Step1.tsv')
-            os.remove(testFolder+'Step0-Step1-ddic.tsv')
-            os.remove(testFolder+'Step0-Step1-ddic.vtk')
-            os.remove(testFolder+'merged.tsv')
+            #os.remove(testFolder+'Step0.tif')
+            #os.remove(testFolder+'Lab0.tif')
+            #os.remove(testFolder+'Step1.tif')
+            #os.remove(testFolder+'Step0-Step1-registration.tsv')
+            #os.remove(testFolder+'Step0-Step1.tsv')
+            #os.remove(testFolder+'Step0-Step1-ddic.tsv')
+            #os.remove(testFolder+'Step0-Step1-ddic.vtk')
+            #os.remove(testFolder+'merged.tsv')
             os.remove(testFolder+'TSV_getDisplacementFromNeighbours.tsv')
             os.remove(testFolder+'TSV_mergeRegistrationAndDiscreteFields.tsv')
 
@@ -734,142 +734,6 @@ class testAll(unittest.TestCase):
         self.assertTrue(numpy.sum(decomposedF['dev'][tetsWithPointToMove] != 0) == len(tetsWithPointToMove))
 
 
-    def _test_merge(self):
-        #######################################################
-        ### We're using the DDIC test from scripts here, lightly modified
-        #######################################################
-        #First we need to create some data using DEM dataset
-        pixelSize = 0.0001
-        blurSTD = 0.8
-        noiseSTD = 0.01
-        boxSizeDEM, centres, radii = spam.datasets.loadUniformDEMboxsizeCentreRadius()