vv/vvMeshReader.cxx

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