tools/clitkTransformLandmarks.cxx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ===========================================================================*/
  18 #include "clitkTransformLandmarks_ggo.h"
  19
  20 #include "clitkTransformUtilities.h"
  21 #include <string>
  22 #include <vector>
  23 #include <iostream>
  24 #include <fstream>
  25
  26 #include <vtkSmartPointer.h>
  27 #include <vtkPoints.h>
  28 #include <vtkPolyData.h>
  29 #include "vtkPolyDataReader.h"
  30 #include "vtkPolyDataWriter.h"
  31 #include <vtkTransform.h>
  32 #include <vtkTransformFilter.h>
  33
  34 typedef itk::Matrix<double, 4, 4> MatrixType;
  35 typedef itk::Point<double, 4> PointType;
  36 typedef std::vector<PointType> PointArrayType;
  37 typedef itk::FixedArray<std::string, 2> TxtDataType;
  38 typedef std::vector<TxtDataType> TxtDataArrayType;
  39
  40
  41 void read_points_txt(const std::string& fileName, PointArrayType& points, TxtDataArrayType& data);
  42 void write_points_txt(const std::string& fileName, const PointArrayType& points, const TxtDataArrayType& data);
  43
  44 void read_points_pts(const std::string& fileName, PointArrayType& points);
  45 void write_points_pts(const std::string& fileName, const PointArrayType& points);
  46
  47 void apply_spacing(const PointArrayType& input, const double* spacing, PointArrayType& output);
  48 void transform_points(const PointArrayType& input, const MatrixType& matrix, PointArrayType& output);
  49
  50
  51 bool verbose = false;
  52
  53 int main(int argc, char** argv)
  54 {
  55   GGO(clitkTransformLandmarks, args_info);
  56   verbose = args_info.verbose_flag;
  57
  58   TxtDataArrayType data;
  59   PointArrayType inputPoints;
  60   if (strcmp(args_info.type_arg, "txt") == 0) {
  61     read_points_txt(args_info.input_arg, inputPoints, data);
  62   }
  63   else if (strcmp(args_info.type_arg, "vtk") == 0) {
  64     vtkSmartPointer<vtkPolyDataReader> reader = vtkSmartPointer<vtkPolyDataReader>::New();
  65     reader->SetFileName(args_info.input_arg);
  66     reader->Update();
  67     vtkSmartPointer<vtkPolyDataWriter> writer = vtkSmartPointer<vtkPolyDataWriter>::New();
  68     writer->SetFileName( args_info.output_arg );
  69
  70     if (args_info.matrix_given) {
  71       vtkSmartPointer<vtkTransformFilter> transformFilter = vtkSmartPointer<vtkTransformFilter>::New();
  72       vtkSmartPointer<vtkTransform> transform = vtkSmartPointer<vtkTransform>::New();
  73       vtkMatrix4x4* matrix = clitk::ReadVTKMatrix3D(args_info.matrix_arg);
  74       vtkSmartPointer<vtkMatrix4x4> matrixT = vtkSmartPointer<vtkMatrix4x4>::New();
  75       vtkMatrix4x4::Invert(matrix, matrixT); //not sure why, but this seems necessary for using the same .mat as when loading file with vv (probably due to the inversion trick performed in the vv reader...)
  76       transform->SetMatrix(matrixT);
  77       transformFilter->SetInputConnection(reader->GetOutputPort());
  78       transformFilter->SetTransform(transform);
  79       writer->SetInputConnection(transformFilter->GetOutputPort());
  80
  81     }
  82     else { //just write the output
  83       writer->SetInputConnection( reader->GetOutputPort() );
  84     }
  85
  86     writer->Write();
  87     return 0;
  88   }
  89   else {
  90     read_points_pts(args_info.input_arg, inputPoints);
  91   }
  92
  93   PointArrayType outputPoints;
  94   PointArrayType spacingPoints;
  95   PointArrayType* workingInputPoints = &inputPoints;
  96   PointArrayType* workingOutputPoints = &outputPoints;
  97   if (args_info.spacing_given) {
  98     if (verbose) std::cout << "Processing spacing..." << std::endl;
  99
 100     apply_spacing(*workingInputPoints, args_info.spacing_arg, spacingPoints);
 101     workingInputPoints = &spacingPoints;
 102     workingOutputPoints = &spacingPoints;
 103   }
 104
 105   MatrixType matrix;
 106   if (args_info.matrix_given) {
 107     matrix = clitk::ReadMatrix3D(args_info.matrix_arg);
 108     transform_points(*workingInputPoints, matrix, outputPoints);
 109     workingOutputPoints = &outputPoints;
 110   }
 111
 112   if (strcmp(args_info.type_arg, "txt") == 0) {
 113     write_points_txt(args_info.output_arg, *workingOutputPoints, data);
 114   }
 115   else {
 116     write_points_pts(args_info.output_arg, *workingOutputPoints);
 117   }
 118
 119   return 0;
 120 }
 121
 122 void read_points_txt(const std::string& fileName, PointArrayType& points, TxtDataArrayType& data)
 123 {
 124   std::ifstream landmarksFile(fileName.c_str());
 125   if (landmarksFile.fail()) {
 126     std::cerr << "ERROR: could not open '" << fileName << "'" << std::endl;
 127     exit(-2);
 128   }
 129
 130   std::string line;
 131   std::getline(landmarksFile, line);
 132   if (line.find("LANDMARKS") == std::string::npos) {
 133     std::cerr << "ERROR: invalid landmarks file '" << fileName << "'" << std::endl;
 134     exit(-3);
 135   }
 136
 137   PointType p;
 138   p.Fill(1);
 139
 140   TxtDataType d;
 141
 142   while (!landmarksFile.eof()) {
 143     // read id, x, y, z, d1, d2
 144     landmarksFile >> d[0] >> p[0] >> p[1] >> p[2] >> d[1];// >> d[2];
 145     if (landmarksFile.fail())
 146       break;
 147
 148     if (verbose){
 149       std::cout << "point " << p << std::endl;
 150       std::cout << "data " << " " << d[0] << " " << d[1] << std::endl;
 151     }
 152
 153     points.push_back(p);
 154     data.push_back(d);
 155   }
 156 }
 157
 158 void write_points_txt(const std::string& fileName, const PointArrayType& points, const TxtDataArrayType& data)
 159 {
 160   std::ofstream landmarksFile(fileName.c_str());
 161
 162   landmarksFile << "LANDMARKS1" << std::endl;
 163   for (size_t i = 0; i < points.size(); i++)
 164     landmarksFile << data[i][0] << " " << points[i][0] << " " << points[i][1] << " " << points[i][2] << " " << data[i][1] << " " << std::endl;
 165 }
 166
 167 void read_points_pts(const std::string& fileName, PointArrayType& points)
 168 {
 169   std::ifstream landmarksFile(fileName.c_str());
 170   if (landmarksFile.fail()) {
 171     std::cerr << "ERROR: could not open '" << fileName << "'" << std::endl;
 172     exit(-2);
 173   }
 174
 175   std::string line;
 176   std::getline(landmarksFile, line);
 177   if (line.find("#X") != 0) {
 178     std::cerr << "ERROR: invalid landmarks file '" << fileName << "'" << std::endl;
 179     exit(-3);
 180   }
 181
 182   PointType p;
 183   p.Fill(1);
 184
 185   while (!landmarksFile.eof()) {
 186     // read id, x, y, z, d1, d2
 187     landmarksFile >> p[0] >> p[1] >> p[2];
 188     if (landmarksFile.fail())
 189       break;
 190
 191     if (verbose){
 192       std::cout << "point " << p << std::endl;
 193     }
 194
 195     points.push_back(p);
 196   }
 197 }
 198
 199 void write_points_pts(const std::string& fileName, const PointArrayType& points)
 200 {
 201   std::ofstream landmarksFile(fileName.c_str());
 202
 203   landmarksFile << "#X\tY\tZ" << std::endl;
 204   for (size_t i = 0; i < points.size(); i++)
 205     landmarksFile << points[i][0] << "\t" << points[i][1] << "\t" << points[i][2] << "\t" << std::endl;
 206 }
 207
 208 void apply_spacing(const PointArrayType& input, const double* spacing, PointArrayType& output)
 209 {
 210   PointType out;
 211   out.Fill(1);
 212
 213   for (size_t i = 0; i < input.size(); i++) {
 214     out[0] = input[i][0] * spacing[0];
 215     out[1] = input[i][1] * spacing[1];
 216     out[2] = input[i][2] * spacing[2];
 217     if (verbose){
 218       std::cout << "output " << out << std::endl;
 219     }
 220     output.push_back(out);
 221   }
 222 }
 223
 224 void transform_points(const PointArrayType& input, const MatrixType& matrix, PointArrayType& output)
 225 {
 226   for (size_t i = 0; i < input.size(); i++) {
 227     output.push_back(matrix * input[i]);
 228   }
 229 }
 230