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

[skip-ci] interpolatePhiField now general, spam-filterPhiField starting to work

parent e9d53685
Pipeline #66157 skipped
......@@ -67,6 +67,8 @@ parser = argparse.ArgumentParser(description="spam-filterPhiField "+spam.helpers
# Parse arguments with external helper function
args = spam.helpers.optionsParser.filterPhiField(parser)
if args.PROCESSES is None: args.PROCESSES = multiprocessing.cpu_count()
print("spam-filterPhiField -- Current Settings:")
argsDict = vars(args)
for key in sorted(argsDict):
......@@ -105,6 +107,10 @@ inputPixelSearchCC = PhiFromFile["pixelSearchCC"]
inputDeltaPhiNorm = PhiFromFile["deltaPhiNorm"]
inputIterations = PhiFromFile["iterations"]
inputError = PhiFromFile["error"]
### Empty arrays for
inputGood = numpy.zeros(inputNodePositions.shape[0], dtype=bool)
inputBad = numpy.zeros(inputNodePositions.shape[0], dtype=bool)
inputIgnore = numpy.zeros(inputNodePositions.shape[0], dtype=bool)
# output arrays
outputPhiField = numpy.zeros((inputNodePositions.shape[0], 4, 4))
......@@ -113,8 +119,6 @@ outputDeltaPhiNorm = numpy.ones( (inputNodePositions.shape[0]), dtype=float)*10
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")
......@@ -131,17 +135,95 @@ if args.NEIGHBOUR_RADIUS is None and args.NUMBER_OF_NEIGHBOURS is None:
args.NUMBER_OF_NEIGHBOURS = 27
print("Neither number of neighbours nor neighbour distance set, using default 27 neighbours")
if args.LQC:
###############################################################
### Define IGNORE points:
###############################################################
if args.MASK:
inputIgnore = inputReturnStatus < -4
###############################################################
### Apply threshold to select good and bad points
###############################################################
if args.SRS:
print(f"\n\nSelecting bad points as Return Status <= {args.SRST}")
inputGood = numpy.logical_and(inputReturnStatus > args.SRST, ~inputIgnore)
inputBad = numpy.logical_and(inputReturnStatus <= args.SRST, ~inputIgnore)
if args.SLQC:
print("\tYou passed -slqc but you can only have one selection at a time")
if args.SCC:
print("\tYou passed -scc but you can only have one selection at a time")
elif args.SCC:
print(f"\n\nSelecting bad points with Pixel Search CC <= {args.SCCT}")
inputGood = numpy.logical_and(inputPixelSearchCC > args.SCCT, ~inputIgnore)
inputBad = numpy.logical_and(inputPixelSearchCC <= args.SCCT, ~inputIgnore)
if args.SLQC:
print("\tYou passed -slqc but you can only have one selection at a time")
elif args.SLQC:
print("\n\nCalculate coherency")
LQC = spam.DIC.estimateLocalQuadraticCoherency(inputNodePositions, inputDisplacements, neighbourRadius=args.NEIGHBOUR_RADIUS, nNeighbours=args.NUMBER_OF_NEIGHBOURS, verbose=True)
LQC = spam.DIC.estimateLocalQuadraticCoherency(inputNodePositions[~inputIgnore],
inputDisplacements[~inputIgnore],
neighbourRadius=args.NEIGHBOUR_RADIUS,
nNeighbours=args.NUMBER_OF_NEIGHBOURS,
nProcesses=args.PROCESSES,
verbose=True)
print(LQC.shape)
print(inputGood[~inputIgnore].shape)
inputGood[~inputIgnore] = LQC < 0.1
inputBad[~inputIgnore] = LQC >= 0.1
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)
###############################################################
### Copy over the values for good AND ignore to output
###############################################################
gandi = numpy.logical_or(inputGood, inputIgnore)
outputPhiField[gandi] = inputPhiField[gandi]
outputReturnStatus[gandi] = inputReturnStatus[gandi]
outputDeltaPhiNorm[gandi] = inputDeltaPhiNorm[gandi]
outputIterations[gandi] = inputIterations[gandi]
outputError[gandi] = inputError[gandi]
outputPixelSearchCC[gandi] = inputPixelSearchCC[gandi]
if (args.CINT + args.CLQF) > 0 and numpy.sum(inputBad) == 0:
print("No points to correct, exiting")
exit()
###############################################################
### Correct those bad points
###############################################################
if args.CINT:
print(f"\n\nCorrection based on local interpolation (mode = {args.MODE})")
PhiFieldCorrected = spam.DIC.interpolatePhiField(inputNodePositions[inputGood],
inputPhiField[inputGood],
inputNodePositions[inputBad],
nNeighbours=args.NUMBER_OF_NEIGHBOURS,
neighbourRadius=args.NEIGHBOUR_RADIUS,
mode=args.MODE,
nProcesses=args.PROCESSES,
verbose=True)
outputPhiField[inputBad] = PhiFieldCorrected
outputReturnStatus[inputBad] = 1
if args.CLQF:
print("\tYou asked to correct with local QC fitting with -clqf, but only one correciton mode is supported")
elif args.CLQF:
print("\n\nCorrection based on local quadratic coherency")
dispLQC = spam.DIC.estimateDisplacementFromQuadraticFit(inputNodePositions[inputGood],
inputDisplacements[inputGood],
inputNodePositions[inputBad],
neighbourRadius=args.NEIGHBOUR_RADIUS,
nNeighbours=args.NUMBER_OF_NEIGHBOURS,
nProcesses=args.PROCESSES,
verbose=True)
# pass the displacements
# TODO: what about the other components of F?
outputPhiField[:, 0:3, -1] = dispLQC
outputPhiField[inputBad, 0:3, -1] = dispLQC
outputReturnStatus[inputBad] = 1
#Outputs are:
#- TSV files
......@@ -152,7 +234,7 @@ if args.TSV:
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]),
outMatrix = numpy.array([numpy.arange(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],
......
......@@ -261,105 +261,14 @@ else:
#else:
# TODO: Last case with DDIC in and DDIC out could be with NNEIGHBOURS
# Create the k-d tree of the coordinates of the input F field
inputTree = scipy.spatial.KDTree(goodInputNodePositions)
def interpolateOnePoint(outputNode):
outputNodePosition = outputNodePositions[outputNode]
outputPhi = numpy.zeros((4, 4), dtype='float')
if args.NEIGHBOUR_RADIUS is not None:
# Ball lookup
indices = inputTree.query_ball_point(outputNodePosition, args.NEIGHBOUR_RADIUS)
if len(indices) == 0:
#print("no point!")
return outputNode, numpy.eye(4) * numpy.nan
else:
distances = numpy.linalg.norm(outputNodePosition - goodInputNodePositions[indices], axis=1)
else:
# Number of Neighbour lookup
distances, indices = inputTree.query(outputNodePosition, k=args.NUMBER_OF_NEIGHBOURS)
indices = numpy.array(indices)
distances = numpy.array(distances)
#print(numpy.max(indices))
# Check if there is only one neighbour
if indices.size == 1:
if args.MODE in ["all", "rigid"]: # We need to interpolate all 12 components of Phi
outputPhi = goodInputPhiField[indices].copy()
if args.MODE == "rigid":
outputPhi = spam.deformation.computeRigidPhi(outputPhi)
else: # interpolate only displacements
outputPhi[0:3, 0:3] = numpy.eye(3)
outputPhi[0:3, -1] = goodInputPhiField[indices, 0:3, -1].copy()
return outputNode, outputPhi
# If > 1 neighbour, interpolate based on distance
else:
weights = (1/(distances+tol))
weightsTotal = float(weights.sum())
if args.MODE in ["all", "rigid"]: # We need to interpolate all 12 components of Phi
for neighbourNumber, neighbourIndex in enumerate(indices):
weightInv = weights[neighbourNumber] / weightsTotal
outputPhi[:-1] += goodInputPhiField[neighbourIndex][:-1]*weightInv
if args.MODE == "rigid":
outputPhi = spam.deformation.computeRigidPhi(outputPhi)
else: # interpolate only displacements
outputPhi[0:3, 0:3] = numpy.eye(3)
for neighbourNumber, neighbourIndex in enumerate(indices):
weightInv = weights[neighbourNumber] / weightsTotal
outputPhi[0:3,-1] += goodInputPhiField[neighbourIndex][0:3,-1]*weightInv
return outputNode, outputPhi
if args.PROCESSES is None: args.PROCESSES = multiprocessing.cpu_count()
print("\n\tStarting Phi field interpolation (with {} process{})".format(args.PROCESSES, 'es' if args.PROCESSES > 1 else ''))
widgets = [progressbar.Bar(), ' ', progressbar.AdaptiveETA()]
pbar = progressbar.ProgressBar(widgets=widgets, maxval=outputNumberOfNodes)
pbar.start()
finishedNodes = 0
with multiprocessing.Pool(processes=args.PROCESSES) as pool:
for returns in pool.imap_unordered(interpolateOnePoint, range(outputNumberOfNodes)):
finishedNodes += 1
# Update progres bar if point is not skipped
pbar.update(finishedNodes)
outputPhiField[returns[0]] = returns[1]
pool.close()
pool.join()
pbar.finish()
print("\n")
## Calculate the distance of the point of the current grid to the points of the input F field
## Check if we've hit the same point in the two grids
#if numpy.any(distance == 0):
## Copy F of that point to the F in the current grid point
#PhiField[node] = fieldValues[indices][numpy.where(distance == 0)].copy()
## If not, consider the closest neighbours
#else:
## Compute the "Inverse Distance Weighting" since the closest points should have the major influence
#weightSumInv = sum(1/distance)
## Loop over each neighbour
#for neighbour in range(nNeighbours):
## Calculate it's weight
#weightInv = (1/distance[neighbour])/float(weightSumInv)
## Fill the F of the current grid point with the weighted F components of the ith nearest neighbour in the input F field
#PhiField[node][:-1] += fieldValues[indices[neighbour]][:-1]*weightInv
outputPhiField = spam.DIC.interpolatePhiField(goodInputNodePositions,
goodInputPhiField,
outputNodePositions,
nNeighbours=args.NUMBER_OF_NEIGHBOURS,
neighbourRadius=args.NEIGHBOUR_RADIUS,
mode=args.MODE,
nProcesses=args.PROCESSES,
verbose=True)
#- apply a registration (single Phi) to a series of points:
#- defined on a grid with NS
......
......@@ -517,7 +517,6 @@ class testAll(unittest.TestCase):
exitCode = subprocess.call(["spam-regularStrain",
testFolder + "snow-grid-5a-ldic.tsv",
"-comp", "r",
"-mask",
"-r", "1",
"-od", testFolder, "-pre", "snow-grid-6a",
"-tif"],
......@@ -539,7 +538,7 @@ class testAll(unittest.TestCase):
exitCode = subprocess.call(["spam-regularStrain", "--Q8",
testFolder + "snow-grid-5a-ldic.tsv",
"-comp", "U", "e", "vol", "volss", "dev", "devss",
"-vtk", "-mask",
"-vtk",
"-od", testFolder, "-pre", "snow-grid-6b"],
stdout=FNULL, stderr=FNULL)
self.assertEqual(exitCode, 0)
......
This diff is collapsed.
......@@ -1024,11 +1024,11 @@ def regularStrainParser(parser):
dest='PROCESSES',
help="Number of parallel processes to use. Default = multiprocessing.cpu_count()")
parser.add_argument('-mask',
'--mask',
action="store_true",
parser.add_argument('-nomask',
'--nomask',
action="store_false",
dest='MASK',
help='Mask correlation points according to return status')
help="Don't mask correlation points in background according to return status (i.e., include RS=-5 or less)")
parser.add_argument('-r',
'--neighbourhood-radius-for-strain-calculation',
......@@ -2739,6 +2739,12 @@ def filterPhiField(parser):
dest='PROCESSES',
help="Number of parallel processes to use. Default = multiprocessing.cpu_count()")
parser.add_argument('-nomask',
'--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)")
parser.add_argument('-nr',
'--neighbourhood-radius-px',
type=float,
......@@ -2753,44 +2759,71 @@ def filterPhiField(parser):
dest='NUMBER_OF_NEIGHBOURS',
help="Number of neighbours for field interpolation. Default = None (radius mode is default)")
parser.add_argument('-lqc',
'--local-quadratic-coherency',
parser.add_argument('-srs',
'--select-return-status',
action="store_true",
dest='LQC',
help='Find incoherent points based on local quadratic coherency?')
dest='SRS',
help='Select bad points for correction based on Return Status? This will use -srst as a threshold')
parser.add_argument('-lqcf',
'--local-quadratic-coherency-fit',
parser.add_argument('-srst',
'--select-return-status-threshold',
type=int,
default=1,
dest='SRST',
help='Return Status Threshold for selecting bad points. Default = 1 or less')
parser.add_argument('-scc',
'--select-cc',
action="store_true",
dest='LQCF',
help='Correct bad points based on a local quadratic coherency fit of good neighbours?')
dest='SCC',
help='Select bad points for correction based on Pixel Search CC? This will use -scct as a threshold')
parser.add_argument('-i',
'--interpolate-field',
parser.add_argument('-scct',
'--select-cc-threshold',
type=float,
default=0.99,
dest='SCCT',
help='Pixel Search CC for selecting bad points. Default = 0.99 or less')
parser.add_argument('-slqc',
'--select-local-quadratic-coherency',
action="store_true",
dest='CORRECT_FIELD',
help='Correct bad points based on an inverse distance interpolation of good neighbours?')
dest='SLQC',
help='Select bad points for correction based on local quadratic coherency? Threshold = 0.1 or more')
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('-cint',
'--correct-by-interpolation',
action="store_true",
dest='CINT',
help='Correct with a local interpolation with weights equal to the inverse of the distance?')
parser.add_argument('-cintm',
'--correct-by-interpolation-mode',
type=str,
default='all',
dest='MODE',
help="What do you want to interpolate? Options: 'all': the full Phi, 'rigid': Rigid body motion, 'disp': Only displacements (faster). Default = 'all'")
parser.add_argument('-clqf',
'--correct-by-local-quadratic-fit',
action="store_true",
dest='CLQF',
help='Correct with a local interpolation with weights equal to the inverse of the distance?')
parser.add_argument('-lqcf',
'--local-quadratic-coherency-fit',
action="store_true",
dest='LQCF',
help='Correct bad points based on a fitting of the local quadratic coherency?')
#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('-mf',
'--median-filter',
......
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