1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
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
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.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
25 #include "clitkCommon.h"
31 #include <itksys/SystemTools.hxx>
34 #include <vtkVersion.h>
35 #include <vtkSmartPointer.h>
36 #include <vtkFloatArray.h>
37 #include <vtkPointData.h>
38 #include <vtkPolyData.h>
39 #include <vtkPolyDataReader.h>
40 #include <vtkOBJReader.h>
41 #include <vtkImageData.h>
42 #include <vtkImageStencil.h>
43 #include <vtkLinearExtrusionFilter.h>
44 #include <vtkPolyDataToImageStencil.h>
45 #include <vtkMarchingCubes.h>
46 #include <vtkMetaImageWriter.h>
53 void vvMesh::AddMesh(vtkPolyData* p)
55 vtkPolyData * mesh=vtkPolyData::New();
57 meshes.push_back(mesh);
60 void vvMesh::ReadFromVTK(const char * filename)
62 std::string extension=itksys::SystemTools::GetFilenameLastExtension(std::string(filename));
63 if (extension == ".vtk" || extension== ".VTK") {
64 assert(GetNumberOfMeshes() == 0); ///We assume the object is empty
65 vtkSmartPointer<vtkPolyDataReader> r=vtkSmartPointer<vtkPolyDataReader>::New();
66 r->SetFileName(filename);
68 assert(r->GetOutput());
69 AddMesh(r->GetOutput());
70 } else if (extension == ".obj" || extension== ".OBJ") {
71 assert(GetNumberOfMeshes() == 0); ///We assume the object is empty
72 vtkSmartPointer<vtkOBJReader> r=vtkSmartPointer<vtkOBJReader>::New();
73 r->SetFileName(filename);
75 assert(r->GetOutput());
76 AddMesh(r->GetOutput());
78 assert (false) ; //shouldn't happen
80 assert(GetNumberOfMeshes() != 0); ///We assume the object is empty
81 structure_name=filename;
84 void vvMesh::RemoveMeshes()
86 for (std::vector<vtkPolyData*>::const_iterator i=meshes.begin(); i!=meshes.end(); i++)
88 meshes.erase(meshes.begin(),meshes.end());
91 void vvMesh::AddMask(vtkImageData* im)
93 assert(im->GetScalarType() == VTK_UNSIGNED_CHAR);
94 vtkImageData* image=vtkImageData::New();
95 image->ShallowCopy(im);
96 masks.push_back(image);
99 void vvMesh::RemoveMasks()
101 for (std::vector<vtkImageData*>::const_iterator i=masks.begin(); i!=masks.end(); i++)
103 masks.erase(masks.begin(),masks.end());
112 void vvMesh::CopyInformation(vvMesh::Pointer input)
117 structure_name=input->structure_name;
118 slice_spacing=input->slice_spacing;
121 void vvMesh::Print() const
123 std::cout << this << " : " << structure_name << std::endl << "RGB: " << r << "," << g << "," << b << std::endl;
124 for (std::vector<vtkPolyData*>::const_iterator i=meshes.begin(); i!=meshes.end(); i++) {
125 std::cout << (*i)->GetNumberOfPoints() << " points, " << (*i)->GetNumberOfCells() << " cells." << std::endl;
126 DDV((*i)->GetBounds(),6);
128 std::cout << "-------------------------" << std::endl << std::endl;
131 void vvMesh::ComputeMasks(vtkImageData* sample,bool extrude)
134 for (std::vector<vtkPolyData*>::iterator i=meshes.begin(); i!=meshes.end(); i++) {
135 vtkPolyData* mesh=*i;
136 double *bounds=mesh->GetBounds();
138 vtkSmartPointer<vtkImageData> binary_image=vtkSmartPointer<vtkImageData>::New();
139 #if VTK_MAJOR_VERSION <= 5
140 binary_image->SetScalarTypeToUnsignedChar();
141 ///Use the smallest mask in which the mesh fits
142 // Add two voxels on each side to make sure the mesh fits
143 double * samp_origin=sample->GetOrigin();
144 double * spacing=sample->GetSpacing();
145 binary_image->SetSpacing(spacing);
147 /// Put the origin on a voxel to avoid small skips
148 binary_image->SetOrigin(floor((bounds[0]-samp_origin[0])/spacing[0]-2)*spacing[0]+samp_origin[0],
149 floor((bounds[2]-samp_origin[1])/spacing[1]-2)*spacing[1]+samp_origin[1],
150 floor((bounds[4]-samp_origin[2])/spacing[2]-2)*spacing[2]+samp_origin[2]);
151 double * origin=binary_image->GetOrigin();
152 binary_image->SetExtent(0,ceil((bounds[1]-origin[0])/spacing[0]+4),
153 0,ceil((bounds[3]-origin[1])/spacing[1]+4),
154 0,ceil((bounds[5]-origin[2])/spacing[2])+4);
155 binary_image->AllocateScalars();
157 ///Use the smallest mask in which the mesh fits
158 // Add two voxels on each side to make sure the mesh fits
159 double * samp_origin=sample->GetOrigin();
160 double * spacing=sample->GetSpacing();
161 binary_image->SetSpacing(spacing);
163 /// Put the origin on a voxel to avoid small skips
164 binary_image->SetOrigin(floor((bounds[0]-samp_origin[0])/spacing[0]-2)*spacing[0]+samp_origin[0],
165 floor((bounds[2]-samp_origin[1])/spacing[1]-2)*spacing[1]+samp_origin[1],
166 floor((bounds[4]-samp_origin[2])/spacing[2]-2)*spacing[2]+samp_origin[2]);
167 double * origin=binary_image->GetOrigin();
168 binary_image->SetExtent(0,ceil((bounds[1]-origin[0])/spacing[0]+4),
169 0,ceil((bounds[3]-origin[1])/spacing[1]+4),
170 0,ceil((bounds[5]-origin[2])/spacing[2])+4);
171 binary_image->AllocateScalars(VTK_UNSIGNED_CHAR, 1);
173 memset(binary_image->GetScalarPointer(),0,binary_image->GetDimensions()[0]*binary_image->GetDimensions()[1]*binary_image->GetDimensions()[2]*sizeof(unsigned char));
176 vtkSmartPointer<vtkPolyDataToImageStencil> sts=vtkSmartPointer<vtkPolyDataToImageStencil>::New();
177 //The following line is extremely important
178 //http://www.nabble.com/Bug-in-vtkPolyDataToImageStencil--td23368312.html#a23370933
179 sts->SetTolerance(0);
180 sts->SetInformationInput(binary_image);
183 vtkSmartPointer<vtkLinearExtrusionFilter> extrude=vtkSmartPointer<vtkLinearExtrusionFilter>::New();
184 #if VTK_MAJOR_VERSION <= 5
185 extrude->SetInput(mesh);
187 extrude->SetInputData(mesh);
189 ///We extrude in the -slice_spacing direction to respect the FOCAL convention
190 extrude->SetVector(0, 0, -slice_spacing);
191 #if VTK_MAJOR_VERSION <= 5
192 sts->SetInput(extrude->GetOutput());
194 sts->SetInputConnection(extrude->GetOutputPort());
197 #if VTK_MAJOR_VERSION <= 5
200 sts->SetInputData(mesh);
204 vtkSmartPointer<vtkImageStencil> stencil=vtkSmartPointer<vtkImageStencil>::New();
205 #if VTK_MAJOR_VERSION <= 5
206 stencil->SetStencil(sts->GetOutput());
207 stencil->SetInput(binary_image);
209 stencil->SetStencilData(sts->GetOutput());
210 stencil->SetInputData(binary_image);
213 this->AddMask(stencil->GetOutput());
216 vtkSmartPointer<vtkMetaImageWriter> w = vtkSmartPointer<vtkMetaImageWriter>::New();
217 w->SetInput(stencil->GetOutput());
218 w->SetFileName("binary.mhd");
224 void vvMesh::ComputeMeshes()
226 this->RemoveMeshes();
227 for (std::vector<vtkImageData*>::iterator i=masks.begin(); i!=masks.end(); i++) {
228 vtkSmartPointer<vtkMarchingCubes> marching = vtkSmartPointer<vtkMarchingCubes>::New();
229 #if VTK_MAJOR_VERSION <= 5
230 marching->SetInput(*i);
232 marching->SetInputData(*i);
234 marching->SetValue(0,0.5);
236 this->AddMesh(marching->GetOutput());
240 void vvMesh::propagateContour(vvImage::Pointer vf)
242 assert(this->GetNumberOfMeshes() == 1);
243 std::vector<vtkImageData*> sgrids=vf->GetVTKImages();
244 vtkSmartPointer<vtkPolyData> reference_mesh = vtkSmartPointer<vtkPolyData>::New();
245 reference_mesh->ShallowCopy(this->GetMesh(0));
246 this->RemoveMeshes();
248 for (std::vector<vtkImageData*>::iterator i=sgrids.begin();
249 i!=sgrids.end(); i++) {
250 vtkPolyData* new_mesh=vtkPolyData::New();
251 new_mesh->DeepCopy(reference_mesh);
252 double Ox=vf->GetOrigin()[0];
253 double Oy=vf->GetOrigin()[1];
254 double Oz=vf->GetOrigin()[2];
255 double Sx=vf->GetSpacing()[0];
256 double Sy=vf->GetSpacing()[1];
257 double Sz=vf->GetSpacing()[2];
258 int *dims=vf->GetVTKImages()[0]->GetDimensions();
259 assert((*i)->GetScalarType() == VTK_FLOAT); //vfs are assumed to be of float type
260 assert((*i)->GetNumberOfScalarComponents() == 3);
261 float * vector_data=reinterpret_cast<float*>((*i)->GetScalarPointer());
262 for (int j=0; j<new_mesh->GetNumberOfPoints(); j++) {
263 double* old=new_mesh->GetPoint(j);
264 int ix=(old[0]-Ox)/Sx;
265 int iy=(old[1]-Oy)/Sy;
266 int iz=(old[2]-Oz)/Sz;
267 float* vector=vector_data+(ix+iy*vf->GetSize()[0]+iz*vf->GetSize()[0]*vf->GetSize()[1])*3;
268 if (ix>=0 && ix < dims[0]
269 && iy>=0 && iy < dims[1]
270 && iz>=0 && iz < dims[2])
271 new_mesh->GetPoints()->SetPoint(j,old[0]+vector[0],old[1]+vector[1],old[2]+vector[2]);
273 this->AddMesh(new_mesh);
275 if (GetNumberOfMasks()) { //If the input mesh has a mask, use it to compute the warped meshes' masks
276 vtkSmartPointer<vtkImageData> ref_mask = vtkSmartPointer<vtkImageData>::New();
277 ref_mask->ShallowCopy(GetMask(0));
278 this->ComputeMasks(ref_mask);