1 /*=========================================================================
5 Author : Joel Schaerer (joel.schaerer@insa-lyon.fr)
8 Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
9 CREATIS-LRMN http://www.creatis.insa-lyon.fr
11 This program is free software: you can redistribute it and/or modify
12 it under the terms of the GNU General Public License as published by
13 the Free Software Foundation, version 3 of the License.
15 This program is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU General Public License for more details.
20 You should have received a copy of the GNU General Public License
21 along with this program. If not, see <http://www.gnu.org/licenses/>.
23 =========================================================================*/
28 #include <vtkSmartPointer.h>
29 #include <vtkFloatArray.h>
30 #include <vtkPointData.h>
31 #include <vtkPolyData.h>
32 #include <vtkPolyDataReader.h>
33 #include <vtkImageData.h>
34 #include "clitkCommon.h"
36 #include <vtkImageStencil.h>
37 #include <vtkLinearExtrusionFilter.h>
38 #include <vtkPolyDataToImageStencil.h>
39 #include <vtkMarchingCubes.h>
41 #include <vtkMetaImageWriter.h>
48 void vvMesh::AddMesh(vtkPolyData* p)
50 vtkPolyData * mesh=vtkPolyData::New();
52 meshes.push_back(mesh);
55 void vvMesh::ReadFromVTK(const char * filename)
57 assert(GetNumberOfMeshes() == 0); ///We assume the object is empty
58 vtkSmartPointer<vtkPolyDataReader> r=vtkSmartPointer<vtkPolyDataReader>::New();
59 r->SetFileName(filename);
61 AddMesh(r->GetOutput());
62 structure_name=filename;
65 void vvMesh::RemoveMeshes()
67 for (std::vector<vtkPolyData*>::const_iterator i=meshes.begin();i!=meshes.end();i++)
69 meshes.erase(meshes.begin(),meshes.end());
72 void vvMesh::AddMask(vtkImageData* im)
74 assert(im->GetScalarType() == VTK_UNSIGNED_CHAR);
75 vtkImageData* image=vtkImageData::New();
76 image->ShallowCopy(im);
77 masks.push_back(image);
80 void vvMesh::RemoveMasks()
82 for (std::vector<vtkImageData*>::const_iterator i=masks.begin();i!=masks.end();i++)
84 masks.erase(masks.begin(),masks.end());
93 void vvMesh::CopyInformation(vvMesh::Pointer input)
98 structure_name=input->structure_name;
99 slice_spacing=input->slice_spacing;
102 void vvMesh::Print() const
104 std::cout << this << " : " << structure_name << std::endl << "RGB: " << r << "," << g << "," << b << std::endl;
105 for (std::vector<vtkPolyData*>::const_iterator i=meshes.begin();i!=meshes.end();i++)
107 std::cout << (*i)->GetNumberOfPoints() << " points, " << (*i)->GetNumberOfCells() << " cells." << std::endl;
108 DDV((*i)->GetBounds(),6);
110 std::cout << "-------------------------" << std::endl << std::endl;
113 void vvMesh::ComputeMasks(vtkImageData* sample,bool extrude)
116 for (std::vector<vtkPolyData*>::iterator i=meshes.begin();i!=meshes.end();i++)
118 vtkPolyData* mesh=*i;
119 double *bounds=mesh->GetBounds();
121 vtkSmartPointer<vtkImageData> binary_image=vtkSmartPointer<vtkImageData>::New();
122 binary_image->SetScalarTypeToUnsignedChar();
123 ///Use the smallest mask in which the mesh fits
124 // Add two voxels on each side to make sure the mesh fits
125 double * samp_origin=sample->GetOrigin();
126 double * spacing=sample->GetSpacing();
127 binary_image->SetSpacing(spacing);
128 /// Put the origin on a voxel to avoid small skips
129 binary_image->SetOrigin(floor((bounds[0]-samp_origin[0])/spacing[0]-2)*spacing[0]+samp_origin[0],
130 floor((bounds[2]-samp_origin[1])/spacing[1]-2)*spacing[1]+samp_origin[1],
131 floor((bounds[4]-samp_origin[2])/spacing[2]-2)*spacing[2]+samp_origin[2]);
132 double * origin=binary_image->GetOrigin();
133 binary_image->SetExtent(0,ceil((bounds[1]-origin[0])/spacing[0]+4),
134 0,ceil((bounds[3]-origin[1])/spacing[1]+4),
135 0,ceil((bounds[5]-origin[2])/spacing[2])+4);
136 binary_image->AllocateScalars();
137 memset(binary_image->GetScalarPointer(),0,binary_image->GetDimensions()[0]*binary_image->GetDimensions()[1]*binary_image->GetDimensions()[2]*sizeof(unsigned char));
140 vtkSmartPointer<vtkPolyDataToImageStencil> sts=vtkSmartPointer<vtkPolyDataToImageStencil>::New();
141 //The following line is extremely important
142 //http://www.nabble.com/Bug-in-vtkPolyDataToImageStencil--td23368312.html#a23370933
143 sts->SetTolerance(0);
144 sts->SetInformationInput(binary_image);
148 vtkSmartPointer<vtkLinearExtrusionFilter> extrude=vtkSmartPointer<vtkLinearExtrusionFilter>::New();
149 extrude->SetInput(mesh);
150 ///We extrude in the -slice_spacing direction to respect the FOCAL convention
151 extrude->SetVector(0, 0, -slice_spacing);
152 sts->SetInput(extrude->GetOutput());
157 vtkSmartPointer<vtkImageStencil> stencil=vtkSmartPointer<vtkImageStencil>::New();
158 stencil->SetStencil(sts->GetOutput());
159 stencil->SetInput(binary_image);
161 this->AddMask(stencil->GetOutput());
162 //vtkSmartPointer<vtkMetaImageWriter> w = vtkSmartPointer<vtkMetaImageWriter>::New();
163 //w->SetInput(stencil->GetOutput());
164 //w->SetFileName("binary.mhd");
169 void vvMesh::ComputeMeshes()
171 this->RemoveMeshes();
172 for (std::vector<vtkImageData*>::iterator i=masks.begin();i!=masks.end();i++)
174 vtkSmartPointer<vtkMarchingCubes> marching = vtkSmartPointer<vtkMarchingCubes>::New();
175 marching->SetInput(*i);
176 marching->SetValue(0,0.5);
178 this->AddMesh(marching->GetOutput());
182 void vvMesh::propagateContour(vvImage::Pointer vf)
184 assert(this->GetNumberOfMeshes() == 1);
185 std::vector<vtkImageData*> sgrids=vf->GetVTKImages();
186 vtkSmartPointer<vtkPolyData> reference_mesh = vtkSmartPointer<vtkPolyData>::New();
187 reference_mesh->ShallowCopy(this->GetMesh(0));
188 this->RemoveMeshes();
190 for (std::vector<vtkImageData*>::iterator i=sgrids.begin();
193 vtkPolyData* new_mesh=vtkPolyData::New();
194 new_mesh->DeepCopy(reference_mesh);
195 double Ox=vf->GetOrigin()[0];
196 double Oy=vf->GetOrigin()[1];
197 double Oz=vf->GetOrigin()[2];
198 double Sx=vf->GetSpacing()[0];
199 double Sy=vf->GetSpacing()[1];
200 double Sz=vf->GetSpacing()[2];
201 int *dims=vf->GetVTKImages()[0]->GetDimensions();
202 assert((*i)->GetScalarType() == VTK_FLOAT); //vfs are assumed to be of float type
203 assert((*i)->GetNumberOfScalarComponents() == 3);
204 float * vector_data=reinterpret_cast<float*>((*i)->GetScalarPointer());
205 for (int j=0;j<new_mesh->GetNumberOfPoints();j++)
207 double* old=new_mesh->GetPoint(j);
208 int ix=(old[0]-Ox)/Sx;
209 int iy=(old[1]-Oy)/Sy;
210 int iz=(old[2]-Oz)/Sz;
211 float* vector=vector_data+(ix+iy*vf->GetSize()[0]+iz*vf->GetSize()[0]*vf->GetSize()[1])*3;
212 if (ix>=0 && ix < dims[0]
213 && iy>=0 && iy < dims[1]
214 && iz>=0 && iz < dims[2])
215 new_mesh->GetPoints()->SetPoint(j,old[0]+vector[0],old[1]+vector[1],old[2]+vector[2]);
217 this->AddMesh(new_mesh);
219 if (GetNumberOfMasks()) //If the input mesh has a mask, use it to compute the warped meshes' masks
221 vtkSmartPointer<vtkImageData> ref_mask = vtkSmartPointer<vtkImageData>::New();
222 ref_mask->ShallowCopy(GetMask(0));
223 this->ComputeMasks(ref_mask);