Skip to content

GitLab

Commit 9771b4e7 authored by Edward Andò's avatar Edward Andò
Browse files

spam-filterPhiField finished

parent ae8acaf3
Pipeline #66174 passed with stages
in 13 minutes and 47 seconds
Hide whitespace changes
Inline Side-by-side
scripts/spam-filterPhiField
View file @ 9771b4e7
......@@ -52,6 +52,7 @@ import spam.label
import numpy
import multiprocessing
import scipy.spatial
import scipy.ndimage
import progressbar
import argparse
import tifffile
......@@ -107,7 +108,7 @@ inputPixelSearchCC = PhiFromFile["pixelSearchCC"]
inputDeltaPhiNorm = PhiFromFile["deltaPhiNorm"]
inputIterations = PhiFromFile["iterations"]
inputError = PhiFromFile["error"]
### Empty arrays for
### Empty arrays for masking points
inputGood = numpy.zeros(inputNodePositions.shape[0], dtype=bool)
inputBad = numpy.zeros(inputNodePositions.shape[0], dtype=bool)
inputIgnore = numpy.zeros(inputNodePositions.shape[0], dtype=bool)
......@@ -225,6 +226,60 @@ elif args.CLQF:
if args.FILTER_MEDIAN:
if discrete:
print("Median filter for discrete mode not implemented... does it even make sense?")
else:
# Filter ALL POINTS
# if asked, apply a median filter of a specific size in the Phi field
print("\nApplying median filter...")
filterPointsRadius = int(args.FILTER_MEDIAN_RADIUS)
if args.MASK:
inputPhiField[inputIgnore] = numpy.nan
if args.MODE == 'rigid':
print("Rigid mode not well defined for overall median filtering, exiting")
exit()
if args.MODE == 'all':
# Filter F components
print("Filtering F components...")
print("\t1/9")
outputPhiField[:, 0, 0] = scipy.ndimage.generic_filter(inputPhiField[:, 0, 0].reshape(inputNodesDim), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
print("\t2/9")
outputPhiField[:, 1, 0] = scipy.ndimage.generic_filter(inputPhiField[:, 1, 0].reshape(inputNodesDim), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
print("\t3/9")
outputPhiField[:, 2, 0] = scipy.ndimage.generic_filter(inputPhiField[:, 2, 0].reshape(inputNodesDim), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
print("\t4/9")
outputPhiField[:, 0, 1] = scipy.ndimage.generic_filter(inputPhiField[:, 0, 1].reshape(inputNodesDim), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
print("\t5/9")
outputPhiField[:, 1, 1] = scipy.ndimage.generic_filter(inputPhiField[:, 1, 1].reshape(inputNodesDim), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
print("\t6/9")
outputPhiField[:, 2, 1] = scipy.ndimage.generic_filter(inputPhiField[:, 2, 1].reshape(inputNodesDim), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
print("\t7/9")
outputPhiField[:, 0, 2] = scipy.ndimage.generic_filter(inputPhiField[:, 0, 2].reshape(inputNodesDim), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
print("\t8/9")
outputPhiField[:, 1, 2] = scipy.ndimage.generic_filter(inputPhiField[:, 1, 2].reshape(inputNodesDim), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
print("\t9/9")
outputPhiField[:, 2, 2] = scipy.ndimage.generic_filter(inputPhiField[:, 2, 2].reshape(inputNodesDim), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
if args.MODE == 'disp':
for n in range(inputNodePositions.shape[0]):
outputPhiField[n] = numpy.eye(4)
print("Filtering displacements...")
print("\t1/3")
outputPhiField[:, 0, -1] = scipy.ndimage.generic_filter(inputPhiField[:, 0, -1].reshape(inputNodesDim), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
print("\t2/3")
outputPhiField[:, 1, -1] = scipy.ndimage.generic_filter(inputPhiField[:, 1, -1].reshape(inputNodesDim), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
print("\t3/3")
outputPhiField[:, 2, -1] = scipy.ndimage.generic_filter(inputPhiField[:, 2, -1].reshape(inputNodesDim), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
if args.MASK:
outputPhiField[inputIgnore] = numpy.nan
#Outputs are:
#- TSV files
#- (optional) VTK files for visualisation
......
tools/DIC/kinematics.py
View file @ 9771b4e7
......@@ -16,6 +16,7 @@ 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 spam.deformation
import numpy
import tifffile
import scipy.spatial
......@@ -239,6 +240,8 @@ def estimateDisplacementFromQuadraticFit(fieldCoords, displacements, pointsToEst
# check if neighbours are selected based on radius or based on number, default is radius
if (nNeighbours is None) == (neighbourRadius is None):
print("spam.DIC.estimateDisplacementFromQuadraticFit(): One and only one of nNeighbours and neighbourRadius must be passed")
ball = None
if nNeighbours is not None:
ball = False
elif neighbourRadius is not None:
......@@ -419,12 +422,13 @@ def interpolatePhiField(fieldCoords, PhiField, pointsToInterpolate, nNeighbours=
def interpolateOnePoint(pointNumber):
pointToInterpolate = pointsToInterpolate[pointNumber]
outputPhi = numpy.zeros((4, 4), dtype=PhiField.dtype)
outputPhi[-1] = [0, 0, 0, 1]
if ball:
# Ball lookup
indices = treeCoord.query_ball_point(pointToInterpolate, neighbourRadius)
if len(indices) == 0:
#print("no point!")
# No point!
return pointNumber, numpy.eye(4) * numpy.nan
else:
distances = numpy.linalg.norm(pointToInterpolate - fieldCoords[indices], axis=1)
......
tools/helpers/optionsParser.py
View file @ 9771b4e7
......@@ -2792,19 +2792,19 @@ def filterPhiField(parser):
help='Select bad points for correction based on local quadratic coherency? Threshold = 0.1 or more')
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?')
help='Correct with a local interpolation with weights equal to the inverse of the distance? -mode applies')
parser.add_argument('-cintm',
'--correct-by-interpolation-mode',
parser.add_argument('-m',
'-mode',
'--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'")
help="What do you want to interpolate/filter? Options: 'all': the full Phi, 'rigid': Rigid body motion, 'disp': Only displacements (faster). Default = 'all'.")
parser.add_argument('-clqf',
'--correct-by-local-quadratic-fit',
......@@ -2825,18 +2825,18 @@ def filterPhiField(parser):
#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',
parser.add_argument('-fm',
'--filter-median',
action="store_true",
dest='CORRECT_MEDIAN_FILTER',
help="Activates an overall median filter on the input displacement field")
dest='FILTER_MEDIAN',
help="Activates an overall median filter on the input Phi Field. -mode 'all' or 'disp' can be applied")
parser.add_argument('-mfr',
'--median-filter-radius',
parser.add_argument('-fmr',
'--filter-median-radius',
type=int,
default=2,
dest='MEDIAN_FILTER_RADIUS',
help="Radius (in pixels) of median filter. Default = 2")
default=1,
dest='FILTER_MEDIAN_RADIUS',
help="Radius (in pixels) of median filter. Default = 1")
parser.add_argument('-od',
'--out-dir',
......@@ -2894,6 +2894,10 @@ def filterPhiField(parser):
if not os.path.isdir(args.OUT_DIR):
raise
if (args.SRS + args.SCC + args.SLQC + args.CINT + args.CLQF) > 1 and args.FILTER_MEDIAN:
print("WARINING: you can't ask for an overall median filter and a correction")
exit()
# Output file name prefix
if args.PREFIX is None:
args.PREFIX = os.path.splitext(os.path.basename(args.PHIFILE.name))[0] + "-filtered"
......
tools/tests/test_DIC.py
View file @ 9771b4e7
......@@ -104,62 +104,6 @@ class TestFunctionDVC(unittest.TestCase):
self.assertEqual(twoDtest is None, True)
def test_correctPhiField(self):
# case 0: wrong tsv name
f0 = spam.DIC.correctPhiField(fileName="./sapmPhiFieldCF.tsv")
# case 0b: don't pass all separate arrays
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.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.DIC.correctPhiField(fileName=None,
fieldCoords=fieldCoords, fieldValues=PhiField.copy(),
fieldRS=RS.copy(), fieldDPhi=DPhi.copy(),
fieldPixelSearchCC=numpy.zeros_like(RS), fieldIT=IT,
ignoreBadPoints=True, correctBadPoints=False, saveFile=False)
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.DIC.correctPhiField(fileName=None,
fieldCoords=fieldCoords, fieldValues=PhiField,
fieldRS=RS, fieldDPhi=numpy.full((fieldCoords.shape[0]), 100),
fieldPixelSearchCC=numpy.ones_like(RS), fieldIT=IT,
ignoreBadPoints=True, correctBadPoints=False, saveFile=False)
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.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.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.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.DIC.correctPhiField(fieldCoords=fieldCoords, fieldValues=PhiField.copy(),
fieldRS=RS, fieldDPhi=numpy.full((fieldCoords.shape[0]), 100),
fieldPixelSearchCC=PSCC, fieldIT=IT,
correctBadPoints=True, ignoreBackGround=True,
verbose=True, saveFile=True)
self.assertEqual(numpy.nansum(f3b[:, 0:3, -1]), (fieldCoords.shape[0] - 2) * 3)
os.remove("spam-corrected-N12.tsv")
os.remove("spamPhiFieldCF.tsv")
def test_applyPhiField(self):
a = numpy.random.rand(10, 10, 10)
Phi = spam.deformation.computePhi({'t': [0.0, 3.0, 3.0]})
......@@ -658,10 +602,21 @@ class TestFunctionDVC(unittest.TestCase):
self.assertTrue(len(numpy.where(LQC3 > 0.1)[0]), 13)
# correct disp with radius
dispModFit2 = spam.DIC.estimateDisplacementFromQuadraticFit(points, dispMod, LQC2, neighbourRadius=20)
dispModFit2corrected = spam.DIC.estimateDisplacementFromQuadraticFit(points[ LQC2 < 0.1],
dispMod[LQC2 < 0.1],
points[ LQC2 >= 0.1],
neighbourRadius=20)
dispModFit2 = dispMod.copy()
dispModFit2[LQC2 >= 0.1] = dispModFit2corrected
self.assertTrue(numpy.allclose(dispModFit2, disp))
# correct disp with neighbours
dispModFit3 = spam.DIC.estimateDisplacementFromQuadraticFit(points, dispMod, LQC3, nNeighbours=35)
dispModFit3corrected = spam.DIC.estimateDisplacementFromQuadraticFit(points[ LQC3 < 0.1],
dispMod[LQC3 < 0.1],
points[ LQC3 >= 0.1],
nNeighbours=35)
dispModFit3 = dispMod.copy()
dispModFit3[LQC3 >= 0.1] = dispModFit3corrected
self.assertTrue(numpy.allclose(dispModFit3, disp))
......
tools/tests/test_deformation.py
View file @ 9771b4e7
......@@ -739,31 +739,61 @@ class testAll(unittest.TestCase):
# z = 0 --> transformation below:
transformation = {'t': [5., 0., 0.],
'r': [0., 0., 0.]}
Phi = spam.deformation.computePhi(transformation)
#print(Phi)
centre = [50, 50, 50]
# z = 100 --> nothing
fieldCoords = numpy.array([[ 0, 0, 0],
[ 0, 0,100],
[ 0,100, 0],
[ 0,100,100],
[100, 0, 0],
[100, 0, 0],
[100, 0,100],
[100,100, 0],
[100,100,100]])
fieldValues = numpy.zeros([8,4,4])
fieldValues[0:4] = spam.deformation.computePhi(transformation)
fieldValues[4:8] = numpy.eye(4)
fieldValues = numpy.zeros([fieldCoords.shape[0],4,4])
#fieldValues[0:4] = spam.deformation.computePhi(transformation)
#fieldValues[4:8] = numpy.eye(4)
for n in range(fieldCoords.shape[0]):
fieldValues[n] = numpy.eye(4)
fieldValues[n, 0:3, 0:3] = Phi[0:3, 0:3]
fieldValues[n, 0:3, -1] = spam.deformation.decomposePhi(Phi, PhiPoint=fieldCoords[n], PhiCentre=centre)['t']
interpCoords = numpy.array([[50.,50.,50.]])
interpolatedPhi = spam.DIC.interpolatePhiField(fieldCoords, fieldValues, interpCoords, nNeighbours=1, verbose=True)
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, -1], transformation['t'], atol=0.01))
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, 0:3], Phi[0:3, 0:3], atol=0.01))
interpolatedPhi = spam.DIC.interpolatePhiField(fieldCoords, fieldValues, interpCoords, nNeighbours=1, mode='rigid')
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, -1], transformation['t'], atol=0.01))
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, 0:3], Phi[0:3, 0:3], atol=0.01))
interpolatedPhi = spam.DIC.interpolatePhiField(fieldCoords, fieldValues, interpCoords, nNeighbours=1, mode='disp')
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, -1], transformation['t'], atol=0.01))
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, 0:3], numpy.eye(3), atol=0.01))
interpolatedPhi = spam.DIC.interpolatePhiField(fieldCoords, fieldValues, interpCoords, nNeighbours=2)
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, -1], transformation['t'], atol=0.01))
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, 0:3], Phi[0:3, 0:3], atol=0.01))
interpCoods = numpy.array([[50.,50.,50.]])
interpolatedPhi = spam.DIC.interpolatePhiField(fieldCoords, fieldValues, interpCoords, nNeighbours=2, mode='rigid')
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, -1], transformation['t'], atol=0.01))
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, 0:3], Phi[0:3, 0:3], atol=0.01))
interpolatedPhi = spam.DIC.interpolatePhiField(fieldCoords, fieldValues, interpCoods)
interpolatedPhi = spam.DIC.interpolatePhiField(fieldCoords, fieldValues, interpCoords, nNeighbours=2, mode='disp')
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, -1], transformation['t'], atol=0.01))
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, 0:3], numpy.eye(3), atol=0.01))
decomposedInterpolatedPhi = spam.deformation.decomposePhi(interpolatedPhi[0])
interpolatedPhi = spam.DIC.interpolatePhiField(fieldCoords, fieldValues, interpCoords, neighbourRadius=10)
self.assertTrue(numpy.isnan(interpolatedPhi[0]).sum() == 16)
for key in transformation.keys():
for i in range(3):
self.assertAlmostEqual(decomposedInterpolatedPhi[key][i],
transformation[key][i]/2.0 ,places=3)
interpolatedPhi = spam.DIC.interpolatePhiField(fieldCoords, fieldValues, interpCoords, neighbourRadius=100)
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, -1], transformation['t'], atol=0.01))
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, 0:3], numpy.eye(3), atol=0.01))
def test_computeRigidPhi(self):
......
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