1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
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
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.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
23 #include <QApplication>
27 #if GDCM_MAJOR_VERSION == 2
28 #include <gdcmReader.h>
30 #include <gdcmAttribute.h>
33 #include <gdcmSQItem.h>
37 #include <vtkSmartPointer.h>
38 #include <vtkAppendPolyData.h>
39 #include <vtkCellArray.h>
40 #include <vtkMetaImageWriter.h>
41 #include <vtkMetaImageReader.h>
42 #include <vtkPolyDataWriter.h>
43 #include <vtkImageData.h>
44 #include <vtkImageCast.h>
45 #include <vtkImageGaussianSmooth.h>
48 #include "clitkCommon.h"
49 #include "clitkCommon.h"
50 #include "clitkDicomRT_StructureSet.h"
51 #include "vvMeshReader.h"
52 #include "vvProgressDialog.h"
54 //------------------------------------------------------------------------------
55 vvMeshReader::vvMeshReader() :
58 //------------------------------------------------------------------------------
61 //------------------------------------------------------------------------------
62 void vvMeshReader::Update()
64 DD("vvMeshReader::Update");
66 //Show a progress bar only when opening a DC-struct (ie. multiple contours)
67 vvProgressDialog progress("Opening " + filename,(!vtk_mode) && (selected_contours.size()>1));
69 while (this->isRunning()) {
70 progress.SetProgress(result.size(),selected_contours.size());
72 qApp->processEvents();
75 //------------------------------------------------------------------------------
78 //------------------------------------------------------------------------------
79 void vvMeshReader::run()
81 ///Verify the important stuff has been set
82 assert(filename != "");
83 if (vtk_mode) { //Read vtkPolyData
84 vvMesh::Pointer m=vvMesh::New();
85 m->ReadFromVTK(filename.c_str());
86 if (vf) m->propagateContour(vf);
87 m->ComputeMasks(image->GetVTKImages()[0],false); //don't extrude the contour
89 } else { //Read a Dicom-struct file
90 assert(selected_contours.size() > 0);
92 std::vector<vvMesh::Pointer> contour_stacks=readSelectedContours();
93 for (std::vector<vvMesh::Pointer>::iterator i=contour_stacks.begin();
94 i!=contour_stacks.end(); i++) {
95 (*i)->ComputeMasks(image->GetVTKImages()[0],true); //Remesh the contour
96 (*i)->ComputeMeshes();
97 if (vf) (*i)->propagateContour(vf);
102 //------------------------------------------------------------------------------
105 //------------------------------------------------------------------------------
106 std::vector<std::pair<int,std::string> > vvMeshReader::GetROINames()
108 assert(filename!="");
109 std::vector<std::pair<int, std::string> > roi_names;
110 DD(GDCM_MAJOR_VERSION);
111 DD(CLITK_USE_SYSTEM_GDCM);
113 #if GDCM_MAJOR_VERSION == 2
115 // Read RT-struct data
117 vtkSmartPointer<vtkGDCMPolyDataReader> reader = vtkGDCMPolyDataReader::New();
118 reader->SetFileName(filename.c_str());
122 // get info on roi names
123 vtkRTStructSetProperties * p = reader->GetRTStructSetProperties();
124 DD(p->GetNumberOfStructureSetROIs());
125 DD(p->GetStructureSetROIName(0));
126 DD(p->GetStructureSetROINumber(0));
127 int n = p->GetNumberOfStructureSetROIs();
130 for(unsigned int i=0; i<n; i++) {
131 std::string name = p->GetStructureSetROIName(i);
132 int nb = p->GetStructureSetROINumber(i);
133 roi_names.push_back(make_pair(nb,name));
136 // ====================
139 // duplicate code from clitk::DicomRT_StructureSet::Read
140 gdcm::Reader * reader = new gdcm::Reader;
141 reader->SetFileName( filename.c_str() );
143 DD("after gdcm read");
145 const gdcm::DataSet &ds = reader->GetFile().GetDataSet();
151 //Verify if the file is a RT-Structure-Set dicom file
152 gdcm::File * mFile = &(reader->GetFile());
153 gdcm::MediaStorage ms;
154 ms.SetFromFile(*mFile);
155 if( ms != gdcm::MediaStorage::RTStructureSetStorage )
157 std::cerr << "Error. the file " << filename
158 << " is not a Dicom Struct ? (must have a SOP Class UID [0008|0016] = 1.2.840.10008.5.1.4.1.1.481.3 ==> [RT Structure Set Storage])"
163 gdcm::Attribute<0x8,0x60> modality;
164 modality.SetFromDataSet( ds );
165 if( modality.GetValue() != "RTSTRUCT" )
167 std::cerr << "Error. the file " << filename
168 << " is not a Dicom Struct ? (must have 0x0008,0x0060 = RTSTRUCT [RT Structure Set Storage])"
173 gdcm::Attribute<0x20,0x10> studyid;
174 studyid.SetFromDataSet( ds );
175 DD(studyid.GetValue());
177 gdcm::Tag tssroisq(0x3006,0x0020);
178 // 0x3006,0x0020 = [ Structure Set ROI Sequence ]
179 bool b = ds.FindDataElement(tssroisq);
181 clitkExceptionMacro("Error: tag 0x3006,0x0020 [ Structure Set ROI Sequence ] not found");
184 const gdcm::DataElement &ssroisq = ds.GetDataElement( tssroisq );
186 gdcm::SmartPointer<gdcm::SequenceOfItems> roi_seq = ssroisq.GetValueAsSQ();
187 assert(roi_seq); // FIXME error message
189 DD(roi_seq->GetNumberOfItems());
191 for(unsigned int ridx = 0; ridx < roi_seq->GetNumberOfItems(); ++ridx)
194 gdcm::Item & item = roi_seq->GetItem( ridx + 1); // Item starts at 1
196 const gdcm::Item & sitem = roi_seq->GetItem(ridx+1); // Item start at #1
198 const gdcm::DataSet& snestedds = sitem.GetNestedDataSet();
199 const gdcm::DataSet& nestedds = item.GetNestedDataSet();
201 DD(nestedds.IsEmpty());
203 if( snestedds.FindDataElement( gdcm::Tag(0x3006,0x22) ) )
206 // const gdcm::DataElement & a = nestedds.GetDataElement(gdcm::Tag(0x3006,0x26));
209 gdcm::Attribute<0x3006,0x26> roiname;
210 roiname.SetFromDataSet( snestedds );
211 std::string name = roiname.GetValue(); // 0x3006,0x0026 = [ROI Name]
213 gdcm::Attribute<0x3006,0x0022> roinumber;
214 roinumber.SetFromDataSet( snestedds );
215 int nb = roinumber.GetValue(); // 0x3006,0x0022 = [ROI Number]
218 roi_names.push_back(make_pair(nb,name));
228 reader.SetFileName(filename.c_str());
229 reader.SetMaxSizeLoadEntry(16384);
232 gdcm::SeqEntry * roi_info=reader.GetSeqEntry(0x3006,0x0020);
236 for (gdcm::SQItem* i=roi_info->GetFirstSQItem(); i!=0; i=roi_info->GetNextSQItem())
237 if (i->GetEntryValue(0x3006,0x0022)!= gdcm::GDCM_UNFOUND)
238 roi_names.push_back(make_pair(atoi(i->GetEntryValue(0x3006,0x0022).c_str()),i->GetEntryValue(0x3006,0x0026)));
242 //------------------------------------------------------------------------------
245 //------------------------------------------------------------------------------
246 std::vector<vvMesh::Pointer> vvMeshReader::readSelectedContours()
248 std::vector<vvMesh::Pointer> result;
249 #if GDCM_MAJOR_VERSION == 2
251 reader.SetFileName(filename.c_str());
254 const gdcm::DataSet &ds = reader.GetFile().GetDataSet();
256 gdcm::SmartPointer<gdcm::SequenceOfItems> rois = ds.GetDataElement(gdcm::Tag(0x3006,0x39)).GetValueAsSQ();
257 gdcm::SmartPointer<gdcm::SequenceOfItems> roi_info = ds.GetDataElement(gdcm::Tag(0x3006,0x20)).GetValueAsSQ();
258 assert(rois); // TODO error message
259 assert(roi_info); // TODO error message
260 assert(rois->GetNumberOfItems() == roi_info->GetNumberOfItems());
262 for (unsigned ridx = 0; ridx < rois->GetNumberOfItems(); ++ridx)
264 vtkSmartPointer<vtkAppendPolyData> append=vtkSmartPointer<vtkAppendPolyData>::New();
265 const gdcm::DataSet& ds_rois = rois->GetItem( ridx + 1).GetNestedDataSet();
266 const gdcm::DataSet& ds_roi_info = roi_info->GetItem( ridx + 1).GetNestedDataSet();
268 gdcm::Attribute<0x3006,0x84> roinumber;
269 roinumber.SetFromDataSet(ds_rois);
270 if (std::find(selected_contours.begin(), selected_contours.end(), roinumber.GetValue()) != selected_contours.end()) //Only read selected ROIs
272 gdcm::Attribute<0x3006,0x2a> trgb;
273 trgb.SetFromDataSet(ds_rois);
274 vvMesh::Pointer current_roi=vvMesh::New();
275 current_roi->r = trgb[0] / 255.0;
276 current_roi->g = trgb[1] / 255.0;
277 current_roi->b = trgb[2] / 255.0;
279 gdcm::Attribute<0x3006,0x26> tstructure_name;
280 tstructure_name.SetFromDataSet(ds_roi_info);
281 current_roi->structure_name = tstructure_name.GetValue();
283 gdcm::SmartPointer<gdcm::SequenceOfItems> roi_seq = ds_rois.GetDataElement(gdcm::Tag(0x3006,0x40)).GetValueAsSQ();
284 double z0=-1; //Used to determine spacing between slices, assumed to be constant
285 for (unsigned j = 0; j < roi_seq->GetNumberOfItems(); ++j)
287 gdcm::Item & item_roi_seq = roi_seq->GetItem(j + 1); // Item starts at 1
288 const gdcm::DataSet& ds_roi_seq = item_roi_seq.GetNestedDataSet();
289 gdcm::Attribute<0x3006,0x42> tcontour_type;
290 tcontour_type.SetFromDataSet(ds_roi_seq);
291 std::string contour_type = tcontour_type.GetValue();
292 if (contour_type=="CLOSED_PLANAR ")
294 gdcm::Attribute<0x3006,0x46> tpoint_number;
295 tpoint_number.SetFromDataSet(ds_roi_seq);
296 const gdcm::DataElement & points_data = ds_roi_seq.GetDataElement(gdcm::Tag(0x3006,0x50));
297 gdcm::Attribute<0x3006,0x50> tpoints;
298 tpoints.SetFromDataElement(points_data);
299 assert(tpoints.GetNumberOfValues() == static_cast<unsigned int>(tpoint_number.GetValue()) * 3);
300 const double* points = tpoints.GetValues();
301 if (z0 == -1) //First contour
304 if (current_roi->GetSpacing()==-1 && points[2] != z0 )
305 current_roi->SetSpacing(points[2]-z0);
306 vtkPolyData * contour=vtkPolyData::New();
307 contour->Allocate(); //for cell structures
308 contour->SetPoints(vtkPoints::New());
310 for (unsigned idx = 0; idx < tpoints.GetNumberOfValues(); idx += 3)
312 contour->GetPoints()->InsertNextPoint(points[idx], points[idx+1], points[idx+2]);
314 ids[1] = (ids[0] + 1) % tpoint_number.GetValue(); //0-1,1-2,...,n-1-0
315 contour->GetLines()->InsertNextCell(2, ids);
317 append->AddInput(contour);
320 if (contour_type == "POINT ")
321 ; // silently ignore POINT type since we don't need them at the moment
323 std::cerr << "Warning: contour type " << contour_type << " not handled!" << std::endl;
326 current_roi->AddMesh(append->GetOutput());
327 result.push_back(current_roi);
331 //std::cerr << "Warning: ignoring ROI #" << roi_number << std::endl;
336 reader.SetFileName(filename.c_str());
337 reader.SetMaxSizeLoadEntry(16384);
340 gdcm::SeqEntry * rois=reader.GetSeqEntry(0x3006,0x0039);
341 ///We need to iterate both on the contours themselves, and on the contour info
342 gdcm::SeqEntry * roi_info=reader.GetSeqEntry(0x3006,0x0020);
343 gdcm::SQItem* k=roi_info->GetFirstSQItem();
344 for(gdcm::SQItem* i=rois->GetFirstSQItem(); i!=0; i=rois->GetNextSQItem()) { //loop over ROIS
346 vtkSmartPointer<vtkAppendPolyData> append=vtkSmartPointer<vtkAppendPolyData>::New();
347 std::istringstream ss(i->GetEntryValue(0x3006,0x0084));
350 if (std::find(selected_contours.begin(),selected_contours.end(),roi_number) != selected_contours.end()) { //Only read selected ROIs
351 vvMesh::Pointer current_roi=vvMesh::New();
352 std::vector<double> rgb=clitk::parse_string<double>(i->GetEntryValue(0x3006,0x002a),'\\');
353 assert(rgb.size()==3);
354 current_roi->r=rgb[0]/255;
355 current_roi->g=rgb[1]/255;
356 current_roi->b=rgb[2]/255;
357 current_roi->structure_name=k->GetEntryValue(0x3006,0x0026);
358 gdcm::SeqEntry * contours=i->GetSeqEntry(0x3006,0x0040);
359 double z0=-1; //Used to determine spacing between slices, assumed to be constant
360 for(gdcm::SQItem* j=contours->GetFirstSQItem(); j!=0; j=contours->GetNextSQItem()) { //loop over 2D contours
361 std::string contour_type=j->GetEntryValue(0x3006,0x0042);
362 if (contour_type=="CLOSED_PLANAR ") {
363 int point_number=clitk::parse_value<int>(j->GetEntryValue(0x3006,0x0046));
364 std::vector<float> points=clitk::parse_string<float>(j->GetEntryValue(0x3006,0x0050),'\\');
365 assert(points.size() == static_cast<unsigned int>(point_number)*3);
366 if (z0 == -1) //First contour
368 //2nd contour, spacing not yet set. Need to be sure we are on a different slice,
369 //sometimes there is more than one closed contour per slice
370 else if (current_roi->GetSpacing()==-1 && points[2] != z0 )
371 current_roi->SetSpacing(points[2]-z0);
372 vtkPolyData * contour=vtkPolyData::New();
373 contour->Allocate(); //for cell structures
374 contour->SetPoints(vtkPoints::New());
376 for (unsigned int idx=0; idx<points.size(); idx+=3) {
377 contour->GetPoints()->InsertNextPoint(points[idx],points[idx+1],points[idx+2]);
379 ids[1]=(ids[0]+1)%point_number; //0-1,1-2,...,n-1-0
380 contour->GetLines()->InsertNextCell(2,ids);
382 append->AddInput(contour);
383 } else if (contour_type == "POINT ")
384 ; // silently ignore POINT type since we don't need them at the moment
386 std::cerr << "Warning: contour type " << contour_type << " not handled!" << std::endl;
389 current_roi->AddMesh(append->GetOutput());
390 result.push_back(current_roi);
392 //std::cerr << "Warning: ignoring ROI #" << roi_number << std::endl;
394 k=roi_info->GetNextSQItem(); //increment the second loop variable