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