Vous avez reçu un message "Your GitLab account has been locked ..." ? Pas d'inquiétude : lisez cet article https://docs.gricad-pages.univ-grenoble-alpes.fr/help/unlock/

Commit 399dd662 authored by Edward Andò's avatar Edward Andò
Browse files

adding spam-passPhiField

parent d544b6b5
Pipeline #64502 passed with stages
in 12 minutes and 49 seconds
#!/usr/bin/env python
"""
This script manipulates different Phi fields:
Single Phi values:
- spam-ereg
- spam-reg
- spam-mmr
- spam-mmr-graphical
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:
- apply a registration (single Phi) to a series of points:
- defined on a grid with NS
- or as centres-of-mass of labelled images
- apply an existing Phi-field to a new basis:
- spam-ldic result onto grid with finer NS
- spam-ldic onto centres-of-mass of labels
- spam-ddic result onto grid
- merge fields on different grids
- subtract kinematics field on the same basis
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-passPhiField "+spam.helpers.optionsParser.GLPv3descriptionHeader +\
"This script facilitates the passing of Phi fields onto different bases\n"+
"[reg, grid, discrete] -> [grid, discrete]",
formatter_class=argparse.RawTextHelpFormatter)
# Parse arguments with external helper function
args = spam.helpers.optionsParser.passPhiField(parser)
print("spam-passPhiField -- Current Settings:")
argsDict = vars(args)
for key in sorted(argsDict):
print("\t{}: {}".format(key, argsDict[key]))
###############################################################
### Step 0 define OUTPUT node positions -- either grid or labels:
###############################################################
if args.LAB1 is not None:
lab1 = tifffile.imread(args.LAB1.name).astype(spam.label.labelType)
boundingBoxes = spam.label.boundingBoxes(lab1)
outputNodePositions = spam.label.centresOfMass(lab1, boundingBoxes=boundingBoxes)
outputNumberOfNodes = outputNodePositions.shape[0]
### Otherwise we are in node spacing and half-window size mode
else:
outputNodePositions, outputNodesDim = spam.DIC.makeGrid(args.im1shape, args.NS)
# start setting up
outputNumberOfNodes = outputNodePositions.shape[0]
# Now that we know how many points we want to correlate, initalise PhiField
outputPhiField = numpy.zeros((outputNumberOfNodes, 4, 4))
###############################################################
### Step 1 (mandatory) read input Phi File
###############################################################
PhiFromFile = spam.helpers.readCorrelationTSV(args.PHIFILE.name, fieldBinRatio=args.PHIFILE_BIN_RATIO)
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:
PhiInit = PhiFromFile['PhiField'][0]
print(f"\tI read a registration from a file in binning {args.PHIFILE_BIN_RATIO}")
# In the special case of a registration, initialise output variables:
pixelSearchCC = numpy.zeros((outputNumberOfNodes), dtype=float)
error = numpy.zeros((outputNumberOfNodes), dtype=float)
returnStatus = numpy.ones((outputNumberOfNodes), dtype=int)
deltaPhiNorm = numpy.ones((outputNumberOfNodes), dtype=int)
iterations = numpy.ones((outputNumberOfNodes), dtype=int)
decomposedPhiInit = spam.deformation.decomposePhi(PhiInit)
print("\tTranslations (px)")
print("\t\t", decomposedPhiInit['t'])
print("\tRotations (deg)")
print("\t\t", decomposedPhiInit['r'])
print("\tZoom")
print("\t\t", decomposedPhiInit['z'])
del decomposedPhiInit
# We have a registration to apply to all points.
# This is done in 2 steps:
# 1. by copying the registration's little F to the Fs of all points
# 2. by calling the decomposePhi function to compute the translation of each point
if args.APPLY_F is "all":
for node in range(outputNumberOfNodes):
outputPhiField[node] = PhiInit
elif args.APPLY_F is "rigid":
PhiInitRigid = spam.deformation.computeRigidPhi(PhiInit)
for node in range(outputNumberOfNodes):
outputPhiField[node] = PhiInitRigid
else: # args.APPLY_F is "no"
for node in range(outputNumberOfNodes):
outputPhiField[node] = numpy.eye(4)
# In any case, apply the displacements to the points
for node in range(outputNumberOfNodes):
outputPhiField[node][0:3, -1] = spam.deformation.decomposePhi(PhiInit,
PhiCentre=PhiFromFile["fieldCoords"][0],
PhiPoint=outputNodePositions[node])["t"]
###############################################################
### Input Phi file is a Phi FIELD
###############################################################
else:
# Interpolate these?
pixelSearchCC = numpy.zeros((outputNumberOfNodes), dtype=float)
error = numpy.zeros((outputNumberOfNodes), dtype=float)
returnStatus = numpy.ones((outputNumberOfNodes), dtype=int)
deltaPhiNorm = numpy.ones((outputNumberOfNodes), dtype=int)
iterations = numpy.ones((outputNumberOfNodes), dtype=int)
# We don't trust this completely, re-interpolate it onto the grid
# Read the coordinates and values of the input F field
inputNodePositions = PhiFromFile["fieldCoords"]
inputPhiField = PhiFromFile["PhiField"]
goodPointsMask = numpy.where(PhiFromFile["returnStatus"] > args.RETURN_STATUS_THRESHOLD)
# 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 args.LAB1 is None:
args.NEIGHBOUR_RADIUS = 2*int(numpy.mean(args.NS)//1)
print(f"Neither number of neighbours nor neighbour distance set, using default distance of 2*mean(NS) = {args.NEIGHBOUR_RADIUS}")
else:
# Come up with a good default radius size
args.NEIGHBOUR_RADIUS = 5*numpy.mean(spam.label.equivalentRadii(lab1, boundingBoxes=boundingBoxes)[1:])
print(f"Neither number of neighbours nor neighbour distance set, using default distance of 5*mean particle radius = {args.NEIGHBOUR_RADIUS}")
#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
goodInputNodePositions = inputNodePositions[goodPointsMask]
inputTree = scipy.spatial.KDTree(goodInputNodePositions)
# Loop over each point of the current grid
for outputNode in range(outputNumberOfNodes):
outputNodePosition = outputNodePositions[outputNode]
if args.APPLY_F in ["all", "rigid"]:
# reset F to zero, since we'll be doing an additive interpolation at the bottom here
outputPhiField[outputNode][0:3, 0:3] = 0
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!")
pass
else:
# average input and fill in output
#print(outputNodePosition)
#print(goodInputNodePositions[indices])
distances = numpy.linalg.norm(outputNodePosition - goodInputNodePositions[indices], axis=1)
#print(distances)
weights = (1/(distances+tol))
weightsTotal = float(weights.sum())
# TODO:
#print(weights)
#print(weights)
#print(weightsTotal)
for neighbourNumber, neighbourIndex in enumerate(indices):
weightInv = weights[neighbourNumber] / weightsTotal
outputPhiField[outputNode][:-1] += inputPhiField[neighbourIndex][:-1]*weightInv
else:
# Number of Neighbour lookup
#distance, indices = inputTree.query(outputNodePosition, k=args.NUMBER_OF_NEIGHBOURS)
pass
else:
# Leave F alone and just interpolate displacements
outputPhiField[outputNode][0:3, 0:3] = numpy.eye(3)
pass
## 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
#- apply a registration (single Phi) to a series of points:
#- defined on a grid with NS
#- or as centres-of-mass of labelled images
#- apply an existing Phi-field to a new basis:
#- spam-ldic result onto grid with finer NS
#- spam-ldic onto centres-of-mass of labels
#- spam-ddic result onto grid
#- merge fields on different grids
#- subtract kinematics field on the same basis
#Outputs are:
#- TSV files
#- (optional) VTK files for visualisation
#- (optional) TIF files in the case of gridded data
if args.LAB1 is not None:
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(range(outputNumberOfNodes)),
outputNodePositions[:, 0], outputNodePositions[:, 1], outputNodePositions[:, 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],
pixelSearchCC,
returnStatus,
error,
deltaPhiNorm,
iterations]).T
numpy.savetxt(args.OUT_DIR+"/"+args.PREFIX+".tsv",
outMatrix,
fmt='%.7f',
delimiter='\t',
newline='\n',
comments='',
header=TSVheader)
if args.TIFF:
if args.LAB1 == None:
if outputNodesDim[0] != 1:
tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-Zdisp.tif", outputPhiField[:, 0, -1].astype('<f4').reshape(outputNodesDim))
tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-Ydisp.tif", outputPhiField[:, 1, -1].astype('<f4').reshape(outputNodesDim))
tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-Xdisp.tif", outputPhiField[:, 2, -1].astype('<f4').reshape(outputNodesDim))
#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 and args.LAB1 is None:
cellData = {}
cellData['displacements'] = outputPhiField[:, :-1, 3].reshape((outputNodesDim[0], outputNodesDim[1], outputNodesDim[2], 3))
#cellData['pixelSearchCC'] = pixelSearchCC.reshape(outputNodesDim)
# 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
args.HWS = numpy.array(args.NS)/2
spam.helpers.writeStructuredVTK(origin=outputNodePositions[0]-args.HWS, aspectRatio=args.NS, cellData=cellData, fileName=args.OUT_DIR+"/"+args.PREFIX+".vtk")
elif args.VTK and args.LAB1 is not None:
# Redundant output for VTK visualisation
magDisp = numpy.zeros(outputNumberOfNodes)
for node in range(outputNumberOfNodes):
magDisp[node] = numpy.linalg.norm(outputPhiField[node][0:3,-1])
VTKglyphDict = {'displacements': outputPhiField[:, 0:3, -1],
'mag(displacements)': magDisp,
#'pixelSearchCC': pixelSearchCC
}
spam.helpers.writeGlyphsVTK(outputNodePositions, VTKglyphDict, fileName=args.OUT_DIR + "/" + args.PREFIX + ".vtk")
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