sketch of spam-filterPhiField

scripts/spam-filterPhiField

+#!/usr/bin/env python
+
+"""
+This script manipulates a Phi field:
+
+Gridded Phi values:
+  - spam-pixelSearch
+  - spam-pixelSearchPropagate
+  - spam-ldic
+
+Phis defined at labels centres:
+  - spam-pixelSearch
+  - spam-pixelSearchPropagate
+  - spam-ddic
+
+
+This script allows you to:
+  - correct bad points inside a PhiField based on RS, or CC
+  - correct incoherent points inside a PhiField based on LQC
+  - apply a median filter to the PhiField
+
+Outputs are:
+  - TSV files
+  - (optional) VTK files for visualisation
+  - (optional) TIF files in the case of gridded data
+
+Copyright (C) 2020 SPAM Contributors
+
+This program is free software: you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the Free
+Software Foundation, either version 3 of the License, or (at your option)
+any later version.
+
+This program is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
+more details.
+
+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 spam.DIC
+import spam.deformation
+import spam.helpers
+#import spam.mesh
+import spam.label
+
+import numpy
+import multiprocessing
+import scipy.spatial
+#import progressbar
+import argparse
+import tifffile
+
+tol = 1e-6
+
+# Define argument parser object
+parser = argparse.ArgumentParser(description="spam-filterPhiField "+spam.helpers.optionsParser.GLPv3descriptionHeader +\
+                                             "This script process Phi fields by\n"+
+                                             "correcting bad or incoherent points or filtering",
+                                 formatter_class=argparse.RawTextHelpFormatter)
+
+# Parse arguments with external helper function
+args = spam.helpers.optionsParser.filterPhiField(parser)
+
+print("spam-filterPhiField -- Current Settings:")
+argsDict = vars(args)
+for key in sorted(argsDict):
+    print("\t{}: {}".format(key, argsDict[key]))
+
+###############################################################
+### Step 1 (mandatory) read input Phi File
+###############################################################
+PhiFromFile = spam.helpers.readCorrelationTSV(args.PHIFILE.name, readConvergence=True, readPixelSearchCC=True, readError=True)
+if PhiFromFile is None:
+    print(f"\tFailed to read your TSV file passed with -pf {args.PHIFILE.name}")
+    exit()
+
+# 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")
+    exit()
+
+# Check if it is a discrete or gridded field
+grid = True
+discrete = False
+if PhiFromFile["numberOfLabels"] != 0:
+    discrete = True
+    grid = False
+
+###############################################################
+### Input Phi file is a Phi FIELD
+###############################################################
+inputNodesDim      = PhiFromFile["fieldDims"]
+inputNodePositions = PhiFromFile["fieldCoords"]
+inputPhiField      = PhiFromFile["PhiField"]
+inputDisplacements = PhiFromFile["PhiField"][:, 0:3, -1]
+inputReturnStatus  = PhiFromFile["returnStatus"]
+inputPixelSearchCC = PhiFromFile["pixelSearchCC"]
+inputDeltaPhiNorm  = PhiFromFile["deltaPhiNorm"]
+inputIterations    = PhiFromFile["iterations"]
+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
+outputIterations    = numpy.zeros((inputNodePositions.shape[0]), dtype=float)
+outputError         = numpy.ones((inputNodePositions.shape[0]), dtype=float)*100
+outputPixelSearchCC = numpy.zeros((inputNodePositions.shape[0]), dtype=float)
+
+
+# Check neighbour inputs, either args.NEIGHBOUR_RADIUS or args.NUMBER_OF_NEIGHBOURS should be set.
+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
+        nodeSpacing = numpy.array([numpy.unique(inputNodePositions[:, i])[1] - numpy.unique(inputNodePositions[:, i])[0] if len(numpy.unique(inputNodePositions[:, i])) > 1 else numpy.unique(inputNodePositions[:, i])[0] for i in range(3)])
+        args.NEIGHBOUR_RADIUS = 4*int(numpy.mean(nodeSpacing))
+        print(f"Neither number of neighbours nor neighbour distance set, using default distance of 4*mean(nodeSpacing) = {args.NEIGHBOUR_RADIUS}")
+    else:
+        # Discrete input field
+        args.NUMBER_OF_NEIGHBOURS = 27
+        print("Neither number of neighbours nor neighbour distance set, using default 27 neighbours")
+
+if args.LQC:
+    print("\n\nCalculate coherency")
+    LQC = spam.DIC.estimateLocalQuadraticCoherency(inputNodePositions, inputDisplacements, neighbourRadius=args.NEIGHBOUR_RADIUS, nNeighbours=args.NUMBER_OF_NEIGHBOURS, verbose=True)
+
+if args.LQCF:
+    print("\n\nCorrection based on local quadratic coherency")
+    dispLQC = spam.DIC.estimateDisplacementFromQuadraticFit(inputNodePositions, inputDisplacements, LQC, neighbourRadius=args.NEIGHBOUR_RADIUS, nNeighbours=args.NUMBER_OF_NEIGHBOURS, verbose=True)
+
+    # pass the displacements
+    # TODO: what about the other components of F?
+    outputPhiField[:, 0:3, -1] = dispLQC
+
+#Outputs are:
+  #- TSV files
+  #- (optional) VTK files for visualisation
+  #- (optional) TIF files in the case of gridded data
+if args.TSV:
+    if discrete:
+        TSVheader = "Label\tZpos\tYpos\tXpos\tFzz\tFzy\tFzx\tZdisp\tFyz\tFyy\tFyx\tYdisp\tFxz\tFxy\tFxx\tXdisp\tpixelSearchCC\treturnStatus\terror\tdeltaPhiNorm\titerations"
+    else:
+        TSVheader = "NodeNumber\tZpos\tYpos\tXpos\tFzz\tFzy\tFzx\tZdisp\tFyz\tFyy\tFyx\tYdisp\tFxz\tFxy\tFxx\tXdisp\tpixelSearchCC\treturnStatus\terror\tdeltaPhiNorm\titerations"
+    outMatrix = numpy.array([numpy.array(inputNodePositions.shape[0]),
+                            inputNodePositions[:, 0],  inputNodePositions[:, 1],  inputNodePositions[:, 2],
+                            outputPhiField[:, 0, 0],    outputPhiField[:, 0, 1],    outputPhiField[:, 0, 2],  outputPhiField[:, 0, 3],
+                            outputPhiField[:, 1, 0],    outputPhiField[:, 1, 1],    outputPhiField[:, 1, 2],  outputPhiField[:, 1, 3],
+                            outputPhiField[:, 2, 0],    outputPhiField[:, 2, 1],    outputPhiField[:, 2, 2],  outputPhiField[:, 2, 3],
+                            outputPixelSearchCC,
+                            outputReturnStatus,
+                            outputError,
+                            outputDeltaPhiNorm,
+                            outputIterations]).T
+
+    numpy.savetxt(args.OUT_DIR+"/"+args.PREFIX+".tsv",
+                outMatrix,
+                fmt='%.7f',
+                delimiter='\t',
+                newline='\n',
+                comments='',
+                header=TSVheader)
+
+if args.TIFF:
+    if grid:
+        if inputNodesDim[0] != 1:
+            tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-Zdisp.tif", outputPhiField[:, 0, -1].astype('<f4').reshape(inputNodesDim))
+        tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-Ydisp.tif",     outputPhiField[:, 1, -1].astype('<f4').reshape(inputNodesDim))
+        tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-Xdisp.tif",     outputPhiField[:, 2, -1].astype('<f4').reshape(inputNodesDim))
+        #tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-CC.tif",                     pixelSearchCC.astype('<f4').reshape(nodesDim))
+        #tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-returnStatus.tif",           returnStatus.astype('<f4').reshape(nodesDim))
+    else:
+        # Think about relabelling grains here automatically?
+        pass
+
+
+# Collect data into VTK output
+if args.VTK:
+    if grid:
+        cellData = {}
+        cellData['displacements'] = outputPhiField[:, :-1, 3].reshape((inputNodesDim[0], inputNodesDim[1], inputNodesDim[2], 3))
+
+        # Overwrite nans and infs with 0, rubbish I know
+        cellData['displacements'][numpy.logical_not(numpy.isfinite(cellData['displacements']))] = 0
+
+        # This is perfect in the case where NS = 2xHWS, these cells will all be in the right place
+        #   In the case of overlapping of under use of data, it should be approximately correct
+        # If you insist on overlapping, then perhaps it's better to save each point as a cube glyph
+        #   and actually *have* overlapping
+        # HACK assume HWS is half node spacing
+        nodeSpacing = numpy.array([numpy.unique(inputNodePositions[:, i])[1] - numpy.unique(inputNodePositions[:, i])[0] if len(numpy.unique(inputNodePositions[:, i])) > 1 else numpy.unique(inputNodePositions[:, i])[0] for i in range(3)])
+        HWS = nodeSpacing/2
+        spam.helpers.writeStructuredVTK(origin=inputNodePositions[0]-HWS, aspectRatio=nodeSpacing, cellData=cellData, fileName=args.OUT_DIR+"/"+args.PREFIX+".vtk")
+
+    else:
+        disp = outputPhiField[:, 0:3, -1]
+        disp[numpy.logical_not(numpy.isfinite(disp))] = 0
+
+        magDisp = numpy.linalg.norm(disp, axis=1)
+
+        VTKglyphDict = {'displacements': outputPhiField[:, 0:3, -1],
+                        'mag(displacements)': magDisp
+                        }
+
+        spam.helpers.writeGlyphsVTK(inputNodePositions, VTKglyphDict, fileName=args.OUT_DIR + "/" + args.PREFIX + ".vtk")

setup.py

@@ -157,6 +157,7 @@ scripts = [
            'scripts/spam-mmr-graphical',
            'scripts/spam-moveLabels',
            'scripts/spam-passPhiField',
+           'scripts/spam-filterPhiField',
            'scripts/spam-pixelSearch',
            'scripts/spam-pixelSearchPropagate',
            'scripts/spam-reg',

tools/helpers/optionsParser.py

@@ -2674,3 +2674,158 @@ def passPhiField(parser):
         exit()
 
     return args
+
+def filterPhiField(parser):
+    parser.add_argument('-pf',
+                        '-phiFile',
+                        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")
+
+    parser.add_argument('-pfb',
+                        '--phiFile-bin-ratio',
+                        type=int,
+                        default=1,
+                        dest='PHIFILE_BIN_RATIO',
+                        help="Ratio of binning level between 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('-np',
+                        '--number-of-processes',
+                        default=None,
+                        type=int,
+                        dest='PROCESSES',
+                        help="Number of parallel processes to use. Default = multiprocessing.cpu_count()")
+
+    parser.add_argument('-nr',
+                        '--neighbourhood-radius-px',
+                        type=float,
+                        default=None,
+                        dest='NEIGHBOUR_RADIUS',
+                        help="Radius (in pixels) inside which to select neighbours for field interpolation. Excludes -nn option")
+
+    parser.add_argument('-nn',
+                        '--number-of-neighbours',
+                        type=int,
+                        default=None,
+                        dest='NUMBER_OF_NEIGHBOURS',
+                        help="Number of neighbours for field interpolation. Default = None (radius mode is default)")
+
+    parser.add_argument('-lqc',
+                        '--local-quadratic-coherency',
+                        action="store_true",
+                        dest='LQC',
+                        help='Find incoherent points based on local quadratic coherency?')
+
+    parser.add_argument('-lqcf',
+                        '--local-quadratic-coherency-fit',
+                        action="store_true",
+                        dest='LQCF',
+                        help='Correct bad points based on a local quadratic coherency fit of good neighbours?')
+
+    parser.add_argument('-i',
+                        '--interpolate-field',
+                        action="store_true",
+                        dest='CORRECT_FIELD',
+                        help='Correct bad points based on an inverse distance interpolation of good neighbours?')
+
+    parser.add_argument('-rst',
+                        '--return-status-threshold',
+                        type=int,
+                        default=-4,
+                        dest='RETURN_STATUS_THRESHOLD',
+                        help='Lowest return status value to preserve in input PhiField. Default = -4')
+
+    parser.add_argument('-dpt',
+                        '--delta-phi-norm-threshold',
+                        type=numpy.float,
+                        default=0.001,
+                        dest='DELTA_PHI_NORM_THRESHOLD',
+                        help="Delta Phi norm threshold BELOW which to consider the point good. Only for a point with return status = 1 . Default = 0.001")
+
+    parser.add_argument('-pscct',
+                        '--pixel-search-cc-threshold',
+                        type=numpy.float,
+                        default=0,
+                        dest='PIXEL_SEARCH_CC_THRESHOLD',
+                        help="Pixel search correlation coefficient threshold ABOVE which to consider the point good. Default = 0.9")
+
+    parser.add_argument('-mf',
+                        '--median-filter',
+                        action="store_true",
+                        dest='CORRECT_MEDIAN_FILTER',
+                        help="Activates an overall median filter on the input displacement field")
+
+    parser.add_argument('-mfr',
+                        '--median-filter-radius',
+                        type=int,
+                        default=2,
+                        dest='MEDIAN_FILTER_RADIUS',
+                        help="Radius (in pixels) of median filter. Default = 2")
+
+    parser.add_argument('-od',
+                        '--out-dir',
+                        type=str,
+                        default=None,
+                        dest='OUT_DIR',
+                        help='Output directory, default is the dirname of input file')
+
+    parser.add_argument('-pre',
+                        '--prefix',
+                        type=str,
+                        default=None,
+                        dest='PREFIX',
+                        help='Prefix for output files (without extension). Default is basename of input file')
+
+    parser.add_argument('-tif',
+                        '-tiff',
+                        action="store_true",
+                        dest='TIFF',
+                        help='Activate TIFF output format. Default = False')
+
+    parser.add_argument('-notsv',
+                        '-noTSV',
+                        action="store_false",
+                        dest='TSV',
+                        help='Disactivate TSV output format?')
+
+    parser.add_argument('-vtk',
+                        '--VTKout',
+                        action="store_true",
+                        dest='VTK',
+                        help='Activate VTK output format. Default = False')
+
+    args = parser.parse_args()
+
+    if args.PHIFILE is None:
+        print("This function definitely needs a TSV Phi file input")
+        exit()
+
+
+    # If we have no out dir specified, deliver on our default promise -- this can't be done inline before since parser.parse_args() has not been run at that stage.
+    if args.OUT_DIR is None:
+        args.OUT_DIR = os.path.dirname(args.PHIFILE.name)
+        # However if we have no dir, notice this and make it the current directory.
+        if args.OUT_DIR == "":
+            args.OUT_DIR = "./"
+    else:
+        # Check existence of output directory
+        try:
+            if args.OUT_DIR:
+                os.makedirs(args.OUT_DIR)
+            else:
+                args.DIR_out = os.path.dirname(args.PHIFILE.name)
+        except OSError:
+            if not os.path.isdir(args.OUT_DIR):
+                raise
+
+    # Output file name prefix
+    if args.PREFIX is None:
+        args.PREFIX = os.path.splitext(os.path.basename(args.PHIFILE.name))[0] + "-filtered"
+    else:
+        args.PREFIX += "-filtered"
+
+    if args.LQCF:
+        args.PREFIX += "-LQC"
+
+    return args
\ No newline at end of file