]> Creatis software - clitk.git/blob - vv/vvMeshReader.cxx
Merge branch 'master' of git.creatis.insa-lyon.fr:clitk
[clitk.git] / 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://www.centreleonberard.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 <gdcmFile.h>
23 #if GDCM_MAJOR_VERSION == 2
24   #include <gdcmReader.h>
25   #include <gdcmTag.h>
26   #include <gdcmAttribute.h>
27 #else
28   #include <gdcm.h>
29   #include <gdcmSQItem.h>
30 #endif
31
32 #include <vtkSmartPointer.h>
33 #include <vtkAppendPolyData.h>
34 #include <vtkCellArray.h>
35 #include <vtkMetaImageWriter.h>
36 #include <vtkMetaImageReader.h>
37 #include <vtkPolyDataWriter.h>
38 #include <vtkImageData.h>
39 #include <vtkImageCast.h>
40 #include <vtkImageGaussianSmooth.h>
41
42 #include "clitkCommon.h"
43 #include "vvMeshReader.h"
44 #include "vvProgressDialog.h"
45 #include <clitkCommon.h>
46
47 vvMeshReader::vvMeshReader() :
48   vtk_mode(false)
49 {}
50
51 void vvMeshReader::Update()
52 {
53   //Show a progress bar only when opening a DC-struct (ie. multiple contours)
54   vvProgressDialog progress("Opening " + filename,(!vtk_mode) && (selected_contours.size()>1));
55   this->start();
56   while (this->isRunning()) {
57     progress.SetProgress(result.size(),selected_contours.size());
58     this->wait(50);
59     qApp->processEvents();
60   }
61 }
62
63 void vvMeshReader::run()
64 {
65   ///Verify the important stuff has been set
66   assert(filename != "");
67   if (vtk_mode) { //Read vtkPolyData
68     vvMesh::Pointer m=vvMesh::New();
69     m->ReadFromVTK(filename.c_str());
70     if (vf) m->propagateContour(vf);
71     m->ComputeMasks(image->GetVTKImages()[0],false); //don't extrude the contour
72     result.push_back(m);
73   } else { //Read a Dicom-struct file
74     assert(selected_contours.size() > 0);
75     assert(image);
76     std::vector<vvMesh::Pointer> contour_stacks=readSelectedContours();
77     for (std::vector<vvMesh::Pointer>::iterator i=contour_stacks.begin();
78          i!=contour_stacks.end(); i++) {
79       (*i)->ComputeMasks(image->GetVTKImages()[0],true); //Remesh the contour
80       (*i)->ComputeMeshes();
81       if (vf) (*i)->propagateContour(vf);
82       result.push_back(*i);
83     }
84   }
85 }
86
87 std::vector<std::pair<int,std::string> > vvMeshReader::GetROINames()
88 {
89   assert(filename!="");
90   std::vector<std::pair<int, std::string> > roi_names;
91 #if GDCM_MAJOR_VERSION == 2
92   // duplicate code from  clitk::DicomRT_StructureSet::Read
93   gdcm::Reader reader;
94   reader.SetFileName( filename.c_str() );
95   reader.Read();
96
97   const gdcm::DataSet &ds = reader.GetFile().GetDataSet();
98
99   gdcm::Tag tssroisq(0x3006,0x0020);
100   const gdcm::DataElement &ssroisq = ds.GetDataElement( tssroisq );
101   gdcm::SmartPointer<gdcm::SequenceOfItems> roi_seq = ssroisq.GetValueAsSQ();
102   assert(roi_seq); // TODO error message
103   for(unsigned int ridx = 0; ridx < roi_seq->GetNumberOfItems(); ++ridx)
104     {
105     gdcm::Item & item = roi_seq->GetItem( ridx + 1); // Item starts at 1
106     const gdcm::DataSet& nestedds = item.GetNestedDataSet();
107     if( nestedds.FindDataElement( gdcm::Tag(0x3006,0x22) ) )
108       {
109       gdcm::Attribute<0x3006,0x26> roiname;
110       roiname.SetFromDataSet( nestedds );
111       std::string name = roiname.GetValue();      // 0x3006,0x0026 = [ROI Name]
112       gdcm::Attribute<0x3006,0x0022> roinumber;
113       roinumber.SetFromDataSet( nestedds );
114       int nb = roinumber.GetValue();  // 0x3006,0x0022 = [ROI Number]
115
116       roi_names.push_back(make_pair(nb,name));
117       }
118     }
119 #else
120   gdcm::File reader;
121   reader.SetFileName(filename.c_str());
122   reader.SetMaxSizeLoadEntry(16384);
123   reader.Load();
124
125   gdcm::SeqEntry * roi_info=reader.GetSeqEntry(0x3006,0x0020);
126   assert(roi_info);
127   // DD("ici");
128   //int n=0;
129   for (gdcm::SQItem* i=roi_info->GetFirstSQItem(); i!=0; i=roi_info->GetNextSQItem())
130     if (i->GetEntryValue(0x3006,0x0022)!= gdcm::GDCM_UNFOUND)
131       roi_names.push_back(make_pair(atoi(i->GetEntryValue(0x3006,0x0022).c_str()),i->GetEntryValue(0x3006,0x0026)));
132 #endif
133   return roi_names;
134 }
135
136 std::vector<vvMesh::Pointer> vvMeshReader::readSelectedContours()
137 {
138   std::vector<vvMesh::Pointer> result;
139 #if GDCM_MAJOR_VERSION == 2
140   gdcm::Reader reader;
141   reader.SetFileName(filename.c_str());
142   reader.Read();
143
144   const gdcm::DataSet &ds = reader.GetFile().GetDataSet();
145
146   gdcm::SmartPointer<gdcm::SequenceOfItems> rois = ds.GetDataElement(gdcm::Tag(0x3006,0x39)).GetValueAsSQ();
147   gdcm::SmartPointer<gdcm::SequenceOfItems> roi_info = ds.GetDataElement(gdcm::Tag(0x3006,0x20)).GetValueAsSQ();
148   assert(rois); // TODO error message
149   assert(roi_info); // TODO error message
150   assert(rois->GetNumberOfItems() == roi_info->GetNumberOfItems());
151
152   for (unsigned ridx = 0; ridx < rois->GetNumberOfItems(); ++ridx)
153   {
154     vtkSmartPointer<vtkAppendPolyData> append=vtkSmartPointer<vtkAppendPolyData>::New();
155     const gdcm::DataSet& ds_rois = rois->GetItem( ridx + 1).GetNestedDataSet();
156     const gdcm::DataSet& ds_roi_info = roi_info->GetItem( ridx + 1).GetNestedDataSet();
157
158     gdcm::Attribute<0x3006,0x84> roinumber;
159     roinumber.SetFromDataSet(ds_rois);
160     if (std::find(selected_contours.begin(), selected_contours.end(), roinumber.GetValue()) != selected_contours.end()) //Only read selected ROIs
161     {
162       gdcm::Attribute<0x3006,0x2a> trgb;
163       trgb.SetFromDataSet(ds_rois);
164       vvMesh::Pointer current_roi=vvMesh::New();
165       current_roi->r = trgb[0] / 255;
166       current_roi->g = trgb[1] / 255;
167       current_roi->b = trgb[2] / 255;
168
169       gdcm::Attribute<0x3006,0x26> tstructure_name;
170       tstructure_name.SetFromDataSet(ds_roi_info);
171       current_roi->structure_name = tstructure_name.GetValue();
172
173       gdcm::SmartPointer<gdcm::SequenceOfItems> roi_seq = ds_rois.GetDataElement(gdcm::Tag(0x3006,0x40)).GetValueAsSQ();
174       double z0=-1; //Used to determine spacing between slices, assumed to be constant
175       for (unsigned j = 0; j < roi_seq->GetNumberOfItems(); ++j)
176       {
177         gdcm::Item & item_roi_seq = roi_seq->GetItem(j + 1); // Item starts at 1
178         const gdcm::DataSet& ds_roi_seq = item_roi_seq.GetNestedDataSet();
179         gdcm::Attribute<0x3006,0x42> tcontour_type;
180         tcontour_type.SetFromDataSet(ds_roi_seq);
181         std::string contour_type = tcontour_type.GetValue();
182         if (contour_type=="CLOSED_PLANAR ")
183         {
184           gdcm::Attribute<0x3006,0x46> tpoint_number;
185           tpoint_number.SetFromDataSet(ds_roi_seq);
186           const gdcm::DataElement & points_data = ds_roi_seq.GetDataElement(gdcm::Tag(0x3006,0x50));
187           gdcm::Attribute<0x3006,0x50> tpoints;
188           tpoints.SetFromDataElement(points_data);
189           assert(tpoints.GetNumberOfValues() == static_cast<unsigned int>(tpoint_number.GetValue()) * 3);
190           const double* points = tpoints.GetValues();
191           if (z0 == -1) //First contour
192             z0=points[2];
193           else
194             if (current_roi->GetSpacing()==-1 && points[2] != z0 )
195               current_roi->SetSpacing(points[2]-z0);
196           vtkPolyData * contour=vtkPolyData::New();
197           contour->Allocate(); //for cell structures
198           contour->SetPoints(vtkPoints::New());
199           vtkIdType ids[2];
200           for (unsigned idx = 0; idx < tpoints.GetNumberOfValues(); idx += 3)
201           {
202             contour->GetPoints()->InsertNextPoint(points[idx], points[idx+1], points[idx+2]);
203             ids[0] = idx / 3;
204             ids[1] = (ids[0] + 1) % tpoint_number.GetValue(); //0-1,1-2,...,n-1-0
205             contour->GetLines()->InsertNextCell(2, ids);
206           }
207           append->AddInput(contour);
208         }
209         else
210           if (contour_type == "POINT ")
211           ; // silently ignore POINT type since we don't need them at the moment
212           else
213             std::cerr << "Warning: contour type " << contour_type << " not handled!" << std::endl;
214       }
215       append->Update();
216       current_roi->AddMesh(append->GetOutput());
217       result.push_back(current_roi);
218     }
219     else
220     {
221     //std::cerr << "Warning: ignoring ROI #" << roi_number << std::endl;
222     }
223   }
224 #else
225   gdcm::File reader;
226   reader.SetFileName(filename.c_str());
227   reader.SetMaxSizeLoadEntry(16384);
228   reader.Load();
229
230   gdcm::SeqEntry * rois=reader.GetSeqEntry(0x3006,0x0039);
231   ///We need to iterate both on the contours themselves, and on the contour info
232   gdcm::SeqEntry * roi_info=reader.GetSeqEntry(0x3006,0x0020);
233   gdcm::SQItem* k=roi_info->GetFirstSQItem();
234   for(gdcm::SQItem* i=rois->GetFirstSQItem(); i!=0; i=rois->GetNextSQItem()) { //loop over ROIS
235     assert(k!=0);
236     vtkSmartPointer<vtkAppendPolyData> append=vtkSmartPointer<vtkAppendPolyData>::New();
237     std::istringstream ss(i->GetEntryValue(0x3006,0x0084));
238     int roi_number;
239     ss >> roi_number;
240     if (std::find(selected_contours.begin(),selected_contours.end(),roi_number) != selected_contours.end()) { //Only read selected ROIs
241       vvMesh::Pointer current_roi=vvMesh::New();
242       std::vector<double> rgb=clitk::parse_string<double>(i->GetEntryValue(0x3006,0x002a),'\\');
243       assert(rgb.size()==3);
244       current_roi->r=rgb[0]/255;
245       current_roi->g=rgb[1]/255;
246       current_roi->b=rgb[2]/255;
247       current_roi->structure_name=k->GetEntryValue(0x3006,0x0026);
248       gdcm::SeqEntry * contours=i->GetSeqEntry(0x3006,0x0040);
249       double z0=-1; //Used to determine spacing between slices, assumed to be constant
250       for(gdcm::SQItem* j=contours->GetFirstSQItem(); j!=0; j=contours->GetNextSQItem()) { //loop over 2D contours
251         std::string contour_type=j->GetEntryValue(0x3006,0x0042);
252         if (contour_type=="CLOSED_PLANAR ") {
253           int point_number=clitk::parse_value<int>(j->GetEntryValue(0x3006,0x0046));
254           std::vector<float> points=clitk::parse_string<float>(j->GetEntryValue(0x3006,0x0050),'\\');
255           assert(points.size() == static_cast<unsigned int>(point_number)*3);
256           if (z0 == -1) //First contour
257             z0=points[2];
258           //2nd contour, spacing not yet set. Need to be sure we are on a different slice,
259           //sometimes there is more than one closed contour per slice
260           else if (current_roi->GetSpacing()==-1 && points[2] != z0 )
261             current_roi->SetSpacing(points[2]-z0);
262           vtkPolyData * contour=vtkPolyData::New();
263           contour->Allocate(); //for cell structures
264           contour->SetPoints(vtkPoints::New());
265           vtkIdType ids[2];
266           for (unsigned int idx=0; idx<points.size(); idx+=3) {
267             contour->GetPoints()->InsertNextPoint(points[idx],points[idx+1],points[idx+2]);
268             ids[0]=idx/3;
269             ids[1]=(ids[0]+1)%point_number; //0-1,1-2,...,n-1-0
270             contour->GetLines()->InsertNextCell(2,ids);
271           }
272           append->AddInput(contour);
273         } else if (contour_type == "POINT ")
274           ; // silently ignore POINT type since we don't need them at the moment
275         else
276           std::cerr << "Warning: contour type " << contour_type << " not handled!" << std::endl;
277       }
278       append->Update();
279       current_roi->AddMesh(append->GetOutput());
280       result.push_back(current_roi);
281     } else {
282       //std::cerr << "Warning: ignoring ROI #" << roi_number << std::endl;
283     }
284     k=roi_info->GetNextSQItem(); //increment the second loop variable
285   }
286 #endif
287   return result;
288 }
289