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 390043dc authored by Olga Stamati's avatar Olga Stamati
Browse files

multiprocessing in applyPhiField()

parent 201a8ecf
Pipeline #50861 passed with stages
in 24 minutes and 1 second
......@@ -162,7 +162,7 @@ def applyPhiPython(im, Phi=None, PhiPoint=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):
def applyPhiField(im, fieldName=None, fieldCoords=None, fieldValues=None, fieldBinRatio=1.0, neighbours=8, interpolationOrder=3, verbose=False, numberOfThreads=1):
"""
Deform a 3D image using a field of deformation functions "Phi" coming from a regularGrid,
applied using scipy.ndimage.map_coordinates.
......@@ -197,12 +197,22 @@ def applyPhiField(im, fieldName=None, fieldCoords=None, fieldValues=None, fieldB
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
numberOfThreads : int, optional
Number of Threads for multiprocessing
Default = 1
Returns
-------
imDef : 3D array
deformed greylevels by a field of deformation functions "Phi"
"""
import spam.deformation
import multiprocessing
import progressbar
# Create the grid of the input image
imSize = im.shape
......@@ -216,9 +226,7 @@ def applyPhiField(im, fieldName=None, fieldCoords=None, fieldValues=None, fieldB
coordIn[:, 1] = coordinates_mgrid[1].ravel()
coordIn[:, 2] = coordinates_mgrid[2].ravel()
numberofPoints = imSize[0] * imSize[1] * imSize[2]
# Copy initial coordinates to the deformed coordinates
# Initialise deformed coordinates
coordDef = coordIn.copy()
# Read input PhiField, usually the result of a regularGrid correlation
......@@ -236,48 +244,96 @@ def applyPhiField(im, fieldName=None, fieldCoords=None, fieldValues=None, fieldB
tree = KDTree(fieldCoords)
# Loop over each point of the grid of the input image
for point in range(coordIn.shape[0]):
if verbose:
print("\rWorking on point {} {}%".format(point, (point / float(numberofPoints)) * 100), end='')
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()
# 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)
#print("coordNew", coordNew)
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")
# Calculate the distance of the current point to the points of the input Phi field
distance, indices = tree.query(coordIn[point][0:3], k=neighbours)
#print("\nMaster: Launching workers")
for thread in range(numberOfThreads):
p = multiprocessing.Process(target=worker, args=(thread, qJobs, qResults, ))
p.start()
# Check if we've hit the same point
if numpy.any(distance == 0):
if verbose:
pbar = progressbar.ProgressBar(maxval=numberofPoints).start()
# Deform the coordinates of the current point
# by subtracting the translation part of the deformation function Phi
coordDef[point][:3] -= fieldValues[indices][numpy.where(distance == 0)][0][0:3, -1].copy()
finishedThreads = 0
finishedJobs = 0
# Check if we have asked only for the closest neighbour
elif neighbours == 1:
while finishedThreads < numberOfThreads:
result = qResults.get()
# Deform the coordinates of the current point
# by subtracting the translation part of the deformation function Phi
# applied on the current point
coordDef[point][:3] -= spam.deformation.decomposePhi(fieldValues[indices].copy(),
PhiCentre=fieldCoords[indices],
PhiPoint=coordIn[point][:3])["t"]
if result == "STOP":
finishedThreads += 1
#print("\tNumber of finished threads = ", finishedThreads)
# 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
coordDef[point][:3] -= numpy.dot(spam.deformation.decomposePhi(fieldValues[indices][neighbour].copy(),
PhiCentre=fieldCoords[indices][neighbour],
PhiPoint=coordIn[point][:3])["t"],
weightInv)
#print("Master: got {}".format(result))
finishedJobs += 1
coordDef[result[0], :3] = result[1]
if verbose:
pbar.update(finishedJobs)
if verbose:
pbar.finish()
# Deform the image
imDef = numpy.zeros_like(im, dtype='<f4')
......
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