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