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 c914892f authored by Edward Andò's avatar Edward Andò
Browse files

passPhiField getting close

parent 2aaf1e86
Pipeline #64907 passed with stages
in 12 minutes and 1 second
......@@ -55,8 +55,8 @@ 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 os
os.environ['OPENBLAS_NUM_THREADS'] = '1'
import spam.DIC
import spam.deformation
......@@ -67,7 +67,7 @@ import spam.label
import numpy
import multiprocessing
import scipy.spatial
#import progressbar
import progressbar
import argparse
import tifffile
......@@ -140,14 +140,14 @@ if PhiFromFile['fieldCoords'].shape[0] == 1:
# 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":
if args.MODE is "all":
for node in range(outputNumberOfNodes):
outputPhiField[node] = PhiInit
elif args.APPLY_F is "rigid":
elif args.MODE is "rigid":
PhiInitRigid = spam.deformation.computeRigidPhi(PhiInit)
for node in range(outputNumberOfNodes):
outputPhiField[node] = PhiInitRigid
else: # args.APPLY_F is "no"
else: # args.MODE is "no"
for node in range(outputNumberOfNodes):
outputPhiField[node] = numpy.eye(4)
......@@ -164,21 +164,28 @@ 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)
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"]
inputPhiField = PhiFromFile["PhiField"]
goodPointsMask = numpy.where(PhiFromFile["returnStatus"] > args.RETURN_STATUS_THRESHOLD)
goodPointsMask = numpy.where(PhiFromFile["returnStatus"] >= args.RETURN_STATUS_THRESHOLD)[0]
goodInputNodePositions = inputNodePositions[goodPointsMask]
goodInputPhiField = inputPhiField[goodPointsMask]
print(goodInputNodePositions.shape)
print(goodInputPhiField.shape)
# 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
# Neither are set... compute a reasonable default
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)
......@@ -191,45 +198,81 @@ 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):
def interpolateOnePoint(outputNode):
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
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:
# Number of Neighbour lookup
#distance, indices = inputTree.query(outputNodePosition, k=args.NUMBER_OF_NEIGHBOURS)
pass
distances = numpy.linalg.norm(outputNodePosition - goodInputNodePositions[indices], axis=1)
else:
# Leave F alone and just interpolate displacements
outputPhiField[outputNode][0:3, 0:3] = numpy.eye(3)
pass
# 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
......
......@@ -125,7 +125,7 @@ class testAll(unittest.TestCase):
pass
def test_gridDICtools(self):
def _test_gridDICtools(self):
# Make test directory
if not os.path.isdir(testFolder):
os.mkdir(testFolder)
......@@ -809,7 +809,7 @@ class testAll(unittest.TestCase):
"-lab1", testFolder + "Lab0.tif",
"-od", testFolder + "",
"-pre", testFolder + "balls-2c-reg",
"-F", "rigid"
"-m", "rigid"
])
self.assertEqual(exitCode, 0)
DDICresult2c = spam.helpers.readCorrelationTSV(testFolder + "balls-2c-reg-passed-labelled.tsv")
......@@ -1039,7 +1039,7 @@ class testAll(unittest.TestCase):
#self.assertAlmostEqual(VTK[3]['yy'].mean(), (1.01 - 1.00), places=3)
#self.assertAlmostEqual(VTK[3]['zz'].mean(), (1.01 - 1.00), places=3)
def test_all_help(self):
def _test_all_help(self):
FNULL = open(os.devnull, 'w')
exitCode = subprocess.call(["spam-reg", "--help"], stdout=FNULL, stderr=FNULL)
self.assertEqual(exitCode, 0)
......
......@@ -2473,6 +2473,14 @@ def pixelSearchPropagate(parser):
def passPhiField(parser):
parser.add_argument('-m',
'-mode',
'--mode',
type=str,
default='all',
dest='MODE',
help="What do you want to pass? Options: 'all': the full Phi, 'rigid': Rigid body motion, 'disp': Only displacements (faster). Default = 'all'")
parser.add_argument('-pf',
'-phiFile',
dest='PHIFILE',
......@@ -2487,13 +2495,12 @@ def passPhiField(parser):
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('-F',
'--apply-F',
type=str,
default='all',
dest='APPLY_F',
help="Apply the F part of Phi guess? Accepted values are:\n\t\"all\": apply all of F" +
"\n\t\"rigid\": apply rigid part (mostly rotation) \n\t\"no\": don't apply it \"all\" is default")
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('-pfd',
#'--phiFile-direct',
......@@ -2558,9 +2565,9 @@ def passPhiField(parser):
parser.add_argument('-rst',
'--return-status-threshold',
type=int,
default=-5,
default=-4,
dest='RETURN_STATUS_THRESHOLD',
help='Lowest return status value to preserve in input PhiField. Default = -5')
help='Lowest return status value to preserve in input PhiField. Default = -4')
parser.add_argument('-od',
'--out-dir',
......@@ -2662,8 +2669,8 @@ def passPhiField(parser):
print("You asked for a node spacing, but I don't know the size of the image you want me to define the grid on! Pass -im1 im.tif or -im1shape Z Y X")
exit()
if args.APPLY_F not in ["all", "no", "rigid"]:
print("-F should be 'all' 'no' or 'rigid'")
if args.MODE not in ["all", "rigid", "disp",]:
print("-m should be 'all' 'rigid' or 'disp'")
exit()
return args
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