Resolves #206, default RS trusted in regularStrains is now 2. Filter F option...

Resolves #206, default RS trusted in regularStrains is now 2. Filter F option in filterPhiField now properly 'no'

scripts/spam-filterPhiField

@@ -215,7 +215,7 @@ if args.CINT:
         print("\tYou asked to correct with local QC fitting with -clqf, but only one correciton mode is supported")
 
 elif args.CLQF:
-    if args.FILTER_F is not 'disp':
+    if args.FILTER_F != 'no':
         print("WARNING: non-displacement quadratic coherency correction not implemented, only doing displacements, and returning F=eye(3)\n")
 
     print("\n\nCorrection based on local quadratic coherency")
@@ -271,7 +271,7 @@ if args.FILTER_MEDIAN:
             print("\t9/9")
             outputPhiField[:, 2,  2] = scipy.ndimage.generic_filter(inputPhiField[:, 2,  2].reshape(inputNodesDim), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
 
-        if args.FILTER_F == 'disp':
+        if args.FILTER_F == 'no':
             for n in range(inputNodePositions.shape[0]):
                 outputPhiField[n] = numpy.eye(4)
 

scripts/spam-regularStrain

@@ -24,8 +24,9 @@ os.environ['OMP_NUM_THREADS'] = '1'
 
 import argparse
 import tifffile
-import numpy
 import multiprocessing
+import numpy
+numpy.seterr(all='ignore')
 
 import spam.helpers
 import spam.DIC
@@ -76,21 +77,29 @@ dispFlat = f["PhiField"][:,:3,-1]
 
 # Check if a mask of points (based on return status of the correlation) is asked
 if args.MASK:
-    mask = f["returnStatus"] < -4
+    mask = f["returnStatus"] < args.RETURN_STATUS_THRESHOLD
+    print(f"\nspam-regularStrain: excluding points based on return threshold < {args.RETURN_STATUS_THRESHOLD} (excluded {100*(1-numpy.mean(mask)):2.1f}%)")
     dispFlat[mask] = numpy.nan
 
 disp = dispFlat.reshape(dims[0], dims[1], dims[2], 3)
 
 print("\nspam-regularStrain: Computing F=I+du/dx")
-if not args.Q8:
-    Ffield = spam.deformation.FfieldRegularGeers(disp,
+if args.Q8:
+    Ffield = spam.deformation.FfieldRegularQ8(   disp,
                                                  nodeSpacing=nodeSpacing,
-                                                 neighbourRadius=args.STRAIN_NEIGHBOUR_RADIUS,
                                                  nProcesses=args.PROCESSES,
                                                  verbose=True)
+elif args.RAW:
+    # Just take it straight form the file
+    Ffield = f["PhiField"][:,:3,:3]
+    if args.MASK:
+        Ffield[mask] = numpy.nan
+    Ffield = Ffield.reshape(dims[0], dims[1], dims[2], 3, 3)
+
 else:
-    Ffield = spam.deformation.FfieldRegularQ8(   disp,
+    Ffield = spam.deformation.FfieldRegularGeers(disp,
                                                  nodeSpacing=nodeSpacing,
+                                                 neighbourRadius=args.STRAIN_NEIGHBOUR_RADIUS,
                                                  nProcesses=args.PROCESSES,
                                                  verbose=True)
 
@@ -104,30 +113,28 @@ decomposedFfield = spam.deformation.decomposeFfield(Ffield,
 
 
 # Define base fileName
-if not args.Q8:
-    fileNameBase = args.OUT_DIR+"/"+args.PREFIX+"-strain-Geers"
-    mode = "Geers"
-else:
+if args.Q8:
     if twoD:
         fileNameBase = args.OUT_DIR+"/"+args.PREFIX+"-strain-Q4"
     else:
         fileNameBase = args.OUT_DIR+"/"+args.PREFIX+"-strain-Q8"
     mode = "Q8"
+elif args.RAW:
+    fileNameBase = args.OUT_DIR+"/"+args.PREFIX+"-strain-raw"
+    mode = "raw"
+else:
+    fileNameBase = args.OUT_DIR+"/"+args.PREFIX+"-strain-Geers"
+    mode = "Geers"
+
 
 # Save strain fields
 print("\nspam-regularStrain: Saving strain fields...")
 if args.TSV:
-    # Positions for Geers are the measurement points
-    if not args.Q8:
-        if twoD:
-            outputPositions = fieldCoords.copy().reshape(1, dims[1], dims[2], 3)
-        else:
-            outputPositions = fieldCoords.copy().reshape(dims[0], dims[1], dims[2], 3)
 
     # Positions for the centres of the Q8 elements are between the measurement points
     #   (so there is one number fewer compared to measurement points
     #   so we strip off last node points -- not Z ones in twoD for Q8 mode
-    else:
+    if args.Q8:
         if twoD:
             outputPositions = fieldCoords.copy().reshape(1, dims[1], dims[2], 3)[:, 0:-1, 0:-1, :]
         else:
@@ -136,6 +143,12 @@ if args.TSV:
         outputPositions[:, :, :, 0] += nodeSpacing[0] / 2.0
         outputPositions[:, :, :, 1] += nodeSpacing[1] / 2.0
         outputPositions[:, :, :, 2] += nodeSpacing[2] / 2.0
+    else:
+        # Positions for Geers and "raw" are the measurement points
+        if twoD:
+            outputPositions = fieldCoords.copy().reshape(1, dims[1], dims[2], 3)
+        else:
+            outputPositions = fieldCoords.copy().reshape(dims[0], dims[1], dims[2], 3)
 
     # Here we want to pass an Nx3 matrix of poitions:
     spam.helpers.writeStrainTSV(fileNameBase+".tsv", outputPositions.reshape(-1,3), decomposedFfield, firstColumn="StrainPointNumber")

tests/test_scripts.py

@@ -560,6 +560,28 @@ class testAll(unittest.TestCase):
         #for component in ["U", "e", "vol", "volss", "dev", "devss"]:
             #self.assertTrue(numpy.nanmean(strain6b[component]), numpy.nanmean(VTK6b[component]), places=2)
 
+        ### Step 6b Now let's calculate small and large strains with Q8 elements and TSV output
+        exitCode = subprocess.call(["spam-regularStrain", "-raw",
+                                    testFolder + "snow-grid-5a-ldic.tsv",
+                                    "-comp", "U", "e", "vol", "volss", "dev", "devss",
+                                    "-vtk",
+                                    "-od", testFolder, "-pre", "snow-grid-6b"],
+                                   stdout=FNULL, stderr=FNULL)
+        self.assertEqual(exitCode, 0)
+        strain6b = spam.helpers.readStrainTSV(testFolder + "snow-grid-6b-strain-raw.tsv")
+        dims6b = [strain6b['fieldDims'][0], strain6b['fieldDims'][1], strain6b['fieldDims'][2], 3, 3]
+        # Check small strains, should be ~0
+        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))
+
+        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.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))
+
     def test_ITKwatershed(self):
         # Run on snow data
         # load 3D image

tools/helpers/optionsParser.py

@@ -1556,7 +1556,14 @@ def regularStrainParser(parser):
                         '--nomask',
                         action="store_false",
                         dest='MASK',
-                        help="Don't mask correlation points in background according to return status (i.e., include RS=-5 or less)")
+                        help="Don't mask correlation points according to return status (use everything)")
+
+    parser.add_argument('-rst',
+                        '--return-status-threshold',
+                        type=int,
+                        default=2,
+                        dest='RETURN_STATUS_THRESHOLD',
+                        help='Lowest return status value to preserve in input PhiField. Default = 2')
 
     parser.add_argument('-r',
                         '--neighbourhood-radius-for-strain-calculation',
@@ -1566,12 +1573,19 @@ def regularStrainParser(parser):
                         help="Radius (in units of nodeSpacing) inside which to select neighbours for displacement gradient calculation. Ignored if -cub is set. Default = 1.5")
 
     parser.add_argument('-cub',
+                        '-Q8',
                         '--cubic-element',
                         '--Q8',
                         action="store_true",
                         dest='Q8',
                         help='Use Q8 element interpolation? More noisy and strain values not centred on displacement points')
 
+    parser.add_argument('-raw',
+                        '--no-shape-function',
+                        action="store_true",
+                        dest='RAW',
+                        help='Just use F straight from the correlation windows instead of computing it from the displacement field.')
+
     parser.add_argument('-od',
                         '--out-dir',
                         type=str,
@@ -1579,6 +1593,8 @@ def regularStrainParser(parser):
                         dest='OUT_DIR',
                         help='Output directory, default is the dirname of input file')
 
+
+
     parser.add_argument('-pre',
                         '--prefix',
                         type=str,
@@ -1632,6 +1648,10 @@ def regularStrainParser(parser):
     if args.PREFIX is None:
         args.PREFIX = os.path.splitext(os.path.basename(args.inFile.name))[0]
 
+    if args.RAW and args.Q8:
+        print("You can't ask for both F-from-correlation and F-from-Q8")
+        exit()
+
     # Make sure at least one output format has been asked for
     if args.VTK + args.TIFF + args.TSV == 0:
         print("#############################################################")
@@ -2832,6 +2852,10 @@ def filterPhiField(parser):
         print("WARINING: you can't ask for an overall median filter and a correction")
         exit()
 
+    if args.FILTER_F not in ['all', 'rigid', 'no']:
+        print("-F option must be either 'all', 'rigid' or 'no'")
+        exit()
+
     # Output file name prefix
     if args.PREFIX is None:
         args.PREFIX = os.path.splitext(os.path.basename(args.PHIFILE.name))[0] + "-filtered"