Revert chirality of ill-transformed Mesh.

"""Transforms mesh volume from scanner coordinates to

image space and saves. """

#TODO: Merge _03_meshSurf_toImageSpace and _04_meshVol_toImagespace?

import numpy as numpy

import os

import declarations as dcl

from pom_funkce_VTK import readVTK, writeVTK

#

# If  result  does not  match, adjust  matrix_poz  or

# matrix_neg  and  shift  along  z  axis.

#

if __name__ == '__main__':

    import declarations as dcl

    mesh_folder = dcl.mesh_folder # or change to e.g.: "mesh2/"

    series_tuple = dcl.series_tuple

    endo_meshVol = 'Mesh_inImageSpace.vtk'

    vti_folder_prefix = dcl.vti_folder_prefix

    base = dcl.base

    matrix_file = dcl.matrix_file

else:

    import declarations as dcl

    mesh_folder = dcl.mesh_folder # never change this block

    series_tuple = dcl.series_tuple

    endo_meshVol = dcl.endo_meshVol

    vti_folder_prefix = dcl.vti_folder_prefix

    base = dcl.base

    matrix_file = dcl.matrix_file

    # from declarations import * # shorter alternative

info_folder = vti_folder_prefix + str(series_tuple[0]) + os.sep

offset_file = "offset.txt"

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

o = numpy.loadtxt( info_folder + offset_file )

T = numpy.loadtxt( info_folder + matrix_file)

#T_orig_wrong  =  [T[1],T[0],T[2]]   #  if  matrix_poz.txt initially

T_orig_wrong  =  [T[1],T[0],-T[2]]   #  if  matrix_neg.txt initially

T_orig_wrong_inv  =  numpy.linalg.inv(T_orig_wrong)

mesh = readVTK( mesh_folder + endo_meshVol)

points = mesh.GetPoints()

dataArray = points.GetData()

nodes = numpy.zeros((mesh.GetNumberOfPoints(),3))

for i in range (0,mesh.GetNumberOfPoints()):

    for j in range(0,3):

        nodes[i,j] = dataArray.GetComponent(i,j)

#nodes = nodes - o

#  Transform back to  physical  space.

for i in range(0,mesh.GetNumberOfPoints()):

    nodes[i,:] = T_orig_wrong_inv.dot(nodes[i,:])

#   Transform  to  the correct  image  space.

for i in range(0,mesh.GetNumberOfPoints()):

    nodes[i,:] = T.dot(nodes[i,:])

nodes =  nodes  + [0, 0, 2.5*(36 - 1)]  # shift  along the  z  axis

for i in range (0,mesh.GetNumberOfPoints()):

    for j in range(0,3):

        dataArray.SetComponent(i,j, nodes[i,j])

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

writeVTK(mesh, mesh_folder + endo_meshVol)

print( mesh_folder + endo_meshVol)