generateImages.py

#coding=utf8

########################################################################
###                                                                  ###
### Created by Martin Genet, 2012-2016                               ###
###                                                                  ###
### University of California at San Francisco (UCSF), USA            ###
### Swiss Federal Institute of Technology (ETH), Zurich, Switzerland ###
### École Polytechnique, Palaiseau, France                           ###
###                                                                  ###
########################################################################

import glob
import math
import numpy
import os
import random
import vtk

import myVTKPythonLibrary as myVTK

########################################################################

#class ImagesInfo():
    #def __init__(self, n_dim, L, n_voxels, n_integration, T, n_frames, data_type, images_folder, images_basename):
        #assert (n_dim in (1,2,3))
        #self.n_dim = n_dim

        #if (type(L) == float):
            #assert (L>0)
            #self.L = numpy.array([L]*self.n_dim)
        #elif (type(L) == int):
            #assert (L>0)
            #self.L = numpy.array([float(L)]*self.n_dim)
        #else:
            #assert (len(L) == self.n_dim)
            #self.L = numpy.array(L)
            #assert ((self.L>0).all())

        #if (type(n_voxels) == int):
            #assert (n_voxels>0)
            #self.n_voxels = numpy.array([n_voxels]*self.n_dim)
        #else:
            #assert (len(n_voxels) == self.n_dim)
            #self.n_voxels = numpy.array(n_voxels)
            #assert ((self.n_voxels>0).all())

        #if (type(n_integration) == int):
            #assert (n_integration>0)
            #self.n_integration = numpy.array([n_integration]*self.n_dim)
        #else:
            #assert (len(n_integration) == self.n_dim)
            #self.n_integration = numpy.array(n_integration)
            #assert ((self.n_integration>0).all())

        #assert (T>0.)
        #self.T = T

        #assert (n_frames>0)
        #self.n_frames = n_frames

        #assert (data_type in ("int", "float", "unsigned char", "unsigned short", "unsigned int", "unsigned long", "unsigned float" "uint8", "uint16", "uint32", "uint64", "ufloat"))
        #self.data_type = data_type

        #self.images_folder = images_folder
        #self.images_basename = images_basename

#class StructureInfo():
    #def __init__(self, images, type, **kwargs):
        #assert (type in ("no", "heart"))
        #self["type"] = type
        #if (self["type"] == "heart"):
            #self.Ri = kwargs["Ri"]
            #self.Re = kwargs["Re"]
            #if (images.n_dim == 3):
                #self.Zmin = kwargs["Zmin"] if ("Zmin" in kwargs.keys()) else 0.
                #self.Zmax = kwargs["Zmax"] if ("Zmax" in kwargs.keys()) else images.L[2]

#class TextureInfo():
    #def __init__(self, type, **kwargs):
        #assert (type in ("no", "taggX", "taggY", "taggZ"))
        #self["type"] = type

#class NoiseInfo():
    #def __init__(self, type, **kwargs):
        #self["type"] = type

#class DeformationInfo():
    #def __init__(self, type, **kwargs):
        #self["type"] = type

#class EvolutionInfo():
    #def __init__(self, type, **kwargs):
        #self["type"] = type

########################################################################

class Image():
    def __init__(self, images, structure, texture, noise):
        self.L = images["L"]

        # structure
        if (structure["type"] == "no"):
            self.I0_structure = self.I0_structure_no
        elif (structure["type"] == "heart"):
            if (images["n_dim"] == 2):
                self.I0_structure = self.I0_structure_heart_2
                self.R = float()
                self.Ri = structure["Ri"]
                self.Re = structure["Re"]
            elif (images["n_dim"] == 3):
                self.I0_structure = self.I0_structure_heart_3
                self.R = float()
                self.Ri = structure["Ri"]
                self.Re = structure["Re"]
                self.Zmin = structure.Zmin if ("Zmin" in structure.keys()) else 0.
                self.Zmax = structure.Zmax if ("Zmax" in structure.keys()) else images["L"][2]
            else:
                assert (0), "n_dim must be \"2\" or \"3 for \"heart\" type structure. Aborting."
        else:
            assert (0), "structure type must be \"no\" or \"heart\". Aborting."

        # texture
        if (texture["type"] == "no"):
            self.I0_texture = self.I0_texture_no
        elif (texture["type"].startswith("tagging")):
            if   (images["n_dim"] == 1):
                if ("-signed" in texture["type"]):
                    self.I0_texture = self.I0_texture_tagging_signed_X
                else:
                    self.I0_texture = self.I0_texture_tagging_X
            elif (images["n_dim"] == 2):
                if ("-signed" in texture["type"]):
                    self.I0_texture = self.I0_texture_tagging_signed_XY
                else:
                    self.I0_texture = self.I0_texture_tagging_XY
            elif (images["n_dim"] == 3):
                if ("-signed" in texture["type"]):
                    self.I0_texture = self.I0_texture_tagging_signed_XYZ
                else:
                    self.I0_texture = self.I0_texture_tagging_XYZ
            else:
                assert (0), "n_dim must be \"1\", \"2\" or \"3\". Aborting."
            self.s = texture["s"]
        elif (texture["type"].startswith("taggX")):
            self.I0_texture = self.I0_texture_tagging_X
            self.s = texture["s"]
        elif (texture["type"].startswith("taggY")):
            self.I0_texture = self.I0_texture_tagging_Y
            self.s = texture["s"]
        elif (texture["type"].startswith("taggZ")):
            self.I0_texture = self.I0_texture_tagging_Z
            self.s = texture["s"]
        else:
            assert (0), "texture type must be \"no\", \"tagging\", \"taggX\", \"taggY\" or \"taggZ\". Aborting."

        # noise
        if (noise["type"] == "no"):
            self.I0_noise = self.I0_noise_no
        elif (noise["type"] == "normal"):
            self.I0_noise = self.I0_noise_normal
            self.avg = noise.avg if ("avg" in noise.keys()) else 0.
            self.std = noise.std
        else:
            assert (0), "noise type must be \"no\" or \"normal\". Aborting."

    def I0(self, X, i, g=None):
        self.I0_structure(X, i, g)
        self.I0_texture(X, i, g)
        self.I0_noise(i, g)

    def I0_structure_no(self, X, i, g=None):
        i[0] = 1.
        #if (g is not None): g[:] = 1.

    def I0_structure_heart_2(self, X, i, g=None):
        self.R = ((X[0]-self.L[0]/2)**2 + (X[1]-self.L[1]/2)**2)**(1./2)
        if (self.R >= self.Ri) and (self.R <= self.Re):
            i[0] = 1.
            #if (g is not None): g[:] = 1.
        else:
            i[0] = 0.
            #if (g is not None): g[:] = 0.
        #i[0] = 1.

    def I0_structure_heart_3(self, X, i, g=None):
        self.R = ((X[0]-self.L[0]/2)**2 + (X[1]-self.L[1]/2)**2)**(1./2)
        if (self.R >= self.Ri) and (self.R <= self.Re) and (X[2] >= self.Zmin) and (X[2] <= self.Zmax):
            i[0] = 1.
            #if (g is not None): g[:] = 1.
        else:
            i[0] = 0.
            #if (g is not None): g[:] = 0.

    def I0_texture_no(self, X, i, g=None):
        i[0] *= 1.
        #if (g is not None): g[:] *= 0.

    def I0_texture_tagging_signed_X(self, X, i, g=None):
        i[0] *= math.sin(math.pi*X[0]/self.s)
        #if (g is not None):
            #assert (0), "ToDo"

    def I0_texture_tagging_X(self, X, i, g=None):
        i[0] *= abs(math.sin(math.pi*X[0]/self.s))
        #if (g is not None):
            #assert (0), "ToDo"

    def I0_texture_tagging_signed_Y(self, X, i, g=None):
        i[0] *= math.sin(math.pi*X[1]/self.s)
        #if (g is not None):
            #assert (0), "ToDo"

    def I0_texture_tagging_Y(self, X, i, g=None):
        i[0] *= abs(math.sin(math.pi*X[1]/self.s))
        #if (g is not None):
            #assert (0), "ToDo"

    def I0_texture_tagging_signed_Z(self, X, i, g=None):
        i[0] *= math.sin(math.pi*X[2]/self.s)
        #if (g is not None):
            #assert (0), "ToDo"

    def I0_texture_tagging_Z(self, X, i, g=None):
        i[0] *= abs(math.sin(math.pi*X[2]/self.s))
        #if (g is not None):
            #assert (0), "ToDo"

    def I0_texture_tagging_signed_XY(self, X, i, g=None):
        i[0] *= (math.sin(math.pi*X[0]/self.s)
              +  math.sin(math.pi*X[1]/self.s))/2
        assert (i[0] >= 0.)
        #if (g is not None):
            #assert (0), "ToDo"

    def I0_texture_tagging_XY(self, X, i, g=None):
        i[0] *= (abs(math.sin(math.pi*X[0]/self.s))
              +  abs(math.sin(math.pi*X[1]/self.s)))/2
        assert (i[0] >= 0.)
        #if (g is not None):
            #assert (0), "ToDo"

    def I0_texture_tagging_signed_XYZ(self, X, i, g=None):
        i[0] *= (math.sin(math.pi*X[0]/self.s)
              +  math.sin(math.pi*X[1]/self.s)
              +  math.sin(math.pi*X[2]/self.s))/3
        #if (g is not None):
            #assert (0), "ToDo"

    def I0_texture_tagging_XYZ(self, X, i, g=None):
        i[0] *= (abs(math.sin(math.pi*X[0]/self.s))
              +  abs(math.sin(math.pi*X[1]/self.s))
              +  abs(math.sin(math.pi*X[2]/self.s)))/3
        #if (g is not None):
            #assert (0), "ToDo"

    def I0_noise_no(self, i, g=None):
        pass

    def I0_noise_normal(self, i, g=None):
        i[0] += random.normalvariate(self.avg, self.std)
        #if (g is not None): g[k] += [2*random.normalvariate(self.avg, self.std) for k in xrange(len(g))]

########################################################################

class Mapping:
    def __init__(self, images, structure, deformation, evolution):
        self.deformation = deformation
        if (self.deformation["type"] == "no"):
            self.init_t = self.init_t_no
            self.X = self.X_no
            self.x = self.x_no
        elif (self.deformation["type"] == "trans"):
            self.init_t = self.init_t_trans
            self.X = self.X_trans
            self.x = self.x_trans
            self.D = numpy.empty(3)
        elif (self.deformation["type"] == "rot"):
            self.init_t = self.init_t_rot
            self.X = self.X_rot
            self.x = self.x_rot
            self.C = numpy.empty(3)
            self.R = numpy.empty((3,3))
            self.Rinv = numpy.empty((3,3))
        elif (self.deformation["type"] == "homogeneous"):
            self.init_t = self.init_t_homogeneous
            self.X = self.X_homogeneous
            self.x = self.x_homogeneous
        elif (self.deformation["type"] == "heart"):
            assert (structure["type"] == "heart"), "structure type must be \"heart\" for \"heart\" type deformation. Aborting."
            self.init_t = self.init_t_heart
            self.X = self.X_heart
            self.x = self.x_heart
            self.x_inplane = numpy.empty(2)
            self.X_inplane = numpy.empty(2)
            self.rt = numpy.empty(2)
            self.RT = numpy.empty(2)
            self.L = images["L"]
            self.Ri = structure["Ri"]
            self.Re = structure["Re"]
            self.R = numpy.empty((3,3))
        else:
            assert (0), "deformation type must be \"no\", \"trans\", \"rot\", \"homogeneous\" or \"heart\". Aborting."

        if (evolution["type"] == "linear"):
            self.phi = self.phi_linear
        elif (evolution["type"] == "sine"):
            self.phi = self.phi_sine
            self.T = evolution["T"]
        else:
            assert (0), "evolution ("+evolution["type"]+") type must be \"linear\" or \"sine\". Aborting."

    def phi_linear(self, t):
        return t

    def phi_sine(self, t):
        return math.sin(math.pi*t/self.T)**2

    def init_t_no(self, t):
        pass

    def init_t_trans(self, t):
        self.D[0] = self.deformation["Dx"]*self.phi(t) if ("Dx" in self.deformation.keys()) else 0.
        self.D[1] = self.deformation["Dy"]*self.phi(t) if ("Dy" in self.deformation.keys()) else 0.
        self.D[2] = self.deformation["Dz"]*self.phi(t) if ("Dz" in self.deformation.keys()) else 0.

    def init_t_rot(self, t):
        self.C[0] = self.deformation["Cx"] if ("Cx" in self.deformation.keys()) else 0.
        self.C[1] = self.deformation["Cy"] if ("Cy" in self.deformation.keys()) else 0.
        self.C[2] = self.deformation["Cz"] if ("Cz" in self.deformation.keys()) else 0.
        Rx = self.deformation["Rx"]*math.pi/180*self.phi(t) if ("Rx" in self.deformation.keys()) else 0.
        Ry = self.deformation["Ry"]*math.pi/180*self.phi(t) if ("Ry" in self.deformation.keys()) else 0.
        Rz = self.deformation["Rz"]*math.pi/180*self.phi(t) if ("Rz" in self.deformation.keys()) else 0.
        RRx = numpy.array([[          1. ,           0. ,           0. ],
                           [          0. , +math.cos(Rx), -math.sin(Rx)],
                           [          0. , +math.sin(Rx), +math.cos(Rx)]])
        RRy = numpy.array([[+math.cos(Ry),           0. , +math.sin(Ry)],
                           [          0. ,           1. ,           0. ],
                           [-math.sin(Ry),           0. , +math.cos(Ry)]])
        RRz = numpy.array([[+math.cos(Rz), -math.sin(Rz),           0. ],
                           [+math.sin(Rz), +math.cos(Rz),           0. ],
                           [          0. ,           0. ,           1. ]])
        self.R[:,:] = numpy.dot(numpy.dot(RRx, RRy), RRz)
        self.Rinv[:,:] = numpy.linalg.inv(self.R)

    def init_t_homogeneous(self, t):
        Exx = self.deformation["Exx"]*self.phi(t) if ("Exx" in self.deformation.keys()) else 0.
        Eyy = self.deformation["Eyy"]*self.phi(t) if ("Eyy" in self.deformation.keys()) else 0.
        Ezz = self.deformation["Ezz"]*self.phi(t) if ("Ezz" in self.deformation.keys()) else 0.
        Exy = self.deformation["Exy"]*self.phi(t) if ("Exy" in self.deformation.keys()) else 0.
        Eyx = self.deformation["Eyx"]*self.phi(t) if ("Eyx" in self.deformation.keys()) else 0.
        Exz = self.deformation["Exz"]*self.phi(t) if ("Exz" in self.deformation.keys()) else 0.
        Ezx = self.deformation["Ezx"]*self.phi(t) if ("Ezx" in self.deformation.keys()) else 0.
        Eyz = self.deformation["Eyz"]*self.phi(t) if ("Eyz" in self.deformation.keys()) else 0.
        Ezy = self.deformation["Ezy"]*self.phi(t) if ("Ezy" in self.deformation.keys()) else 0.
        self.F = numpy.array([[math.sqrt(1.+Exx),              Exy ,              Exz ],
                              [             Eyx , math.sqrt(1.+Eyy),              Eyz ],
                              [             Ezx ,              Ezy , math.sqrt(1.+Ezz)]])
        self.Finv = numpy.linalg.inv(self.F)

    def init_t_heart(self, t):
        self.dRi = self.deformation["dRi"]*self.phi(t) if ("dRi" in self.deformation.keys()) else 0.
        self.dRe = self.deformation["dRi"]*self.phi(t) if ("dRi" in self.deformation.keys()) else 0.
        self.dTi = self.deformation["dTi"]*self.phi(t) if ("dTi" in self.deformation.keys()) else 0.
        self.dTe = self.deformation["dTe"]*self.phi(t) if ("dTe" in self.deformation.keys()) else 0.
        self.A = numpy.array([[1.-(self.dRi-self.dRe)/(self.Re-self.Ri), 0.],
                              [  -(self.dTi-self.dTe)/(self.Re-self.Ri), 1.]])
        self.Ainv = numpy.linalg.inv(self.A)
        self.B = numpy.array([(1.+self.Ri/(self.Re-self.Ri))*self.dRi-self.Ri/(self.Re-self.Ri)*self.dRe,
                              (1.+self.Ri/(self.Re-self.Ri))*self.dTi-self.Ri/(self.Re-self.Ri)*self.dTe])

    def X_no(self, x, X, Finv=None):
        X[:] = x
        #if (Finv is not None): Finv[:,:] = numpy.identity(numpy.sqrt(numpy.size(Finv)))

    def X_trans(self, x, X, Finv=None):
        X[:] = x - self.D
        #if (Finv is not None): Finv[:,:] = numpy.identity(numpy.sqrt(numpy.size(Finv)))

    def X_rot(self, x, X, Finv=None):
        X[:] = numpy.dot(self.Rinv, x - self.C) + self.C
        #if (Finv is not None): Finv[:,:] = self.Rinv

    def X_homogeneous(self, x, X, Finv=None):
        X[:] = numpy.dot(self.Finv, x)
        #if (Finv is not None): Finv[:,:] = self.Finv

    def X_heart(self, x, X, Finv=None):
        #print "x = "+str(x)
        self.x_inplane[0] = x[0] - self.L[0]/2
        self.x_inplane[1] = x[1] - self.L[1]/2
        #print "x_inplane = "+str(self.x_inplane)
        self.rt[0] = numpy.linalg.norm(self.x_inplane)
        self.rt[1] = math.atan2(self.x_inplane[1], self.x_inplane[0])
        #print "rt = "+str(self.rt)
        self.RT[:] = numpy.dot(self.Ainv, self.rt-self.B)
        #print "RT = "+str(self.RT)
        X[0] = self.RT[0] * math.cos(self.RT[1]) + self.L[0]/2
        X[1] = self.RT[0] * math.sin(self.RT[1]) + self.L[1]/2
        X[2] = x[2]
        #print "X = "+str(X)
        #if (Finv is not None):
            #Finv[0,0] = 1.+(self.dRe-self.dRi)/(self.Re-self.Ri)
            #Finv[0,1] = 0.
            #Finv[0,2] = 0.
            #Finv[1,0] = (self.dTe-self.dTi)/(self.Re-self.Ri)*self.rt[0]
            #Finv[1,1] = self.rt[0]/self.RT[0]
            #Finv[1,2] = 0.
            #Finv[2,0] = 0.
            #Finv[2,1] = 0.
            #Finv[2,2] = 1.
            ##print "F = "+str(Finv)
            #Finv[:,:] = numpy.linalg.inv(Finv)
            ##print "Finv = "+str(Finv)
            #self.R[0,0] = +math.cos(self.RT[1])
            #self.R[0,1] = +math.sin(self.RT[1])
            #self.R[0,2] = 0.
            #self.R[1,0] = -math.sin(self.RT[1])
            #self.R[1,1] = +math.cos(self.RT[1])
            #self.R[1,2] = 0.
            #self.R[2,0] = 0.
            #self.R[2,1] = 0.
            #self.R[2,2] = 1.
            ##print "R = "+str(self.R)
            #Finv[:] = numpy.dot(numpy.transpose(self.R), numpy.dot(Finv, self.R))
            ##print "Finv = "+str(Finv)

    def x_no(self, X, x, F=None):
        x[:] = X
        #if (F is not None): F[:,:] = numpy.identity(numpy.sqrt(numpy.size(F)))

    def x_trans(self, X, x, F=None):
        x[:] = X + self.D
        #if (F is not None): F[:,:] = numpy.identity(numpy.sqrt(numpy.size(F)))

    def x_rot(self, X, x, F=None):
        x[:] = numpy.dot(self.R, X - self.C) + self.C
        #if (F is not None): F[:,:] = self.R

    def x_homogeneous(self, X, x, F=None):
        x[:] = numpy.dot(self.F, X)
        #if (F is not None): F[:,:] = self.F

    def x_heart(self, X, x, F=None):
        #print "X = "+str(X)
        self.X_inplane[0] = X[0] - self.L[0]/2
        self.X_inplane[1] = X[1] - self.L[1]/2
        #print "X_inplane = "+str(self.X_inplane)
        self.RT[0] = numpy.linalg.norm(self.X_inplane)
        self.RT[1] = math.atan2(self.X_inplane[1], self.X_inplane[0])
        #print "RT = "+str(self.RT)
        self.rt[:] = numpy.dot(self.A, self.RT) + self.B
        #print "rt = "+str(self.rt)
        x[0] = self.rt[0] * math.cos(self.rt[1]) + self.L[0]/2
        x[1] = self.rt[0] * math.sin(self.rt[1]) + self.L[1]/2
        x[2] = X[2]
        #print "x = "+str(x)
        #if (F is not None):
            #F[0,0] = 1.+(self.dRe-self.dRi)/(self.Re-self.Ri)
            #F[0,1] = 0.
            #F[0,2] = 0.
            #F[1,0] = (self.dTe-self.dTi)/(self.Re-self.Ri)*self.rt[0]
            #F[1,1] = self.rt[0]/self.RT[0]
            #F[1,2] = 0.
            #F[2,0] = 0.
            #F[2,1] = 0.
            #F[2,2] = 1.
            ##print "F = "+str(F)
            #self.R[0,0] = +math.cos(self.RT[1])
            #self.R[0,1] = +math.sin(self.RT[1])
            #self.R[0,2] = 0.
            #self.R[1,0] = -math.sin(self.RT[1])
            #self.R[1,1] = +math.cos(self.RT[1])
            #self.R[1,2] = 0.
            #self.R[2,0] = 0.
            #self.R[2,1] = 0.
            #self.R[2,2] = 1.
            #F[:] = numpy.dot(numpy.transpose(self.R), numpy.dot(F, self.R))
            ##print "F = "+str(F)

########################################################################

def generateImages(
        images,
        structure,
        texture,
        noise,
        deformation,
        evolution,
        generate_image_gradient=False,
        verbose=0):

    myVTK.myPrint(verbose, "*** generateImages ***")

    vtk_image = vtk.vtkImageData()

    if   (images["n_dim"] == 1):
        vtk_image.SetExtent([0, images["n_voxels"][0]-1, 0,                       0, 0,                       0])
    elif (images["n_dim"] == 2):
        vtk_image.SetExtent([0, images["n_voxels"][0]-1, 0, images["n_voxels"][1]-1, 0,                       0])
    elif (images["n_dim"] == 3):
        vtk_image.SetExtent([0, images["n_voxels"][0]-1, 0, images["n_voxels"][1]-1, 0, images["n_voxels"][2]-1])
    else:
        assert (0), "n_dim must be \"1\", \"2\" or \"3\". Aborting."

    spacing = numpy.array(images["L"])/numpy.array(images["n_voxels"])
    if   (images["n_dim"] == 1):
        spacing = numpy.array([images["L"][0]/images["n_voxels"][0], 1., 1.])
    elif (images["n_dim"] == 2):
        spacing = numpy.array([images["L"][0]/images["n_voxels"][0], images["L"][1]/images["n_voxels"][1], 1.])
    elif (images["n_dim"] == 2):
        spacing = numpy.array([images["L"][0]/images["n_voxels"][0], images["L"][1]/images["n_voxels"][1], images["L"][2]/images["n_voxels"][2]])
    vtk_image.SetSpacing(spacing)

    origin = numpy.array(vtk_image.GetSpacing())/2
    if   (images["n_dim"] == 1):
        origin[1] = 0.
        origin[2] = 0.
    elif (images["n_dim"] == 2):
        origin[2] = 0.
    vtk_image.SetOrigin(origin)

    n_points = vtk_image.GetNumberOfPoints()
    vtk_image_scalars = myVTK.createFloatArray(
        name="ImageScalars",
        n_components=1,
        n_tuples=n_points,
        verbose=verbose-1)
    vtk_image.GetPointData().SetScalars(vtk_image_scalars)
    if (generate_image_gradient):
        vtk_image_gradient = myVTK.createFloatArray(
            name="ImageScalarsGradient",
            n_components=images["n_dim"],
            n_tuples=n_points,
            verbose=verbose-1)
        vtk_image.GetPointData().SetVectors(vtk_image_gradient)

    if not os.path.exists(images["folder"]):
        os.mkdir(images["folder"])

    x0   = numpy.empty(3)
    x    = numpy.empty(3)
    X    = numpy.empty(3)
    if (generate_image_gradient):
        F    = numpy.empty((3,3))
        Finv = numpy.empty((3,3))
    else:
        F    = None
        Finv = None
    dx   = spacing[0:images["n_dim"]]/images["n_integration"][0:images["n_dim"]]
    global_min = float("+Inf")
    global_max = float("-Inf")
    I = numpy.empty(1)
    i = numpy.empty(1)
    if (generate_image_gradient):
        G = numpy.empty(images["n_dim"])
        g = numpy.empty(images["n_dim"])
    else:
        G = None
        g = None
    image = Image(images, structure, texture, noise)
    mapping = Mapping(images, structure, deformation, evolution)
    if ("zfill" not in images.keys()):
        images["zfill"] = len(str(images["n_frames"]))
    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)
        mapping.init_t(t)
        for k_point in xrange(n_points):
            vtk_image.GetPoint(k_point, x0)
            #print "x0 = "+str(x0)
            x[:] = x0[:]
            #print "x = "+str(x)
            I[0] = 0.
            #print "I = "+str(I)
            #if (generate_image_gradient): G[:] = 0.
                #print "G = "+str(G)
            if   (images["n_dim"] == 1):
                for k_x in xrange(images["n_integration"][0]):
                    x[0] = x0[0] - dx[0]/2 + (k_x+1./2)*dx[0]/images["n_integration"][0]
                    mapping.X(x, X, Finv)
                    image.I0(X, i, g)
                    I += i
                    #if (generate_image_gradient): G += numpy.dot(g, Finv)
                I /= images["n_integration"][0]
                #if (generate_image_gradient): G /= images["n_integration"][0]
            elif (images["n_dim"] == 2):
                for k_y in xrange(images["n_integration"][1]):
                    x[1] = x0[1] - dx[1]/2 + (k_y+1./2)*dx[1]/images["n_integration"][1]
                    for k_x in xrange(images["n_integration"][0]):
                        x[0] = x0[0] - dx[0]/2 + (k_x+1./2)*dx[0]/images["n_integration"][0]
                        #print "x = "+str(x)
                        mapping.X(x, X, Finv)
                        #print "X = "+str(X)
                        #print "Finv = "+str(Finv)
                        image.I0(X, i, g)
                        #print "i = "+str(i)
                        #print "g = "+str(g)
                        I += i
                        #if (generate_image_gradient): G += numpy.dot(g, Finv)
                I /= images["n_integration"][1]*images["n_integration"][0]
                #if (generate_image_gradient):G /= images["n_integration"][1]*images["n_integration"][0]
            elif (images["n_dim"] == 3):
                for k_z in xrange(images["n_integration"][2]):
                    x[2] = x0[2] - dx[2]/2 + (k_z+1./2)*dx[2]/images["n_integration"][2]
                    for k_y in xrange(images["n_integration"][1]):
                        x[1] = x0[1] - dx[1]/2 + (k_y+1./2)*dx[1]/images["n_integration"][1]
                        for k_x in xrange(images["n_integration"][0]):
                            x[0] = x0[0] - dx[0]/2 + (k_x+1./2)*dx[0]/images["n_integration"][0]
                            mapping.X(x, X, Finv)
                            image.I0(X, i, g)
                            I += i
                            #if (generate_image_gradient): G += numpy.dot(g, Finv)
                I /= images["n_integration"][2]*images["n_integration"][1]*images["n_integration"][0]
                #if (generate_image_gradient): G /= images["n_integration"][2]*images["n_integration"][1]*images["n_integration"][0]
            else:
                assert (0), "n_dim must be \"1\", \"2\" or \"3\". Aborting."
            vtk_image_scalars.SetTuple(k_point, I)
            #if (generate_image_gradient): vtk_image_gradient.SetTuple(k_point, G)
            if (I[0] < global_min): global_min = I[0]
            if (I[0] > global_max): global_max = I[0]
        myVTK.writeImage(
            image=vtk_image,
            filename=images["folder"]+"/"+images["basename"]+"_"+str(k_frame).zfill(images["zfill"])+".vti",
            verbose=verbose-1)

    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")):
        #print "global_min = "+str(global_min)
        #print "global_max = "+str(global_max)
        shifter = vtk.vtkImageShiftScale()
        shifter.SetShift(-global_min)
        if   (images["data_type"] in ("unsigned char", "uint8")):
            shifter.SetScale(float(2**8-1)/(global_max-global_min))
            shifter.SetOutputScalarTypeToUnsignedChar()
        elif (images["data_type"] in ("unsigned short", "uint16")):
            shifter.SetScale(float(2**16-1)/(global_max-global_min))
            shifter.SetOutputScalarTypeToUnsignedShort()
        elif (images["data_type"] in ("unsigned int", "uint32")):
            shifter.SetScale(float(2**32-1)/(global_max-global_min))
            shifter.SetOutputScalarTypeToUnsignedInt()
        elif (images["data_type"] in ("unsigned long", "uint64")):
            shifter.SetScale(float(2**64-1)/(global_max-global_min))
            shifter.SetOutputScalarTypeToUnsignedLong()
        elif (images["data_type"] in ("unsigned float", "ufloat")):
            shifter.SetScale(1./(global_max-global_min))
            shifter.SetOutputScalarTypeToFloat()
        for k_frame in xrange(images["n_frames"]):
            vtk_image = myVTK.readImage(
                filename=images["folder"]+"/"+images["basename"]+"_"+str(k_frame).zfill(images["zfill"])+".vti",
                verbose=verbose-1)
            if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
                shifter.SetInputData(vtk_image)
            else:
                shifter.SetInput(vtk_image)
            shifter.Update()
            vtk_image = shifter.GetOutput()
            myVTK.writeImage(
                image=vtk_image,
                filename=images["folder"]+"/"+images["basename"]+"_"+str(k_frame).zfill(images["zfill"])+".vti",
                verbose=verbose-1)
    else:
        assert (0), "Wrong data type. Aborting."