From 0e13be2d1970221b8086bfeda61bc463e9d65e55 Mon Sep 17 00:00:00 2001
From: Martin Genet <martin.genet@polytechnique.edu>
Date: Tue, 21 May 2019 10:42:37 +0200
Subject: [PATCH] Cleaning k-space resampling
---
downsample_images.py | 51 ++++++++------------
generate_images.py | 109 +++++++++++++++++++------------------------
2 files changed, 69 insertions(+), 91 deletions(-)
diff --git a/downsample_images.py b/downsample_images.py
index 20f0163..162a27a 100644
--- a/downsample_images.py
+++ b/downsample_images.py
@@ -43,25 +43,28 @@ def downsample_images(
images_ndim = myvtk.getImageDimensionality(
image=image,
verbose=0)
+ mypy.my_print(verbose, "images_ndim = "+str(images_ndim))
images_spacing = image.GetSpacing()
- images_origin = image.GetOrigin()
+ mypy.my_print(verbose, "images_spacing = "+str(images_spacing))
+ images_delta = images_spacing[0:images_ndim]
+ mypy.my_print(verbose, "images_delta = "+str(images_delta))
images_nvoxels = myvtk.getImageDimensions(
image=image,
verbose=0)
mypy.my_print(verbose, "images_nvoxels = "+str(images_nvoxels))
- assert (len(images_nvoxels) == images_ndim)
images_npoints = numpy.prod(images_nvoxels)
- # mypy.my_print(verbose, "images_npoints = "+str(images_npoints))
- assert (images_npoints == image.GetNumberOfPoints())
+ mypy.my_print(verbose, "images_npoints = "+str(images_npoints))
images_downsampled_nvoxels = numpy.array(images_nvoxels)/numpy.array(downsampling_factors)
images_downsampled_nvoxels = numpy.ceil(images_downsampled_nvoxels)
images_downsampled_nvoxels = [int(n) for n in images_downsampled_nvoxels]
mypy.my_print(verbose, "images_downsampled_nvoxels = "+str(images_downsampled_nvoxels))
- downsampling_factors = numpy.array(images_nvoxels)/numpy.array(images_downsampled_nvoxels)
+ downsampling_factors = list(numpy.array(images_nvoxels)/numpy.array(images_downsampled_nvoxels))
mypy.my_print(verbose, "downsampling_factors = "+str(downsampling_factors))
downsampling_factor = numpy.prod(downsampling_factors)
mypy.my_print(verbose, "downsampling_factor = "+str(downsampling_factor))
+ images_downsampled_delta = list(numpy.array(images_delta)*numpy.array(downsampling_factors))
+ # mypy.my_print(verbose, "images_downsampled_delta = "+str(images_downsampled_delta))
images_downsampled_npoints = numpy.prod(images_downsampled_nvoxels)
# mypy.my_print(verbose, "images_downsampled_npoints = "+str(images_downsampled_npoints))
@@ -149,32 +152,18 @@ def downsample_images(
else:
image_downsampled = vtk.vtkImageData()
- if (images_ndim == 1):
- extent = [0, images_downsampled_nvoxels[0]-1, 0, 0, 0, 0]
- elif (images_ndim == 2):
- extent = [0, images_downsampled_nvoxels[0]-1, 0, images_downsampled_nvoxels[1]-1, 0, 0]
- elif (images_ndim == 3):
- extent = [0, images_downsampled_nvoxels[0]-1, 0, images_downsampled_nvoxels[1]-1, 0, images_downsampled_nvoxels[2]-1]
- image_downsampled.SetExtent(extent)
- mypy.my_print(verbose, "extent = "+str(extent))
-
- if (images_ndim == 1):
- spacing = [images_spacing[0]*downsampling_factors[0], 1., 1.]
- elif (images_ndim == 2):
- spacing = [images_spacing[0]*downsampling_factors[0], images_spacing[1]*downsampling_factors[1], 1.]
- elif (images_ndim == 3):
- spacing = [images_spacing[0]*downsampling_factors[0], images_spacing[1]*downsampling_factors[1], images_spacing[2]*downsampling_factors[2]]
- image_downsampled.SetSpacing(spacing)
- mypy.my_print(verbose, "spacing = "+str(spacing))
-
- if (images_ndim == 1):
- origin = [spacing[0]/2, 0., 0.]
- elif (images_ndim == 2):
- origin = [spacing[0]/2, spacing[1]/2, 0.]
- elif (images_ndim == 3):
- origin = [spacing[0]/2, spacing[1]/2, spacing[2]/2]
- image_downsampled.SetOrigin(origin)
- mypy.my_print(verbose, "origin = "+str(origin))
+ dimensions_downsampled = images_downsampled_nvoxels+[1]*(3-images_ndim)
+ mypy.my_print(verbose, "dimensions_downsampled = "+str(dimensions_downsampled))
+ image_downsampled.SetDimensions(dimensions_downsampled)
+
+ spacing_downsampled = images_downsampled_delta+[1.]*(3-images_ndim)
+ mypy.my_print(verbose, "spacing_downsampled = "+str(spacing_downsampled))
+ image_downsampled.SetSpacing(spacing_downsampled)
+
+ origin_downsampled = list(numpy.array(images_downsampled_delta)/2)
+ origin_downsampled = origin_downsampled+[0.]*(3-images_ndim)
+ mypy.my_print(verbose, "origin_downsampled = "+str(origin_downsampled))
+ image_downsampled.SetOrigin(origin_downsampled)
image_downsampled_scalars = myvtk.createDoubleArray(
name="ImageScalars",
diff --git a/generate_images.py b/generate_images.py
index cde5c46..21c7de0 100644
--- a/generate_images.py
+++ b/generate_images.py
@@ -28,29 +28,29 @@ def set_I_woGrad(
image,
X,
I,
- vtk_image_scalars,
+ image_upsampled_scalars,
k_point,
G=None,
Finv=None,
- vtk_gradient_vectors=None):
+ image_gradient_upsampled_vectors=None):
image.I0(X, I)
- vtk_image_scalars.SetTuple(k_point, I)
+ image_upsampled_scalars.SetTuple(k_point, I)
def set_I_wGrad(
image,
X,
I,
- vtk_image_scalars,
+ image_upsampled_scalars,
k_point,
G,
Finv,
- vtk_gradient_vectors):
+ image_gradient_upsampled_vectors):
image.I0_wGrad(X, I, G)
- vtk_image_scalars.SetTuple(k_point, I)
+ image_upsampled_scalars.SetTuple(k_point, I)
G = numpy.dot(G, Finv)
- vtk_gradient_vectors.SetTuple(k_point, G)
+ image_gradient_upsampled_vectors.SetTuple(k_point, G)
def generate_images(
images,
@@ -90,61 +90,47 @@ def generate_images(
evolution,
generate_image_gradient)
- vtk_image = vtk.vtkImageData()
-
- delta = numpy.array(images["L"])/numpy.array(images["n_voxels"])
-
- images["n_voxels_upsampled"] = numpy.array(images["n_voxels"])*numpy.array(images["upsampling_factors"])
- delta_upsampled = numpy.array(images["L"])/numpy.array(images["n_voxels_upsampled"])
-
- if (images["n_dim"] == 1):
- extent = [0, images["n_voxels_upsampled"][0]-1, 0, 0, 0, 0]
- elif (images["n_dim"] == 2):
- extent = [0, images["n_voxels_upsampled"][0]-1, 0, images["n_voxels_upsampled"][1]-1, 0, 0]
- elif (images["n_dim"] == 3):
- extent = [0, images["n_voxels_upsampled"][0]-1, 0, images["n_voxels_upsampled"][1]-1, 0, images["n_voxels_upsampled"][2]-1]
- vtk_image.SetExtent(extent)
- mypy.my_print(verbose, "extent = "+str(extent))
-
- if (images["n_dim"] == 1):
- spacing = [delta_upsampled[0], 1., 1.]
- elif (images["n_dim"] == 2):
- spacing = [delta_upsampled[0], delta_upsampled[1], 1.]
- elif (images["n_dim"] == 3):
- spacing = [delta_upsampled[0], delta_upsampled[1], delta_upsampled[2]]
- vtk_image.SetSpacing(spacing)
- mypy.my_print(verbose, "spacing = "+str(spacing))
-
- if (images["n_dim"] == 1):
- origin = [delta[0]/2, 0., 0.]
- elif (images["n_dim"] == 2):
- origin = [delta[0]/2, delta[1]/2, 0.]
- elif (images["n_dim"] == 3):
- origin = [delta[0]/2, delta[1]/2, delta[2]/2]
- vtk_image.SetOrigin(origin)
- mypy.my_print(verbose, "origin = "+str(origin))
-
- n_points = vtk_image.GetNumberOfPoints()
- vtk_image_scalars = myvtk.createFloatArray(
+ image_upsampled = vtk.vtkImageData()
+
+ n_voxels_upsampled = list(numpy.array(images["n_voxels"])*numpy.array(images["upsampling_factors"]))
+
+ dimensions_upsampled = n_voxels_upsampled+[1]*(3-images["n_dim"])
+ mypy.my_print(verbose, "dimensions_upsampled = "+str(dimensions_upsampled))
+ image_upsampled.SetDimensions(dimensions_upsampled)
+
+ delta = list(numpy.array(images["L"])/numpy.array(images["n_voxels"]))
+ delta_upsampled = list(numpy.array(delta)/numpy.array(images["upsampling_factors"]))
+
+ spacing_upsampled = delta_upsampled+[1.]*(3-images["n_dim"])
+ mypy.my_print(verbose, "spacing_upsampled = "+str(spacing_upsampled))
+ image_upsampled.SetSpacing(spacing_upsampled)
+
+ origin_upsampled = list(numpy.array(delta_upsampled)/2)
+ origin_upsampled = origin_upsampled+[0.]*(3-images["n_dim"])
+ mypy.my_print(verbose, "origin_upsampled = "+str(origin_upsampled))
+ image_upsampled.SetOrigin(origin_upsampled)
+
+ n_points_upsampled = image_upsampled.GetNumberOfPoints()
+ image_upsampled_scalars = myvtk.createFloatArray(
name="ImageScalars",
n_components=1,
- n_tuples=n_points,
+ n_tuples=n_points_upsampled,
verbose=verbose-1)
- vtk_image.GetPointData().SetScalars(vtk_image_scalars)
+ image_upsampled.GetPointData().SetScalars(image_upsampled_scalars)
if (generate_image_gradient):
- vtk_gradient = vtk.vtkImageData()
- vtk_gradient.DeepCopy(vtk_image)
+ image_gradient_upsampled = vtk.vtkImageData()
+ image_gradient_upsampled.DeepCopy(image_upsampled)
- vtk_gradient_vectors = myvtk.createFloatArray(
+ image_gradient_upsampled_vectors = myvtk.createFloatArray(
name="ImageScalarsGradient",
n_components=3,
- n_tuples=n_points,
+ n_tuples=n_points_upsampled,
verbose=verbose-1)
- vtk_gradient.GetPointData().SetScalars(vtk_gradient_vectors)
+ image_gradient_upsampled.GetPointData().SetScalars(image_gradient_upsampled_vectors)
else:
- vtk_gradient = None
- vtk_gradient_vectors = None
+ image_gradient_upsampled = None
+ image_gradient_upsampled_vectors = None
x = numpy.empty(3)
X = numpy.empty(3)
@@ -164,23 +150,23 @@ def generate_images(
t = images["T"]*float(k_frame)/(images["n_frames"]-1) if (images["n_frames"]>1) else 0.
mypy.my_print(verbose, "t = "+str(t))
mapping.init_t(t)
- for k_point in xrange(n_points):
- vtk_image.GetPoint(k_point, x)
+ for k_point in xrange(n_points_upsampled):
+ image_upsampled.GetPoint(k_point, x)
#print "x0 = "+str(x)
mapping.X(x, X, Finv)
#print "X = "+str(X)
- set_I(image, X, I, vtk_image_scalars, k_point, G, Finv, vtk_gradient_vectors)
+ set_I(image, X, I, image_upsampled_scalars, k_point, G, Finv, image_gradient_upsampled_vectors)
global_min = min(global_min, I[0])
global_max = max(global_max, I[0])
- #print vtk_image
+ #print image_upsampled
myvtk.writeImage(
- image=vtk_image,
+ image=image_upsampled,
filename=images["folder"]+"/"+images["basename"]+"_"+str(k_frame).zfill(images["zfill"])+"."+images["ext"],
verbose=verbose-1)
if (generate_image_gradient):
- #print vtk_gradient
+ #print image_gradient_upsampled
myvtk.writeImage(
- image=vtk_gradient,
+ image=image_gradient_upsampled,
filename=images["folder"]+"/"+images["basename"]+"-grad"+"_"+str(k_frame).zfill(images["zfill"])+"."+images["ext"],
verbose=verbose-1)
# mypy.my_print(verbose, "global_min = "+str(global_min))
@@ -190,6 +176,7 @@ def generate_images(
downsampling = False
else:
downsampling = True
+
ddic.downsample_images(
images_folder=images["folder"],
images_basename=images["basename"],
@@ -199,8 +186,10 @@ def generate_images(
verbose=verbose)
if (images["data_type"] in ("float")):
- pass
+ normalizing = False
elif (images["data_type"] in ("unsigned char", "unsigned short", "unsigned int", "unsigned long", "unsigned float", "uint8", "uint16", "uint32", "uint64", "ufloat")):
+ normalizing = True
+
ddic.normalize_images(
images_folder=images["folder"],
images_basename=images["basename"]+("_downsampled"*(downsampling)),
--
GitLab