]> Creatis software - clitk.git/blob - vv/vvMesh.cxx
Merge branch 'master' of git.creatis.insa-lyon.fr:clitk
[clitk.git] / vv / vvMesh.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
19 //std 
20 #include <sstream>
21 #include <cassert>
22 #include <vector>
23
24 // clitk
25 #include "clitkCommon.h"
26
27 // vv
28 #include "vvMesh.h"
29
30 // itk
31 #include <itksys/SystemTools.hxx>
32
33 // vtk
34 #include <vtkSmartPointer.h>
35 #include <vtkFloatArray.h>
36 #include <vtkPointData.h>
37 #include <vtkPolyData.h>
38 #include <vtkPolyDataReader.h>
39 #include <vtkOBJReader.h>
40 #include <vtkImageData.h>
41 #include <vtkImageStencil.h>
42 #include <vtkLinearExtrusionFilter.h>
43 #include <vtkPolyDataToImageStencil.h>
44 #include <vtkMarchingCubes.h>
45 #include <vtkMetaImageWriter.h>
46
47 vvMesh::vvMesh() :
48   r(1),g(0),b(0),
49   slice_spacing(-1)
50 {}
51
52 void vvMesh::AddMesh(vtkPolyData* p)
53 {
54   vtkPolyData * mesh=vtkPolyData::New();
55   mesh->ShallowCopy(p);
56   meshes.push_back(mesh);
57 }
58
59 void vvMesh::ReadFromVTK(const char * filename)
60 {
61   std::string extension=itksys::SystemTools::GetFilenameLastExtension(std::string(filename));
62   if (extension == ".vtk" || extension== ".VTK") {
63     assert(GetNumberOfMeshes() == 0); ///We assume the object is empty
64     vtkSmartPointer<vtkPolyDataReader> r=vtkSmartPointer<vtkPolyDataReader>::New();
65     r->SetFileName(filename);
66     r->Update();
67     assert(r->GetOutput());
68     AddMesh(r->GetOutput());
69   } else if (extension == ".obj" || extension== ".OBJ") {
70     assert(GetNumberOfMeshes() == 0); ///We assume the object is empty
71     vtkSmartPointer<vtkOBJReader> r=vtkSmartPointer<vtkOBJReader>::New();
72     r->SetFileName(filename);
73     r->Update();
74     assert(r->GetOutput());
75     AddMesh(r->GetOutput());
76   } else
77     assert (false) ; //shouldn't happen
78
79   assert(GetNumberOfMeshes() != 0); ///We assume the object is empty
80   structure_name=filename;
81 }
82
83 void vvMesh::RemoveMeshes()
84 {
85   for (std::vector<vtkPolyData*>::const_iterator i=meshes.begin(); i!=meshes.end(); i++)
86     (*i)->Delete();
87   meshes.erase(meshes.begin(),meshes.end());
88 }
89
90 void vvMesh::AddMask(vtkImageData* im)
91 {
92   assert(im->GetScalarType() == VTK_UNSIGNED_CHAR);
93   vtkImageData* image=vtkImageData::New();
94   image->ShallowCopy(im);
95   masks.push_back(image);
96 }
97
98 void vvMesh::RemoveMasks()
99 {
100   for (std::vector<vtkImageData*>::const_iterator i=masks.begin(); i!=masks.end(); i++)
101     (*i)->Delete();
102   masks.erase(masks.begin(),masks.end());
103 }
104
105 vvMesh::~vvMesh()
106 {
107   RemoveMeshes();
108   RemoveMasks();
109 }
110
111 void vvMesh::CopyInformation(vvMesh::Pointer input)
112 {
113   r=input->r;
114   g=input->g;
115   b=input->b;
116   structure_name=input->structure_name;
117   slice_spacing=input->slice_spacing;
118 }
119
120 void vvMesh::Print() const
121 {
122   std::cout << this << " : " << structure_name << std::endl << "RGB: " << r << "," << g << "," << b << std::endl;
123   for (std::vector<vtkPolyData*>::const_iterator i=meshes.begin(); i!=meshes.end(); i++) {
124     std::cout << (*i)->GetNumberOfPoints() << " points, " << (*i)->GetNumberOfCells() << " cells." << std::endl;
125     DDV((*i)->GetBounds(),6);
126   }
127   std::cout << "-------------------------" << std::endl << std::endl;
128 }
129
130 void vvMesh::ComputeMasks(vtkImageData* sample,bool extrude)
131 {
132   this->RemoveMasks();
133   for (std::vector<vtkPolyData*>::iterator i=meshes.begin(); i!=meshes.end(); i++) {
134     vtkPolyData* mesh=*i;
135     double *bounds=mesh->GetBounds();
136
137     vtkSmartPointer<vtkImageData> binary_image=vtkSmartPointer<vtkImageData>::New();
138     binary_image->SetScalarTypeToUnsignedChar();
139     ///Use the smallest mask in which the mesh fits
140     // Add two voxels on each side to make sure the mesh fits
141     double * samp_origin=sample->GetOrigin();
142     double * spacing=sample->GetSpacing();
143     binary_image->SetSpacing(spacing);
144
145     /// Put the origin on a voxel to avoid small skips
146     binary_image->SetOrigin(floor((bounds[0]-samp_origin[0])/spacing[0]-2)*spacing[0]+samp_origin[0],
147                             floor((bounds[2]-samp_origin[1])/spacing[1]-2)*spacing[1]+samp_origin[1],
148                             floor((bounds[4]-samp_origin[2])/spacing[2]-2)*spacing[2]+samp_origin[2]);
149     double * origin=binary_image->GetOrigin();
150     binary_image->SetExtent(0,ceil((bounds[1]-origin[0])/spacing[0]+4),
151                             0,ceil((bounds[3]-origin[1])/spacing[1]+4),
152                             0,ceil((bounds[5]-origin[2])/spacing[2])+4);
153     binary_image->AllocateScalars();
154     memset(binary_image->GetScalarPointer(),0,binary_image->GetDimensions()[0]*binary_image->GetDimensions()[1]*binary_image->GetDimensions()[2]*sizeof(unsigned char));
155
156
157     vtkSmartPointer<vtkPolyDataToImageStencil> sts=vtkSmartPointer<vtkPolyDataToImageStencil>::New();
158     //The following line is extremely important
159     //http://www.nabble.com/Bug-in-vtkPolyDataToImageStencil--td23368312.html#a23370933
160     sts->SetTolerance(0);
161     sts->SetInformationInput(binary_image);
162
163     if (extrude) {
164       vtkSmartPointer<vtkLinearExtrusionFilter> extrude=vtkSmartPointer<vtkLinearExtrusionFilter>::New();
165       extrude->SetInput(mesh);
166       ///We extrude in the -slice_spacing direction to respect the FOCAL convention
167       extrude->SetVector(0, 0, -slice_spacing);
168       sts->SetInput(extrude->GetOutput());
169     } else
170       sts->SetInput(mesh);
171
172     vtkSmartPointer<vtkImageStencil> stencil=vtkSmartPointer<vtkImageStencil>::New();
173     stencil->SetStencil(sts->GetOutput());
174     stencil->SetInput(binary_image);
175     stencil->Update();
176     this->AddMask(stencil->GetOutput());
177
178     /*
179       vtkSmartPointer<vtkMetaImageWriter> w = vtkSmartPointer<vtkMetaImageWriter>::New();
180       w->SetInput(stencil->GetOutput());
181       w->SetFileName("binary.mhd");
182       w->Write();
183     */
184   }
185 }
186
187 void vvMesh::ComputeMeshes()
188 {
189   this->RemoveMeshes();
190   for (std::vector<vtkImageData*>::iterator i=masks.begin(); i!=masks.end(); i++) {
191     vtkSmartPointer<vtkMarchingCubes> marching = vtkSmartPointer<vtkMarchingCubes>::New();
192     marching->SetInput(*i);
193     marching->SetValue(0,0.5);
194     marching->Update();
195     this->AddMesh(marching->GetOutput());
196   }
197 }
198
199 void vvMesh::propagateContour(vvImage::Pointer vf)
200 {
201   assert(this->GetNumberOfMeshes() == 1);
202   std::vector<vtkImageData*> sgrids=vf->GetVTKImages();
203   vtkSmartPointer<vtkPolyData> reference_mesh = vtkSmartPointer<vtkPolyData>::New();
204   reference_mesh->ShallowCopy(this->GetMesh(0));
205   this->RemoveMeshes();
206
207   for (std::vector<vtkImageData*>::iterator i=sgrids.begin();
208        i!=sgrids.end(); i++) {
209     vtkPolyData* new_mesh=vtkPolyData::New();
210     new_mesh->DeepCopy(reference_mesh);
211     double Ox=vf->GetOrigin()[0];
212     double Oy=vf->GetOrigin()[1];
213     double Oz=vf->GetOrigin()[2];
214     double Sx=vf->GetSpacing()[0];
215     double Sy=vf->GetSpacing()[1];
216     double Sz=vf->GetSpacing()[2];
217     int *dims=vf->GetVTKImages()[0]->GetDimensions();
218     assert((*i)->GetScalarType() == VTK_FLOAT); //vfs are assumed to be of float type
219     assert((*i)->GetNumberOfScalarComponents() == 3);
220     float * vector_data=reinterpret_cast<float*>((*i)->GetScalarPointer());
221     for (int j=0; j<new_mesh->GetNumberOfPoints(); j++) {
222       double* old=new_mesh->GetPoint(j);
223       int ix=(old[0]-Ox)/Sx;
224       int iy=(old[1]-Oy)/Sy;
225       int iz=(old[2]-Oz)/Sz;
226       float* vector=vector_data+(ix+iy*vf->GetSize()[0]+iz*vf->GetSize()[0]*vf->GetSize()[1])*3;
227       if (ix>=0 && ix < dims[0]
228           && iy>=0 && iy < dims[1]
229           && iz>=0 && iz < dims[2])
230         new_mesh->GetPoints()->SetPoint(j,old[0]+vector[0],old[1]+vector[1],old[2]+vector[2]);
231     }
232     this->AddMesh(new_mesh);
233   }
234   if (GetNumberOfMasks()) { //If the input mesh has a mask, use it to compute the warped meshes' masks
235     vtkSmartPointer<vtkImageData> ref_mask = vtkSmartPointer<vtkImageData>::New();
236     ref_mask->ShallowCopy(GetMask(0));
237     this->ComputeMasks(ref_mask);
238   }
239 }