vv/vvMeshReader.cxx

   1 #include <algorithm>
   2
   3 #include <QApplication>
   4
   5 #include <gdcm.h>
   6 #include <gdcmSQItem.h>
   7
   8 #include <vtkSmartPointer.h>
   9 #include <vtkAppendPolyData.h>
  10 #include <vtkCellArray.h>
  11 #include <vtkMetaImageWriter.h>
  12 #include <vtkMetaImageReader.h>
  13 #include <vtkPolyDataWriter.h>
  14 #include <vtkImageData.h>
  15 #include <vtkImageCast.h>
  16 #include <vtkImageGaussianSmooth.h>
  17
  18 #include "clitkCommon.h"
  19 #include "vvMeshReader.h"
  20 #include "vvProgressDialog.h"
  21
  22 vvMeshReader::vvMeshReader() :
  23     vtk_mode(false)
  24 {}
  25
  26 void vvMeshReader::Update()
  27 {
  28     //Show a progress bar only when opening a DC-struct (ie. multiple contours)
  29     vvProgressDialog progress("Opening " + filename,(!vtk_mode) && (selected_contours.size()>1));
  30     this->start();
  31     while (this->isRunning())
  32     {
  33         progress.SetProgress(result.size(),selected_contours.size());
  34         this->wait(50);
  35         qApp->processEvents();
  36     }
  37 }
  38
  39 void vvMeshReader::run()
  40 {
  41     ///Verify the important stuff has been set
  42     assert(filename != "");
  43     if (vtk_mode) //Read vtkPolyData
  44     {
  45         vvMesh::Pointer m=vvMesh::New();
  46         m->ReadFromVTK(filename.c_str());
  47         if (vf) m->propagateContour(vf);
  48         m->ComputeMasks(image->GetVTKImages()[0],true);
  49         result.push_back(m);
  50     }
  51     else //Read a Dicom-struct file
  52     {
  53         assert(selected_contours.size() > 0);
  54         assert(image);
  55         std::vector<vvMesh::Pointer> contour_stacks=readSelectedContours();
  56         for (std::vector<vvMesh::Pointer>::iterator i=contour_stacks.begin();
  57                 i!=contour_stacks.end();i++)
  58         {
  59             (*i)->ComputeMasks(image->GetVTKImages()[0],true); //Remesh the contour
  60             (*i)->ComputeMeshes();
  61             if (vf) (*i)->propagateContour(vf);
  62             result.push_back(*i);
  63         }
  64     }
  65 }
  66
  67 template<class ElementType>
  68 ElementType parse_value(std::string str)
  69 {
  70     std::istringstream parser(str);
  71     ElementType value;
  72     parser >> value;
  73     assert(!parser.fail());
  74     return value;
  75 }
  76
  77 template<class ElementType>
  78 std::vector<ElementType> parse_string(std::string str,char delim)
  79 {
  80     std::istringstream ss(str);
  81     std::string token;
  82     std::vector<ElementType> result;
  83     while (getline(ss,token,delim))
  84     {
  85         result.push_back(parse_value<ElementType>(token));
  86     }
  87     return result;
  88 }
  89
  90 std::vector<std::pair<int,std::string> > vvMeshReader::GetROINames()
  91 {
  92     assert(filename!="");
  93     gdcm::File reader;
  94     reader.SetFileName(filename.c_str());
  95     reader.SetMaxSizeLoadEntry(16384);
  96     reader.Load();
  97
  98     gdcm::SeqEntry * roi_info=reader.GetSeqEntry(0x3006,0x0020);
  99     assert(roi_info);
 100     std::vector<std::pair<int, std::string> > roi_names;
 101     // DD("ici");
 102     //int n=0;
 103     for (gdcm::SQItem* i=roi_info->GetFirstSQItem();i!=0;i=roi_info->GetNextSQItem())
 104         if (i->GetEntryValue(0x3006,0x0022)!= gdcm::GDCM_UNFOUND)
 105             roi_names.push_back(make_pair(atoi(i->GetEntryValue(0x3006,0x0022).c_str()),i->GetEntryValue(0x3006,0x0026)));
 106     return roi_names;
 107 }
 108
 109 std::vector<vvMesh::Pointer> vvMeshReader::readSelectedContours()
 110 {
 111     gdcm::File reader;
 112     reader.SetFileName(filename.c_str());
 113     reader.SetMaxSizeLoadEntry(16384);
 114     reader.Load();
 115
 116     std::vector<vvMesh::Pointer> result;
 117     gdcm::SeqEntry * rois=reader.GetSeqEntry(0x3006,0x0039);
 118     ///We need to iterate both on the contours themselves, and on the contour info
 119     gdcm::SeqEntry * roi_info=reader.GetSeqEntry(0x3006,0x0020);
 120     gdcm::SQItem* k=roi_info->GetFirstSQItem();
 121     for(gdcm::SQItem* i=rois->GetFirstSQItem();i!=0;i=rois->GetNextSQItem()) //loop over ROIS
 122     {
 123         assert(k!=0);
 124         vtkSmartPointer<vtkAppendPolyData> append=vtkSmartPointer<vtkAppendPolyData>::New();
 125         std::istringstream ss(i->GetEntryValue(0x3006,0x0084));
 126         int roi_number;ss >> roi_number;
 127         if (std::find(selected_contours.begin(),selected_contours.end(),roi_number) != selected_contours.end())//Only read selected ROIs
 128         {
 129             vvMesh::Pointer current_roi=vvMesh::New();
 130             std::vector<double> rgb=parse_string<double>(i->GetEntryValue(0x3006,0x002a),'\\');
 131             assert(rgb.size()==3);
 132             current_roi->r=rgb[0]/255; current_roi->g=rgb[1]/255; current_roi->b=rgb[2]/255;
 133             current_roi->structure_name=k->GetEntryValue(0x3006,0x0026);
 134             gdcm::SeqEntry * contours=i->GetSeqEntry(0x3006,0x0040);
 135             double z0=-1; //Used to determine spacing between slices, assumed to be constant
 136             for(gdcm::SQItem* j=contours->GetFirstSQItem();j!=0;j=contours->GetNextSQItem()) //loop over 2D contours
 137             {
 138                 std::string contour_type=j->GetEntryValue(0x3006,0x0042);
 139                 if (contour_type=="CLOSED_PLANAR ")
 140                 {
 141                     int point_number=parse_value<int>(j->GetEntryValue(0x3006,0x0046));
 142                     std::vector<float> points=parse_string<float>(j->GetEntryValue(0x3006,0x0050),'\\');
 143                     assert(points.size() == static_cast<unsigned int>(point_number)*3);
 144                     if (z0 == -1) //First contour
 145                         z0=points[2];
 146                     //2nd contour, spacing not yet set. Need to be sure we are on a different slice,
 147                     //sometimes there is more than one closed contour per slice
 148                     else if (current_roi->GetSpacing()==-1 && points[2] != z0 )
 149                         current_roi->SetSpacing(points[2]-z0);
 150                     vtkPolyData * contour=vtkPolyData::New();
 151                     contour->Allocate(); //for cell structures
 152                     contour->SetPoints(vtkPoints::New());
 153                     vtkIdType ids[2];
 154                     for (unsigned int idx=0;idx<points.size();idx+=3)
 155                     {
 156                         contour->GetPoints()->InsertNextPoint(points[idx],points[idx+1],points[idx+2]);
 157                         ids[0]=idx/3;ids[1]=(ids[0]+1)%point_number; //0-1,1-2,...,n-1-0
 158                         contour->GetLines()->InsertNextCell(2,ids);
 159                     }
 160                     append->AddInput(contour);
 161                 }
 162                 else if (contour_type == "POINT ")
 163                     ; // silently ignore POINT type since we don't need them at the moment
 164                 else
 165                     std::cerr << "Warning: contour type " << contour_type << " not handled!" << std::endl;
 166             }
 167             append->Update();
 168             current_roi->AddMesh(append->GetOutput());
 169             result.push_back(current_roi);
 170         }
 171         else
 172         {
 173             //std::cerr << "Warning: ignoring ROI #" << roi_number << std::endl;
 174         }
 175         k=roi_info->GetNextSQItem(); //increment the second loop variable
 176     }
 177     return result;
 178 }
 179