vv/vvMesh.cxx

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