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

[skip-ci] applyPhiField() with disp interpolation and applyPhi mode working,...

[skip-ci] applyPhiField() with disp interpolation and applyPhi mode working, latter not parallel yet
parent 9d39910c
......@@ -18,6 +18,8 @@ this program. If not, see <http://www.gnu.org/licenses/>.
import numpy
import scipy.ndimage
import multiprocessing
nProcessesDefault = multiprocessing.cpu_count()
#numpy.set_printoptions(precision=3, suppress=True)
###########################################################
......@@ -161,7 +163,7 @@ def applyPhiPython(im, Phi=None, PhiCentre=None, interpolationOrder=3):
###############################################################
# Take a field of Phi and apply it (quite slowly) to an image
###############################################################
def applyPhiField(im, fieldName=None, fieldCoords=None, fieldValues=None, fieldBinRatio=1.0, neighbours=8, interpolationOrder=3, verbose=False, numberOfThreads=1):
def applyPhiField(im, fieldCoordsRef, PhiField, displacementMode="interpolate", nNeighbours=8, interpolationOrder=1, nProcesses=nProcessesDefault, verbose=False):
"""
Deform a 3D image using a field of deformation functions "Phi" coming from a regularGrid,
applied using scipy.ndimage.map_coordinates.
......@@ -171,187 +173,279 @@ def applyPhiField(im, fieldName=None, fieldCoords=None, fieldValues=None, fieldB
im : 3D array
3D array of grey levels to be deformed
fieldName : string, optional
Name of the file containing the deformation functions field
fieldCoords: 2D array, optional
fieldCoordsRef: 2D array, optional
nx3 array of n points coordinates (ZYX)
centre where each deformation function "Phi" has been calculated
fieldValues: 3D array, optional
PhiField: 3D array, optional
nx4x4 array of n points deformation functions
fieldBinRatio : int, optional
If the input field refers to a binned version of the image
`e.g.`, if ``fieldBinRatio = 2`` the ``fieldName`` values have been calculated
for an image half the size of the input image ``im``
Default = 1
displacementMode : string, optional
How do you want to calculate displacements?
With "interpolate" they are just interpolated from the PhiField
With "applyPhi" each neighbour
neighbours : int, optional
Neighbours for field interpolation
If == 1, the nearest neighbour is used, if >1 neighbours are weighted according to distance.
nNeighbours : int, optional
Number of the nearest neighbours to consider
#This OR neighbourRadius must be set.
Default = 8
interpolationOrder : int, optional
Order of image interpolation to use. This value is passed directly to ``scipy.ndimage.map_coordinates`` as "order".
Default = 1
verbose : boolean, optional
If True the evolution of the process is printed
Default = False
nProcesses : integer, optional
Number of processes for multiprocessing
Default = number of CPUs in the system
numberOfThreads : int, optional
Number of Threads for multiprocessing
Default = 1
verbose : boolean, optional
Print progress?
Default = True
Returns
-------
imDef : 3D array
deformed greylevels by a field of deformation functions "Phi"
"""
import spam.DIC
import spam.deformation
import multiprocessing
import progressbar
tol = 1e-6 # OS is responsible for the validity of this magic number
print("making pixel grid")
# Create the grid of the input image
imSize = im.shape
coordinates_mgrid = numpy.mgrid[0:imSize[0],
0:imSize[1],
0:imSize[2]]
coordIn = numpy.ones((imSize[0] * imSize[1] * imSize[2], 4))
coordIn[:, 0] = coordinates_mgrid[0].ravel()
coordIn[:, 1] = coordinates_mgrid[1].ravel()
coordIn[:, 2] = coordinates_mgrid[2].ravel()
pixCoordsDef = numpy.ones((imSize[0] * imSize[1] * imSize[2], 3))
pixCoordsDef[:, 0] = coordinates_mgrid[0].ravel()
pixCoordsDef[:, 1] = coordinates_mgrid[1].ravel()
pixCoordsDef[:, 2] = coordinates_mgrid[2].ravel()
print(pixCoordsDef.shape)
# Initialise deformed coordinates
coordDef = coordIn.copy()
# Read input PhiField, usually the result of a regularGrid correlation
if fieldName:
import spam.helpers.tsvio
PhiFromFile = spam.helpers.tsvio.readCorrelationTSV(fieldName, fieldBinRatio=fieldBinRatio)
fieldCoords = PhiFromFile["fieldCoords"]
fieldValues = PhiFromFile["PhiField"]
else:
fieldCoords = fieldCoords
fieldValues = fieldValues
# Create the k-d tree of the coordinates of the input Phi field
from scipy.spatial import KDTree
tree = KDTree(fieldCoords)
# Loop over each point of the grid of the input image
def worker(workerNumber, qJobs, qResults):
while True:
job = qJobs.get()
if job == "STOP":
qResults.put("STOP")
break
else:
point = job
coordNew = coordIn[point, :3].copy()
#print("\nWorking on point {}".format(point, end=''))
# Calculate the distance of the current point to the points of the input Phi field
distance, indices = tree.query(coordIn[point, :3], k=neighbours)
# Check if we've hit the same point
if numpy.any(distance == 0):
# Deform the coordinates of the current point
# by subtracting the translation part of the deformation function Phi
#coordNew -= fieldValues[indices][numpy.where(distance == 0)][0][0:3, -1].copy()
try:
coordNew += numpy.linalg.inv(fieldValues[indices][numpy.where(distance == 0)][0])[0:3, -1]
except:
coordNew = numpy.zeros(3)
## Check if we have asked only for the closest neighbour
#elif neighbours == 1:
## Deform the coordinates of the current point
## by subtracting the translation part of the deformation function Phi
## applied on the current point
#coordNew -= spam.deformation.decomposePhi(fieldValues[indices].copy(),
#PhiCentre=fieldCoords[indices],
#PhiPoint=coorXdIn[point, :3])["t"]
# Consider the k 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(neighbours):
# Calculate its weight
weightInv = (1/distance[neighbour]) / float(weightSumInv)
# Deform the coordinates of the current point
# by subtracting the translation part of the deformation function Phi
# applied on the current point
# multiplied by the weight of each neighbour
#coordNew -= numpy.dot(spam.deformation.decomposePhi(fieldValues[indices][neighbour].copy(),
#PhiCentre=fieldCoords[indices][neighbour],
#PhiPoint=coordIn[point, :3])["t"],
#weightInv)
try:
coordNew += numpy.dot(spam.deformation.decomposePhi(numpy.linalg.inv(fieldValues[indices][neighbour].copy()),
PhiCentre=fieldCoords[indices][neighbour],
PhiPoint=coordIn[point, :3])["t"],
weightInv)
except:
pass
qResults.put([point, coordNew])
numberofPoints = imSize[0] * imSize[1] * imSize[2]
qJobs = multiprocessing.Queue()
qResults = multiprocessing.Queue()
#print("\nMaster: Adding jobs to queues")
for point in range(numberofPoints):
qJobs.put(point)
for thread in range(numberOfThreads):
qJobs.put("STOP")
#print("\nMaster: Launching workers")
for thread in range(numberOfThreads):
p = multiprocessing.Process(target=worker, args=(thread, qJobs, qResults, ))
p.start()
if verbose:
pbar = progressbar.ProgressBar(maxval=numberofPoints).start()
finishedThreads = 0
finishedJobs = 0
while finishedThreads < numberOfThreads:
result = qResults.get()
if result == "STOP":
finishedThreads += 1
#print("\tNumber of finished threads = ", finishedThreads)
else:
#print("Master: got {}".format(result))
finishedJobs += 1
coordDef[result[0], :3] = result[1]
if verbose:
pbar.update(finishedJobs)
fieldCoordsDef = fieldCoordsRef + PhiField[:, 0:3, -1]
#print(fieldCoordsDef.shape)
#print(fieldCoordsDef)
print("done")
mask = numpy.isfinite(PhiField[:, 0, 0])
if displacementMode == "interpolate":
"""
In this mode we're only taking into account displacements.
We use interpolatePhiField in the deformed configuration, in displacements only,
and we don't feed PhiInv, but only the negative of the displacements
"""
print("sum(mask) = ", mask.sum())
#print(PhiField[mask].shape)
backwardsDisplacementsPhi = numpy.zeros_like(PhiField)
backwardsDisplacementsPhi[:, 0:3, -1] = -1*PhiField[:, 0:3, -1]
pixDispsDef = spam.DIC.interpolatePhiField(fieldCoordsDef[mask],
backwardsDisplacementsPhi[mask],
pixCoordsDef,
nNeighbours=nNeighbours,
interpolateF="no",
nProcesses=nProcesses,
verbose=verbose)
pixCoordsRef = pixCoordsDef + pixDispsDef[:, 0:3, -1]
elif displacementMode == "applyPhi":
"""
In this mode we're NOT interpolating the displacement field.
For each pixel, we're applying the neighbouring Phis and looking
at the resulting displacement of the pixel.
Those different displacements are weighted as a function of distance
and averaged into the point's final displacement.
Obviously if your PhiField is only a displacement field, this changes
nothing from above (except for computation time), but with some stretches
this can become interesting.
"""
print("inversing PhiField")
PhiFieldInv = numpy.zeros_like(PhiField)
for n in range(PhiField.shape[0]):
try:
PhiFieldInv[n] = numpy.linalg.inv(PhiField[n])
except:
mask[n] = False
print("done")
# mask everything
PhiFieldInvMasked = PhiFieldInv[mask]
fieldCoordsRefMasked = fieldCoordsRef[mask]
fieldCoordsDefMasked = fieldCoordsDef[mask]
# build KD-tree
treeCoordDef = scipy.spatial.KDTree(fieldCoordsDefMasked)
pixCoordsRef = numpy.zeros_like(pixCoordsDef, dtype='f4')
# for each group of pixels (divide mask==True into nProcesses chunks):
for pixNumber, pixCoordDef in enumerate(pixCoordsDef):
# get nNeighbours and compute distance weights
distances, indices = treeCoordDef.query(pixCoordDef, k=nNeighbours)
weights = (1/(distances+tol))
displacement = numpy.zeros(3, dtype='float')
#print()
# for each neighbour
for neighbour, index in enumerate(indices):
# apply PhiInv to current point with PhiCentre = fieldCoordsREF <- this is important
# -> this gives a displacement for each neighbour
PhiInv = PhiFieldInvMasked[index]
#print("PhiInv", PhiInv)
translationTmp = PhiInv[0:3, -1].copy()
dist = pixCoordDef - fieldCoordsRefMasked[index]
#print(f"dist = {dist}")
translationTmp -= dist - numpy.dot(PhiInv[0:3, 0:3], dist)
#print(f"translationTmp ({neighbour}): {translationTmp} (weight = {weights[neighbour]})")
displacement += translationTmp*weights[neighbour]
# compute resulting displacement as weights * displacements
# compute pixel coordinates in reference config
#print("pixCoordDef", pixCoordDef)
#print("displacement", displacement)
pixCoordsRef[pixNumber] = pixCoordDef + displacement / weights.sum()
#pixCoordsRef[pixNumber] = pixCoordDef + numpy.sum(displacements*weights[:, numpy.newaxis], axis=0)
#global inter
#def dispOnePoint(defCoord):
## Check if we've hit the same point
#if numpy.any(distance == 0):
## Deform the coordinates of the current point
## by subtracting the translation part of the deformation function Phi
##coordNew -= fieldValues[indices][numpy.where(distance == 0)][0][0:3, -1].copy()
#try:
#coordNew += numpy.linalg.inv(fieldValues[indices][numpy.where(distance == 0)][0])[0:3, -1]
#except:
#coordNew = numpy.zeros(3)
### Check if we have asked only for the closest neighbour
##elif neighbours == 1:
### Deform the coordinates of the current point
### by subtracting the translation part of the deformation function Phi
### applied on the current point
##coordNew -= spam.deformation.decomposePhi(fieldValues[indices].copy(),
##PhiCentre=fieldCoords[indices],
##PhiPoint=coorXdIn[point, :3])["t"]
## Consider the k 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(neighbours):
## Calculate its weight
#weightInv = (1/distance[neighbour]) / float(weightSumInv)
## Deform the coordinates of the current point
## by subtracting the translation part of the deformation function Phi
## applied on the current point
## multiplied by the weight of each neighbour
##coordNew -= numpy.dot(spam.deformation.decomposePhi(fieldValues[indices][neighbour].copy(),
##PhiCentre=fieldCoords[indices][neighbour],
##PhiPoint=pixCoordDef[point, :3])["t"],
##weightInv)
#try:
#coordNew += numpy.dot(spam.deformation.decomposePhi(numpy.linalg.inv(fieldValues[indices][neighbour].copy()),
#PhiCentre=fieldCoords[indices][neighbour],
#PhiPoint=pixCoordDef[point, :3])["t"],
#weightInv)
#except:
#pass
#qResults.put([point, coordNew])
## Iterate through flat field of Fs
#if verbose:
#pbar = progressbar.ProgressBar(maxval=pointsToEstimate.shape[0]).start()
#finishedPoints = 0
## Run multiprocessing filling in FfieldFlatGood, which will then update FfieldFlat
#with multiprocessing.Pool(processes=nProcesses) as pool:
#for returns in pool.imap_unordered(dispOnePoint, range(pointsToEstimate.shape[0])):
#if verbose:
#finishedPoints += 1
#pbar.update(finishedPoints)
#estimatedDisplacements[returns[0]] = returns[1]
#pool.close()
#pool.join()
#if verbose: pbar.finish()
## overwrite bad points displacements
#return estimatedDisplacements
## Loop over each point of the grid of the input image
#def worker(workerNumber, qJobs, qResults):
#while True:
#job = qJobs.get()
#if job == "STOP":
#qResults.put("STOP")
#break
#else:
#numberofPoints = imSize[0] * imSize[1] * imSize[2]
#qJobs = multiprocessing.Queue()
#qResults = multiprocessing.Queue()
##print("\nMaster: Adding jobs to queues")
#for point in range(numberofPoints):
#qJobs.put(point)
#for thread in range(numberOfThreads):
#qJobs.put("STOP")
##print("\nMaster: Launching workers")
#for thread in range(numberOfThreads):
#p = multiprocessing.Process(target=worker, args=(thread, qJobs, qResults, ))
#p.start()
#if verbose:
#pbar = progressbar.ProgressBar(maxval=numberofPoints).start()
#finishedThreads = 0
#finishedJobs = 0
#while finishedThreads < numberOfThreads:
#result = qResults.get()
#if result == "STOP":
#finishedThreads += 1
##print("\tNumber of finished threads = ", finishedThreads)
#else:
##print("Master: got {}".format(result))
#finishedJobs += 1
#coordDef[result[0], :3] = result[1]
#if verbose:
#pbar.update(finishedJobs)
if verbose:
pbar.finish()
#if verbose:
#pbar.finish()
# Deform the image
imDef = numpy.zeros_like(im, dtype='<f4')
#imDef = numpy.zeros_like(im, dtype='<f4')
imDef += scipy.ndimage.map_coordinates(im,
coordDef[:, 0:3].T,
mode="constant",
order=interpolationOrder).reshape(imDef.shape).astype('<f4')
imDef = scipy.ndimage.map_coordinates(im,
pixCoordsRef.T,
mode="constant",
order=interpolationOrder).reshape(im.shape).astype('<f4')
return imDef
......
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