From 4d3758593aaf9a2b7a717bf1c2b006dea6ad70bb Mon Sep 17 00:00:00 2001
From: Martin Genet <martin.genet@polytechnique.edu>
Date: Mon, 20 May 2019 16:32:14 +0200
Subject: [PATCH] k-space resampling with and without discretization change
---
__init__.py | 4 +-
downsample_images.py | 275 ++++++++++++++++++++++++++++++++
generate_images.py | 64 ++++----
generate_undersampled_images.py | 6 +-
normalize_images.py | 18 +--
resample_images.py | 151 ------------------
6 files changed, 327 insertions(+), 191 deletions(-)
create mode 100644 downsample_images.py
delete mode 100644 resample_images.py
diff --git a/__init__.py b/__init__.py
index 089a12d..3ac733e 100644
--- a/__init__.py
+++ b/__init__.py
@@ -1,5 +1,6 @@
#this file was generated by __init__.py.py
+from generate_images_Mapping import *
from fedic import *
from plot_binned_strains_vs_radius import *
from NonlinearSolver import *
@@ -14,7 +15,6 @@ from Energy_Regularization import *
from compute_warped_images import *
from plot_strains import *
from image_expressions_cpp import *
-from resample_images import *
from Energy import *
from compute_warped_mesh import *
from Energy_WarpedImage import *
@@ -32,5 +32,7 @@ from compute_displacements import *
from generate_images import *
from ImageSeries import *
from compute_displacement_error import *
+from downsample_images import *
from project import *
from compute_strains import *
+from generate_images_Image import *
diff --git a/downsample_images.py b/downsample_images.py
new file mode 100644
index 0000000..20f0163
--- /dev/null
+++ b/downsample_images.py
@@ -0,0 +1,275 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Martin Genet, 2016-2019 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+import glob
+import math
+import numpy
+import os
+import random
+import vtk
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+
+################################################################################
+
+def downsample_images(
+ images_folder,
+ images_basename,
+ downsampling_factors,
+ images_ext="vti",
+ keep_resolution=0,
+ write_temp_images=0,
+ verbose=0):
+
+ mypy.my_print(verbose, "*** downsample_images ***")
+
+ images_filenames = glob.glob(images_folder+"/"+images_basename+"_[0-9]*"+"."+images_ext)
+ images_nframes = len(images_filenames)
+ images_zfill = len(images_filenames[0].rsplit("_",1)[-1].split(".",1)[0])
+
+ image = myvtk.readImage(
+ filename=images_folder+"/"+images_basename+"_"+str(0).zfill(images_zfill)+"."+images_ext,
+ verbose=0)
+ images_ndim = myvtk.getImageDimensionality(
+ image=image,
+ verbose=0)
+ images_spacing = image.GetSpacing()
+ images_origin = image.GetOrigin()
+ 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())
+
+ 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)
+ 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_npoints = numpy.prod(images_downsampled_nvoxels)
+ # mypy.my_print(verbose, "images_downsampled_npoints = "+str(images_downsampled_npoints))
+
+ if (images_ext == "vtk"):
+ reader = vtk.vtkImageReader()
+ writer = vtk.vtkImageWriter()
+ if (write_temp_images):
+ writer_fft = vtk.vtkImageWriter()
+ if (keep_resolution):
+ writer_mul = vtk.vtkImageWriter()
+ else:
+ writer_sel = vtk.vtkImageWriter()
+ elif (images_ext == "vti"):
+ reader = vtk.vtkXMLImageDataReader()
+ writer = vtk.vtkXMLImageDataWriter()
+ if (write_temp_images):
+ writer_fft = vtk.vtkXMLImageDataWriter()
+ if (keep_resolution):
+ writer_mul = vtk.vtkXMLImageDataWriter()
+ else:
+ writer_sel = vtk.vtkXMLImageDataWriter()
+ else:
+ assert 0, "\"ext\" must be \".vtk\" or \".vti\". Aborting."
+
+ fft = vtk.vtkImageFFT()
+ fft.SetDimensionality(images_ndim)
+ fft.SetInputConnection(reader.GetOutputPort())
+ if (write_temp_images): writer_fft.SetInputConnection(fft.GetOutputPort())
+
+ if (keep_resolution):
+ image_filename = images_folder+"/"+images_basename+"_"+str(0).zfill(images_zfill)+"."+images_ext
+ mask_image = myvtk.readImage(
+ filename=image_filename,
+ verbose=0)
+ mask_scalars = myvtk.createDoubleArray(
+ name="ImageScalars",
+ n_components=2,
+ n_tuples=images_npoints,
+ verbose=0)
+ mask_image.GetPointData().SetScalars(mask_scalars)
+ # print mask_image.GetScalarType()
+ # print mask_image.GetPointData().GetScalars()
+ if (images_ndim == 1):
+ for k_x in xrange(images_nvoxels[0]):
+ if ((k_x > images_downsampled_nvoxels[0]/2) \
+ and (k_x < images_nvoxels[0] - images_downsampled_nvoxels[0]/2)):
+ mask_scalars.SetTuple(k_x, [0, 0])
+ else:
+ mask_scalars.SetTuple(k_x, [1, 1])
+ if (images_ndim == 2):
+ for k_y in xrange(images_nvoxels[1]):
+ for k_x in xrange(images_nvoxels[0]):
+ k_point = k_y*images_nvoxels[0] + k_x
+ if (((k_x > images_downsampled_nvoxels[0]/2) \
+ and (k_x < images_nvoxels[0] - images_downsampled_nvoxels[0]/2)) \
+ or ((k_y > images_downsampled_nvoxels[1]/2) \
+ and (k_y < images_nvoxels[1] - images_downsampled_nvoxels[1]/2))):
+ mask_scalars.SetTuple(k_point, [0, 0])
+ else:
+ mask_scalars.SetTuple(k_point, [1, 1])
+ if (images_ndim == 3):
+ for k_z in xrange(images_nvoxels[2]):
+ for k_y in xrange(images_nvoxels[1]):
+ for k_x in xrange(images_nvoxels[0]):
+ k_point = k_z*images_nvoxels[1]*images_nvoxels[0] + k_y*images_nvoxels[0] + k_x
+ if (((k_x > images_downsampled_nvoxels[0]/2) \
+ and (k_x < images_nvoxels[0] - images_downsampled_nvoxels[0]/2)) \
+ or ((k_y > images_downsampled_nvoxels[1]/2) \
+ and (k_y < images_nvoxels[1] - images_downsampled_nvoxels[1]/2)) \
+ or ((k_z > images_downsampled_nvoxels[2]/2) \
+ and (k_z < images_nvoxels[2] - images_downsampled_nvoxels[2]/2))):
+ mask_scalars.SetTuple(k_point, [0, 0])
+ else:
+ mask_scalars.SetTuple(k_point, [1, 1])
+ if (write_temp_images): myvtk.writeImage(
+ image=mask_image,
+ filename=images_folder+"/"+images_basename+"_mask"+"."+images_ext,
+ verbose=0)
+
+ mult = vtk.vtkImageMathematics()
+ mult.SetOperationToMultiply()
+ mult.SetInputConnection(0, fft.GetOutputPort())
+ mult.SetInputData(1, mask_image)
+ if (write_temp_images): writer_mul.SetInputConnection(mult.GetOutputPort())
+ 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))
+
+ image_downsampled_scalars = myvtk.createDoubleArray(
+ name="ImageScalars",
+ n_components=2,
+ n_tuples=images_downsampled_npoints,
+ verbose=0)
+ image_downsampled.GetPointData().SetScalars(image_downsampled_scalars)
+ I = numpy.empty(2)
+
+ if (write_temp_images):
+ writer_sel.SetInputData(image_downsampled)
+
+ rfft = vtk.vtkImageRFFT()
+ rfft.SetDimensionality(images_ndim)
+ if (keep_resolution):
+ rfft.SetInputConnection(mult.GetOutputPort())
+ else:
+ rfft.SetInputData(image_downsampled) # MG20190520: Not sure why this does not work.
+
+ writer.SetInputConnection(rfft.GetOutputPort())
+
+ if (keep_resolution):
+ for k_frame in xrange(images_nframes):
+ mypy.my_print(verbose, "k_frame = "+str(k_frame))
+
+ reader.SetFileName(images_folder+"/"+images_basename+"_"+str(k_frame).zfill(images_zfill)+"."+images_ext)
+
+ if (write_temp_images):
+ writer_fft.SetFileName(images_folder+"/"+images_basename+"_fft"+"_"+str(k_frame).zfill(images_zfill)+"."+images_ext)
+ writer_fft.Write()
+
+ if (write_temp_images):
+ writer_mul.SetFileName(images_folder+"/"+images_basename+"_mul"+"_"+str(k_frame).zfill(images_zfill)+"."+images_ext)
+ writer_mul.Write()
+
+ writer.SetFileName(images_folder+"/"+images_basename+"_downsampled"+"_"+str(k_frame).zfill(images_zfill)+"."+images_ext)
+ writer.Write()
+ else:
+ for k_frame in xrange(images_nframes):
+ mypy.my_print(verbose, "k_frame = "+str(k_frame))
+
+ reader.SetFileName(images_folder+"/"+images_basename+"_"+str(k_frame).zfill(images_zfill)+"."+images_ext)
+ reader.Update()
+ # print "reader.GetOutput() = "+str(reader.GetOutput())
+
+ fft.Update()
+ if (write_temp_images):
+ writer_fft.SetFileName(images_folder+"/"+images_basename+"_fft"+"_"+str(k_frame).zfill(images_zfill)+"."+images_ext)
+ writer_fft.Write()
+ # print "fft.GetOutput() = "+str(fft.GetOutput())
+ # print "fft.GetOutput().GetScalarType() = "+str(fft.GetOutput().GetScalarType())
+ # print "fft.GetOutput().GetPointData().GetScalars() = "+str(fft.GetOutput().GetPointData().GetScalars())
+
+ image_scalars = fft.GetOutput().GetPointData().GetScalars()
+ if (images_ndim == 1):
+ for k_x_downsampled in xrange(images_downsampled_nvoxels[0]):
+ k_x = k_x_downsampled if (k_x_downsampled <= images_downsampled_nvoxels[0]/2) else k_x_downsampled+(images_nvoxels[0]-images_downsampled_nvoxels[0])
+ image_scalars.GetTuple(k_x, I)
+ I /= downsampling_factor
+ image_downsampled_scalars.SetTuple(k_x_downsampled, I)
+ if (images_ndim == 2):
+ for k_y_downsampled in xrange(images_downsampled_nvoxels[1]):
+ k_y = k_y_downsampled if (k_y_downsampled <= images_downsampled_nvoxels[1]/2) else k_y_downsampled+(images_nvoxels[1]-images_downsampled_nvoxels[1])
+ for k_x_downsampled in xrange(images_downsampled_nvoxels[0]):
+ k_x = k_x_downsampled if (k_x_downsampled <= images_downsampled_nvoxels[0]/2) else k_x_downsampled+(images_nvoxels[0]-images_downsampled_nvoxels[0])
+ k_point_downsampled = k_y_downsampled*images_downsampled_nvoxels[0] + k_x_downsampled
+ k_point = k_y *images_nvoxels[0] + k_x
+ image_scalars.GetTuple(k_point, I)
+ I /= downsampling_factor
+ image_downsampled_scalars.SetTuple(k_point_downsampled, I)
+ if (images_ndim == 3):
+ for k_z_downsampled in xrange(images_downsampled_nvoxels[2]):
+ k_z = k_z_downsampled if (k_z_downsampled <= images_downsampled_nvoxels[2]/2) else k_z_downsampled+(images_nvoxels[2]-images_downsampled_nvoxels[2])
+ for k_y_downsampled in xrange(images_downsampled_nvoxels[1]):
+ k_y = k_y_downsampled if (k_y_downsampled <= images_downsampled_nvoxels[1]/2) else k_y_downsampled+(images_nvoxels[1]-images_downsampled_nvoxels[1])
+ for k_x_downsampled in xrange(images_downsampled_nvoxels[0]):
+ k_x = k_x_downsampled if (k_x_downsampled <= images_downsampled_nvoxels[0]/2) else k_x_downsampled+(images_nvoxels[0]-images_downsampled_nvoxels[0])
+ k_point_downsampled = k_z_downsampled*images_downsampled_nvoxels[1]*images_downsampled_nvoxels[0] + k_y_downsampled*images_downsampled_nvoxels[0] + k_x_downsampled
+ k_point = k_z *images_nvoxels[1] *images_nvoxels[0] + k_y *images_nvoxels[0] + k_x
+ image_scalars.GetTuple(k_point, I)
+ I /= downsampling_factor
+ image_downsampled_scalars.SetTuple(k_point_downsampled, I)
+ # print "image_downsampled = "+str(image_downsampled)
+ # print "image_downsampled_scalars = "+str(image_downsampled_scalars)
+
+ if (write_temp_images):
+ # writer_sel.SetInputData(image_downsampled)
+ writer_sel.SetFileName(images_folder+"/"+images_basename+"_sel"+"_"+str(k_frame).zfill(images_zfill)+"."+images_ext)
+ writer_sel.Write()
+
+ rfft = vtk.vtkImageRFFT() # MG20190520: Not sure why this is needed.
+ rfft.SetDimensionality(images_ndim) # MG20190520: Not sure why this is needed.
+ rfft.SetInputData(image_downsampled) # MG20190520: Not sure why this is needed.
+ rfft.Update()
+
+ writer.SetInputData(rfft.GetOutput()) # MG20190520: Not sure why this is needed.
+ writer.SetFileName(images_folder+"/"+images_basename+"_downsampled"+"_"+str(k_frame).zfill(images_zfill)+"."+images_ext)
+ writer.Write()
diff --git a/generate_images.py b/generate_images.py
index 63ebfa1..cde5c46 100644
--- a/generate_images.py
+++ b/generate_images.py
@@ -52,7 +52,7 @@ def set_I_wGrad(
G = numpy.dot(G, Finv)
vtk_gradient_vectors.SetTuple(k_point, G)
-def generateImages(
+def generate_images(
images,
structure,
texture,
@@ -62,13 +62,13 @@ def generateImages(
generate_image_gradient=False,
verbose=0):
- mypy.my_print(verbose, "*** generateImages ***")
+ mypy.my_print(verbose, "*** generate_images ***")
assert ("n_integration" not in images),\
- "\"n_integration\" has been deprecated. Use \"resampling\" instead. Aborting."
+ "\"n_integration\" has been deprecated. Use \"upsampling\" instead. Aborting."
- if ("resampling_factors" not in images):
- images["resampling_factors"] = [1]*images["n_dim"]
+ if ("upsampling_factors" not in images):
+ images["upsampling_factors"] = [1]*images["n_dim"]
if ("zfill" not in images):
images["zfill"] = len(str(images["n_frames"]))
if ("ext" not in images):
@@ -92,30 +92,37 @@ def generateImages(
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):
- vtk_image.SetExtent([0, images["n_voxels"][0]*images["resampling_factors"][0]-1, 0, 0, 0, 0])
+ extent = [0, images["n_voxels_upsampled"][0]-1, 0, 0, 0, 0]
elif (images["n_dim"] == 2):
- vtk_image.SetExtent([0, images["n_voxels"][0]*images["resampling_factors"][0]-1, 0, images["n_voxels"][1]*images["resampling_factors"][1]-1, 0, 0])
+ extent = [0, images["n_voxels_upsampled"][0]-1, 0, images["n_voxels_upsampled"][1]-1, 0, 0]
elif (images["n_dim"] == 3):
- vtk_image.SetExtent([0, images["n_voxels"][0]*images["resampling_factors"][0]-1, 0, images["n_voxels"][1]*images["resampling_factors"][1]-1, 0, images["n_voxels"][2]*images["resampling_factors"][2]-1])
- else:
- assert (0), "n_dim must be \"1\", \"2\" or \"3\". Aborting."
+ 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 = numpy.array([images["L"][0]/images["n_voxels"][0]/images["resampling_factors"][0], 1., 1.])
+ spacing = [delta_upsampled[0], 1., 1.]
elif (images["n_dim"] == 2):
- spacing = numpy.array([images["L"][0]/images["n_voxels"][0]/images["resampling_factors"][0], images["L"][1]/images["n_voxels"][1]/images["resampling_factors"][1], 1.])
+ spacing = [delta_upsampled[0], delta_upsampled[1], 1.]
elif (images["n_dim"] == 3):
- spacing = numpy.array([images["L"][0]/images["n_voxels"][0]/images["resampling_factors"][0], images["L"][1]/images["n_voxels"][1]/images["resampling_factors"][1], images["L"][2]/images["n_voxels"][2]/images["resampling_factors"][2]])
+ 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 = numpy.array([spacing[0]/2, 0., 0.])
+ origin = [delta[0]/2, 0., 0.]
elif (images["n_dim"] == 2):
- origin = numpy.array([spacing[0]/2, spacing[1]/2, 0.])
+ origin = [delta[0]/2, delta[1]/2, 0.]
elif (images["n_dim"] == 3):
- origin = numpy.array([spacing[0]/2, spacing[1]/2, spacing[2]/2])
+ 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(
@@ -155,7 +162,7 @@ def generateImages(
for k_frame in xrange(images["n_frames"]):
t = images["T"]*float(k_frame)/(images["n_frames"]-1) if (images["n_frames"]>1) else 0.
- print "t = "+str(t)
+ mypy.my_print(verbose, "t = "+str(t))
mapping.init_t(t)
for k_point in xrange(n_points):
vtk_image.GetPoint(k_point, x)
@@ -176,23 +183,26 @@ def generateImages(
image=vtk_gradient,
filename=images["folder"]+"/"+images["basename"]+"-grad"+"_"+str(k_frame).zfill(images["zfill"])+"."+images["ext"],
verbose=verbose-1)
- print "global_min = "+str(global_min)
- print "global_max = "+str(global_max)
+ # mypy.my_print(verbose, "global_min = "+str(global_min))
+ # mypy.my_print(verbose, "global_max = "+str(global_max))
- if (images["resampling_factors"] == [1]*images["n_dim"]):
- resampling = False
+ if (images["upsampling_factors"] == [1]*images["n_dim"]):
+ downsampling = False
else:
- resampling = True
- ddic.resample_images(
+ downsampling = True
+ ddic.downsample_images(
images_folder=images["folder"],
images_basename=images["basename"],
- resampling_factors=images["resampling_factors"],
- write_temp_images=1)
+ downsampling_factors=images["upsampling_factors"],
+ keep_resolution=0,
+ write_temp_images=0,
+ verbose=verbose)
if (images["data_type"] in ("float")):
pass
elif (images["data_type"] in ("unsigned char", "unsigned short", "unsigned int", "unsigned long", "unsigned float", "uint8", "uint16", "uint32", "uint64", "ufloat")):
ddic.normalize_images(
images_folder=images["folder"],
- images_basename=images["basename"]+("_resampled"*(resampling)),
- images_datatype=images["data_type"])
+ images_basename=images["basename"]+("_downsampled"*(downsampling)),
+ images_datatype=images["data_type"],
+ verbose=verbose)
diff --git a/generate_undersampled_images.py b/generate_undersampled_images.py
index 50418df..5ba1a0d 100644
--- a/generate_undersampled_images.py
+++ b/generate_undersampled_images.py
@@ -39,7 +39,7 @@ def generateUndersampledImages(
images["n_voxels"][2] /= undersampling_level
images["basename"] = images_basename+"-X"
texture["type"] = "taggX"
- fedic.generateImages(
+ ddic.generate_images(
images=images,
structure=structure,
texture=texture,
@@ -55,7 +55,7 @@ def generateUndersampledImages(
images["n_voxels"][2] /= undersampling_level
images["basename"] = images_basename+"-Y"
texture["type"] = "taggY"
- fedic.generateImages(
+ ddic.generate_images(
images=images,
structure=structure,
texture=texture,
@@ -71,7 +71,7 @@ def generateUndersampledImages(
images["n_voxels"][1] /= undersampling_level
images["basename"] = images_basename+"-Z"
texture["type"] = "taggZ"
- fedic.generateImages(
+ ddic.generate_images(
images=images,
structure=structure,
texture=texture,
diff --git a/normalize_images.py b/normalize_images.py
index edc7296..5de154c 100644
--- a/normalize_images.py
+++ b/normalize_images.py
@@ -26,7 +26,10 @@ def normalize_images(
images_folder,
images_basename,
images_datatype,
- images_ext="vti"):
+ images_ext="vti",
+ verbose=0):
+
+ mypy.my_print(verbose, "*** normalize_images ***")
images_filenames = glob.glob(images_folder+"/"+images_basename+"_[0-9]*"+"."+images_ext)
images_nframes = len(images_filenames)
@@ -49,19 +52,19 @@ def normalize_images(
global_min = float("+Inf")
global_max = float("-Inf")
- I = numpy.empty(2)
for k_frame in xrange(images_nframes):
reader.SetFileName(images_folder+"/"+images_basename+"_"+str(k_frame).zfill(images_zfill)+"."+images_ext)
reader.Update()
image_scalars = reader.GetOutput().GetPointData().GetScalars()
+ I = numpy.empty(image_scalars.GetNumberOfComponents())
for k_point in xrange(images_npoints):
image_scalars.GetTuple(k_point, I)
global_min = min(global_min, I[0])
global_max = max(global_max, I[0])
- print "global_min = "+str(global_min)
- print "global_max = "+str(global_max)
+ mypy.my_print(verbose, "global_min = "+str(global_min))
+ mypy.my_print(verbose, "global_max = "+str(global_max))
shifter = vtk.vtkImageShiftScale()
shifter.SetInputConnection(reader.GetOutputPort())
@@ -85,11 +88,8 @@ def normalize_images(
writer.SetInputConnection(shifter.GetOutputPort())
for k_frame in xrange(images_nframes):
- reader.SetFileName(images_folder+"/"+images_basename+"_"+str(k_frame).zfill(images_zfill)+"."+images_ext)
- reader.Update()
-
- shifter.Update()
+ mypy.my_print(verbose, "k_frame = "+str(k_frame))
+ reader.SetFileName(images_folder+"/"+images_basename +"_"+str(k_frame).zfill(images_zfill)+"."+images_ext)
writer.SetFileName(images_folder+"/"+images_basename+"_normalized"+"_"+str(k_frame).zfill(images_zfill)+"."+images_ext)
- writer.Update()
writer.Write()
diff --git a/resample_images.py b/resample_images.py
deleted file mode 100644
index 6488ec7..0000000
--- a/resample_images.py
+++ /dev/null
@@ -1,151 +0,0 @@
-#coding=utf8
-
-################################################################################
-### ###
-### Created by Martin Genet, 2016-2019 ###
-### ###
-### École Polytechnique, Palaiseau, France ###
-### ###
-################################################################################
-
-import glob
-import math
-import numpy
-import os
-import random
-import vtk
-
-import myPythonLibrary as mypy
-import myVTKPythonLibrary as myvtk
-
-import dolfin_dic as ddic
-
-################################################################################
-
-def resample_images(
- images_folder,
- images_basename,
- resampling_factors,
- images_ext="vti",
- write_temp_images=False):
-
- images_filenames = glob.glob(images_folder+"/"+images_basename+"_[0-9]*"+"."+images_ext)
- images_nframes = len(images_filenames)
- images_zfill = len(images_filenames[0].rsplit("_",1)[-1].split(".",1)[0])
-
- image_filename = images_folder+"/"+images_basename+"_"+str(0).zfill(images_zfill)+"."+images_ext
- image = myvtk.readImage(
- filename=image_filename,
- verbose=0)
- images_npoints = image.GetNumberOfPoints()
- images_nvoxels = myvtk.getImageDimensions(
- image=image,
- verbose=0)
- images_ndim = len(images_nvoxels)
-
- if (images_ext == "vtk"):
- reader = vtk.vtkImageReader()
- writer = vtk.vtkImageWriter()
- if (write_temp_images):
- writer_fft = vtk.vtkImageWriter()
- writer_mult = vtk.vtkImageWriter()
- elif (images_ext == "vti"):
- reader = vtk.vtkXMLImageDataReader()
- writer = vtk.vtkXMLImageDataWriter()
- if (write_temp_images):
- writer_fft = vtk.vtkXMLImageDataWriter()
- writer_mult = vtk.vtkXMLImageDataWriter()
- else:
- assert 0, "\"ext\" must be \".vtk\" or \".vti\". Aborting."
-
- fft = vtk.vtkImageFFT()
- fft.SetDimensionality(images_ndim)
- fft.SetInputConnection(reader.GetOutputPort())
- if (write_temp_images): writer_fft.SetInputConnection(fft.GetOutputPort())
-
- image_filename = images_folder+"/"+images_basename+"_"+str(0).zfill(images_zfill)+"."+images_ext
- mask_image = myvtk.readImage(
- filename=image_filename,
- verbose=0)
- mask_scalars = myvtk.createDoubleArray(
- name="ImageScalars",
- n_components=2,
- n_tuples=images_npoints,
- verbose=0)
- mask_image.GetPointData().SetScalars(mask_scalars)
- # print mask_image.GetScalarType()
- # print mask_image.GetPointData().GetScalars()
- if (images_ndim == 1):
- for k_point in xrange(images_npoints):
- if ((k_point > images_nvoxels[0]/resampling_factors[0]/2) \
- and (k_point < images_nvoxels[0] - images_nvoxels[0]/resampling_factors[0]/2)):
- mask_scalars.SetTuple(k_point, [0, 0])
- else:
- mask_scalars.SetTuple(k_point, [1, 1])
- if (images_ndim == 2):
- for k_point in xrange(images_npoints):
- # print "k_point = "+str(k_point)
- k_y = math.floor(k_point/images_nvoxels[0])
- k_x = k_point - k_y*images_nvoxels[0]
- # print "k_x = "+str(k_x)
- # print "k_y = "+str(k_y)
- if (((k_x > images_nvoxels[0]/resampling_factors[0]/2) \
- and (k_x < images_nvoxels[0] - images_nvoxels[0]/resampling_factors[0]/2)) \
- or ((k_y > images_nvoxels[1]/resampling_factors[1]/2) \
- and (k_y < images_nvoxels[1] - images_nvoxels[1]/resampling_factors[1]/2))):
- mask_scalars.SetTuple(k_point, [0, 0])
- else:
- mask_scalars.SetTuple(k_point, [1, 1])
- if (images_ndim == 3):
- for k_point in xrange(images_npoints):
- k_z = math.floor(k_point/images_nvoxels[0]/images_nvoxels[1])
- k_y = math.floor((k_point - k_z*images_nvoxels[0]*images_nvoxels[1])/images_nvoxels[0])
- k_x = k_point - k_z*images_nvoxels[0]*images_nvoxels[1] - k_y*images_nvoxels[0]
- if (((k_x > images_nvoxels[0]/resampling_factors[0]/2) \
- and (k_x < images_nvoxels[0] - images_nvoxels[0]/resampling_factors[0]/2)) \
- or ((k_y > images_nvoxels[1]/resampling_factors[1]/2) \
- and (k_y < images_nvoxels[1] - images_nvoxels[1]/resampling_factors[1]/2)) \
- or ((k_z > images_nvoxels[2]/resampling_factors[2]/2) \
- and (k_z < images_nvoxels[2] - images_nvoxels[2]/resampling_factors[2]/2))):
- mask_scalars.SetTuple(k_point, [0, 0])
- else:
- mask_scalars.SetTuple(k_point, [1, 1])
- if (write_temp_images): myvtk.writeImage(
- image=mask_image,
- filename=images_folder+"/"+images_basename+"_mask"+"."+images_ext,
- verbose=0)
-
- mult = vtk.vtkImageMathematics()
- mult.SetOperationToMultiply()
- mult.SetInputData(0, mask_image)
- mult.SetInputConnection(1, fft.GetOutputPort())
- if (write_temp_images): writer_mult.SetInputConnection(mult.GetOutputPort())
-
- rfft = vtk.vtkImageRFFT()
- rfft.SetDimensionality(images_ndim)
- rfft.SetInputConnection(mult.GetOutputPort())
-
- writer.SetInputConnection(rfft.GetOutputPort())
-
- for k_frame in xrange(images_nframes):
- reader.SetFileName(images_folder+"/"+images_basename+"_"+str(k_frame).zfill(images_zfill)+"."+images_ext)
- reader.Update()
-
- fft.Update()
- if (write_temp_images):
- writer_fft.SetFileName(images_folder+"/"+images_basename+"_fft"+"_"+str(k_frame).zfill(images_zfill)+"."+images_ext)
- writer_fft.Update()
- writer_fft.Write()
- # print fft.GetOutput().GetScalarType()
- # print fft.GetOutput().GetPointData().GetScalars()
-
- mult.Update()
- if (write_temp_images):
- writer_mult.SetFileName(images_folder+"/"+images_basename+"_mult"+"_"+str(k_frame).zfill(images_zfill)+"."+images_ext)
- writer_mult.Update()
- writer_mult.Write()
-
- rfft.Update()
- writer.SetFileName(images_folder+"/"+images_basename+"_resampled"+"_"+str(k_frame).zfill(images_zfill)+"."+images_ext)
- writer.Update()
- writer.Write()
--
GitLab