]> Creatis software - clitk.git/blob - common/clitkDicomRTStruct2ImageFilter.cxx
bf515a423b460ec958a6dd624bea7587103212c2
[clitk.git] / common / clitkDicomRTStruct2ImageFilter.cxx
1 /*=========================================================================
2   Program:         vv http://www.creatis.insa-lyon.fr/rio/vv
3   Main authors :   XX XX XX
4
5   Authors belongs to:
6   - University of LYON           http://www.universite-lyon.fr/
7   - Léon Bérard cancer center    http://www.centreleonberard.fr
8   - CREATIS CNRS laboratory      http://www.creatis.insa-lyon.fr
9
10   This software is distributed WITHOUT ANY WARRANTY; without even
11   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12   PURPOSE.  See the copyright notices for more information.
13
14   It is distributed under dual licence
15   - BSD       http://www.opensource.org/licenses/bsd-license.php
16   - CeCILL-B  http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17
18   =========================================================================*/
19
20 #include <iterator>
21 #include <algorithm>
22
23 // clitk
24 #include "clitkDicomRTStruct2ImageFilter.h"
25 #include "clitkImageCommon.h"
26
27 // vtk
28 #include <vtkVersion.h>
29 #include <vtkPolyDataToImageStencil.h>
30 #include <vtkSmartPointer.h>
31 #include <vtkImageStencil.h>
32 #include <vtkLinearExtrusionFilter.h>
33 #include <vtkMetaImageWriter.h>
34
35
36 //--------------------------------------------------------------------
37 clitk::DicomRTStruct2ImageFilter::DicomRTStruct2ImageFilter()
38 {
39   mROI = NULL;
40   mWriteOutput = false;
41   mCropMask = true;
42 }
43 //--------------------------------------------------------------------
44
45
46 //--------------------------------------------------------------------
47 clitk::DicomRTStruct2ImageFilter::~DicomRTStruct2ImageFilter()
48 {
49
50 }
51 //--------------------------------------------------------------------
52
53
54 //--------------------------------------------------------------------
55 bool clitk::DicomRTStruct2ImageFilter::ImageInfoIsSet() const
56 {
57   return mSize.size() && mSpacing.size() && mOrigin.size();
58 }
59 //--------------------------------------------------------------------
60
61
62 //--------------------------------------------------------------------
63 void clitk::DicomRTStruct2ImageFilter::SetWriteOutputFlag(bool b)
64 {
65   mWriteOutput = b;
66 }
67 //--------------------------------------------------------------------
68
69
70 //--------------------------------------------------------------------
71 void clitk::DicomRTStruct2ImageFilter::SetROI(clitk::DicomRT_ROI * roi)
72 {
73   mROI = roi;
74 }
75 //--------------------------------------------------------------------
76
77
78 //--------------------------------------------------------------------
79 void clitk::DicomRTStruct2ImageFilter::SetCropMaskEnabled(bool b)
80 {
81   mCropMask = b;
82 }
83 //--------------------------------------------------------------------
84
85
86 //--------------------------------------------------------------------
87 void clitk::DicomRTStruct2ImageFilter::SetOutputImageFilename(std::string s)
88 {
89   mOutputFilename = s;
90   mWriteOutput = true;
91 }
92 //--------------------------------------------------------------------
93
94
95 //--------------------------------------------------------------------
96 void clitk::DicomRTStruct2ImageFilter::SetImage(vvImage::Pointer image)
97 {
98   if (image->GetNumberOfDimensions() != 3) {
99     std::cerr << "Error. Please provide a 3D image." << std::endl;
100     exit(0);
101   }
102   mSpacing.resize(3);
103   mOrigin.resize(3);
104   mSize.resize(3);
105   mDirection.resize(3);
106   for(unsigned int i=0; i<3; i++) {
107     mSpacing[i] = image->GetSpacing()[i];
108     mOrigin[i] = image->GetOrigin()[i];
109     mSize[i] = image->GetSize()[i];
110     mDirection[i].resize(3);
111     for(unsigned int j=0; j<3; j++)
112       mDirection[i][j] = image->GetDirection()[i][j];
113   }
114 }
115 //--------------------------------------------------------------------
116
117
118 //--------------------------------------------------------------------
119 void clitk::DicomRTStruct2ImageFilter::SetImageFilename(std::string f)
120 {
121   itk::ImageIOBase::Pointer header = clitk::readImageHeader(f);
122   if (header->GetNumberOfDimensions() < 3) {
123     std::cerr << "Error. Please provide a 3D image instead of " << f << std::endl;
124     exit(0);
125   }
126   if (header->GetNumberOfDimensions() > 3) {
127     std::cerr << "Warning dimension > 3 are ignored" << std::endl;
128   }
129   mSpacing.resize(3);
130   mOrigin.resize(3);
131   mSize.resize(3);
132   mDirection.resize(3);
133   for(unsigned int i=0; i<3; i++) {
134     mSpacing[i] = header->GetSpacing(i);
135     mOrigin[i] = header->GetOrigin(i);
136     mSize[i] = header->GetDimensions(i);
137     mDirection[i].resize(3);
138     for(unsigned int j=0; j<3; j++)
139       mDirection[i][j] = header->GetDirection(i)[j];
140   }
141 }
142 //--------------------------------------------------------------------
143
144
145 //--------------------------------------------------------------------
146 void clitk::DicomRTStruct2ImageFilter::SetOutputOrigin(const double* origin)
147 {
148   std::copy(origin,origin+3,std::back_inserter(mOrigin));
149 }
150 //--------------------------------------------------------------------
151
152
153 //--------------------------------------------------------------------
154 void clitk::DicomRTStruct2ImageFilter::SetOutputSpacing(const double* spacing)
155 {
156   std::copy(spacing,spacing+3,std::back_inserter(mSpacing));
157 }
158 //--------------------------------------------------------------------
159
160
161 //--------------------------------------------------------------------
162 void clitk::DicomRTStruct2ImageFilter::SetOutputSize(const unsigned long* size)
163 {
164   std::copy(size,size+3,std::back_inserter(mSize));
165 }
166 //--------------------------------------------------------------------
167
168
169 //--------------------------------------------------------------------
170 void clitk::DicomRTStruct2ImageFilter::Update()
171 {
172   if (!mROI) {
173     std::cerr << "Error. No ROI set, please use SetROI." << std::endl;
174     exit(0);
175   }
176   if (!ImageInfoIsSet()) {
177     std::cerr << "Error. Please provide image info (spacing/origin) with SetImageFilename" << std::endl;
178     exit(0);
179   }
180
181   // Get Mesh
182   vtkPolyData * mesh = mROI->GetMesh();  
183
184   // Get bounds
185   double *bounds=mesh->GetBounds();
186
187   //Change mOrigin, mSize and mSpacing with respect to the directions
188   // Spacing is influenced by input direction
189   std::vector<double> tempSpacing;
190   tempSpacing.resize(3);
191   for(int i=0; i<3; i++) {
192     tempSpacing[i] = 0.0;
193     for(int j=0; j<3; j++) {
194       tempSpacing[i] += mDirection[i][j] * mSpacing[j];
195     }
196   }
197   for(int i=0; i<3; i++)
198     mSpacing[i] = tempSpacing[i];
199
200   // Size is influenced by affine transform matrix and input direction
201   // Size is converted to double, transformed and converted back to size type.
202   std::vector<double> tempSize;
203   tempSize.resize(3);
204   for(int i=0; i<3; i++) {
205     tempSize[i] = 0.0;
206     for(int j=0; j<3; j++) {
207       tempSize[i] += mDirection[i][j] * mSize[j];
208     }
209   }
210   for(int i=0; i<3; i++) {
211     if (tempSize[i] < 0.0) {
212       tempSize[i] *= -1;
213       mOrigin[i] += mSpacing[i]*(tempSize[i] -1);
214       mSpacing[i] *= -1;
215     }
216     mSize[i] = lrint(tempSize[i]);
217   }
218
219   // Compute origin
220   std::vector<double> origin;
221   origin.resize(3);
222   origin[0] = floor((bounds[0]-mOrigin[0])/mSpacing[0]-2)*mSpacing[0]+mOrigin[0];
223   origin[1] = floor((bounds[2]-mOrigin[1])/mSpacing[1]-2)*mSpacing[1]+mOrigin[1];
224   origin[2] = floor((bounds[4]-mOrigin[2])/mSpacing[2]-2)*mSpacing[2]+mOrigin[2];
225
226   // Compute extend
227   std::vector<double> extend;
228   extend.resize(3);
229   extend[0] = ceil((bounds[1]-origin[0])/mSpacing[0]+4);
230   extend[1] = ceil((bounds[3]-origin[1])/mSpacing[1]+4);
231   extend[2] = ceil((bounds[5]-origin[2])/mSpacing[2]+4);
232
233   // If no crop, set initial image size/origin
234   if (!mCropMask) {
235     for(int i=0; i<3; i++) {
236       origin[i] = mOrigin[i];
237       extend[i] = mSize[i]-1;
238     }
239   }
240
241   // Create new output image
242   mBinaryImage = vtkSmartPointer<vtkImageData>::New();
243 #if VTK_MAJOR_VERSION <= 5
244   mBinaryImage->SetScalarTypeToUnsignedChar();
245 #endif
246   mBinaryImage->SetOrigin(&origin[0]);
247   mBinaryImage->SetSpacing(&mSpacing[0]);
248   mBinaryImage->SetExtent(0, extend[0],
249                           0, extend[1],
250                           0, extend[2]);
251 #if VTK_MAJOR_VERSION <= 5
252   mBinaryImage->AllocateScalars();
253 #else
254   mBinaryImage->AllocateScalars(VTK_UNSIGNED_CHAR, 1);
255 #endif
256
257   memset(mBinaryImage->GetScalarPointer(), 0,
258          mBinaryImage->GetDimensions()[0]*mBinaryImage->GetDimensions()[1]*mBinaryImage->GetDimensions()[2]*sizeof(unsigned char));
259
260   // Extrude
261   vtkSmartPointer<vtkLinearExtrusionFilter> extrude=vtkSmartPointer<vtkLinearExtrusionFilter>::New();
262 #if VTK_MAJOR_VERSION <= 5
263   extrude->SetInput(mesh);
264 #else
265   extrude->SetInputData(mesh);
266 #endif
267   ///We extrude in the -slice_spacing direction to respect the FOCAL convention (NEEDED !)
268   extrude->SetVector(0, 0, -mSpacing[2]);
269
270   // Binarization
271   vtkSmartPointer<vtkPolyDataToImageStencil> sts=vtkSmartPointer<vtkPolyDataToImageStencil>::New();
272   //The following line is extremely important
273   //http://www.nabble.com/Bug-in-vtkPolyDataToImageStencil--td23368312.html#a23370933
274   sts->SetTolerance(0);
275   sts->SetInformationInput(mBinaryImage);
276 #if VTK_MAJOR_VERSION <= 5
277   sts->SetInput(extrude->GetOutput());
278 #else
279   sts->SetInputConnection(extrude->GetOutputPort(0));
280 #endif
281   //sts->SetInput(mesh);
282
283   vtkSmartPointer<vtkImageStencil> stencil=vtkSmartPointer<vtkImageStencil>::New();
284 #if VTK_MAJOR_VERSION <= 5
285   stencil->SetStencil(sts->GetOutput());
286 #else
287   stencil->SetStencilConnection(sts->GetOutputPort(0));
288 #endif
289 #if VTK_MAJOR_VERSION <= 5
290   stencil->SetInput(mBinaryImage);
291 #else
292   stencil->SetInputData(mBinaryImage);
293 #endif
294   stencil->ReverseStencilOn();
295   stencil->Update();
296
297   /*
298   vtkSmartPointer<vtkMetaImageWriter> w = vtkSmartPointer<vtkMetaImageWriter>::New();
299   w->SetInput(stencil->GetOutput());
300   w->SetFileName("binary2.mhd");
301   w->Write();
302   */
303
304   mBinaryImage->ShallowCopy(stencil->GetOutput());
305
306   if (mWriteOutput) {
307     typedef itk::Image<unsigned char, 3> ImageType;
308     typedef itk::VTKImageToImageFilter<ImageType> ConnectorType;
309     ConnectorType::Pointer connector = ConnectorType::New();
310     connector->SetInput(GetOutput());
311     connector->Update();
312     clitk::writeImage<ImageType>(connector->GetOutput(), mOutputFilename);
313   }
314 }
315 //--------------------------------------------------------------------
316
317
318
319 //--------------------------------------------------------------------
320 vtkImageData * clitk::DicomRTStruct2ImageFilter::GetOutput()
321 {
322   assert(mBinaryImage);
323   return mBinaryImage;
324 }
325 //--------------------------------------------------------------------
326