vv/vvMesh.cxx

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