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 <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>
  47
  48 vvMesh::vvMesh() :
  49   r(1),g(0),b(0),
  50   slice_spacing(-1)
  51 {}
  52
  53 void vvMesh::AddMesh(vtkPolyData* p)
  54 {
  55   vtkPolyData * mesh=vtkPolyData::New();
  56   mesh->ShallowCopy(p);
  57   meshes.push_back(mesh);
  58 }
  59
  60 void vvMesh::ReadFromVTK(const char * filename)
  61 {
  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);
  67     r->Update();
  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);
  74     r->Update();
  75     assert(r->GetOutput());
  76     AddMesh(r->GetOutput());
  77   } else
  78     assert (false) ; //shouldn't happen
  79
  80   assert(GetNumberOfMeshes() != 0); ///We assume the object is empty
  81   structure_name=filename;
  82 }
  83
  84 void vvMesh::RemoveMeshes()
  85 {
  86   for (std::vector<vtkPolyData*>::const_iterator i=meshes.begin(); i!=meshes.end(); i++)
  87     (*i)->Delete();
  88   meshes.erase(meshes.begin(),meshes.end());
  89 }
  90
  91 void vvMesh::AddMask(vtkImageData* im)
  92 {
  93   assert(im->GetScalarType() == VTK_UNSIGNED_CHAR);
  94   vtkImageData* image=vtkImageData::New();
  95   image->ShallowCopy(im);
  96   masks.push_back(image);
  97 }
  98
  99 void vvMesh::RemoveMasks()
 100 {
 101   for (std::vector<vtkImageData*>::const_iterator i=masks.begin(); i!=masks.end(); i++)
 102     (*i)->Delete();
 103   masks.erase(masks.begin(),masks.end());
 104 }
 105
 106 vvMesh::~vvMesh()
 107 {
 108   RemoveMeshes();
 109   RemoveMasks();
 110 }
 111
 112 void vvMesh::CopyInformation(vvMesh::Pointer input)
 113 {
 114   r=input->r;
 115   g=input->g;
 116   b=input->b;
 117   structure_name=input->structure_name;
 118   slice_spacing=input->slice_spacing;
 119 }
 120
 121 void vvMesh::Print() const
 122 {
 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);
 127   }
 128   std::cout << "-------------------------" << std::endl << std::endl;
 129 }
 130
 131 void vvMesh::ComputeMasks(vtkImageData* sample,bool extrude)
 132 {
 133   this->RemoveMasks();
 134   for (std::vector<vtkPolyData*>::iterator i=meshes.begin(); i!=meshes.end(); i++) {
 135     vtkPolyData* mesh=*i;
 136     double *bounds=mesh->GetBounds();
 137
 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);
 146
 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();
 156 #else
 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);
 162
 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);
 172 #endif
 173     memset(binary_image->GetScalarPointer(),0,binary_image->GetDimensions()[0]*binary_image->GetDimensions()[1]*binary_image->GetDimensions()[2]*sizeof(unsigned char));
 174
 175
 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);
 181
 182     if (extrude) {
 183       vtkSmartPointer<vtkLinearExtrusionFilter> extrude=vtkSmartPointer<vtkLinearExtrusionFilter>::New();
 184 #if VTK_MAJOR_VERSION <= 5
 185       extrude->SetInput(mesh);
 186 #else
 187       extrude->SetInputData(mesh);
 188 #endif
 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());
 193 #else
 194       sts->SetInputConnection(extrude->GetOutputPort());
 195 #endif
 196     } else {
 197 #if VTK_MAJOR_VERSION <= 5
 198       sts->SetInput(mesh);
 199 #else
 200       sts->SetInputData(mesh);
 201 #endif
 202    }
 203
 204     vtkSmartPointer<vtkImageStencil> stencil=vtkSmartPointer<vtkImageStencil>::New();
 205 #if VTK_MAJOR_VERSION <= 5
 206     stencil->SetStencil(sts->GetOutput());
 207     stencil->SetInput(binary_image);
 208 #else
 209     stencil->SetStencilData(sts->GetOutput());
 210     stencil->SetInputData(binary_image);
 211 #endif
 212     stencil->Update();
 213     this->AddMask(stencil->GetOutput());
 214
 215     /*
 216       vtkSmartPointer<vtkMetaImageWriter> w = vtkSmartPointer<vtkMetaImageWriter>::New();
 217       w->SetInput(stencil->GetOutput());
 218       w->SetFileName("binary.mhd");
 219       w->Write();
 220     */
 221   }
 222 }
 223
 224 void vvMesh::ComputeMeshes()
 225 {
 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);
 231 #else
 232     marching->SetInputData(*i);
 233 #endif
 234     marching->SetValue(0,0.5);
 235     marching->Update();
 236     this->AddMesh(marching->GetOutput());
 237   }
 238 }
 239
 240 void vvMesh::propagateContour(vvImage::Pointer vf)
 241 {
 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();
 247
 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]);
 272     }
 273     this->AddMesh(new_mesh);
 274   }
 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);
 279   }
 280 }