]> Creatis software - clitk.git/blob - tools/clitkTransformLandmarks.cxx
Update Landmark coordinates with transformation matrix
[clitk.git] / 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