1 /*=========================================================================
5 Author : Joel Schaerer (joel.schaerer@insa-lyon.fr)
8 Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
9 CREATIS-LRMN http://www.creatis.insa-lyon.fr
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.
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.
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/>.
23 =========================================================================*/
26 #include <QApplication>
29 #include <gdcmSQItem.h>
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>
41 #include "clitkCommon.h"
42 #include "vvMeshReader.h"
43 #include "vvProgressDialog.h"
45 vvMeshReader::vvMeshReader() :
49 void vvMeshReader::Update()
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));
54 while (this->isRunning())
56 progress.SetProgress(result.size(),selected_contours.size());
58 qApp->processEvents();
62 void vvMeshReader::run()
64 ///Verify the important stuff has been set
65 assert(filename != "");
66 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],true);
74 else //Read a Dicom-struct file
76 assert(selected_contours.size() > 0);
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++)
82 (*i)->ComputeMasks(image->GetVTKImages()[0],true); //Remesh the contour
83 (*i)->ComputeMeshes();
84 if (vf) (*i)->propagateContour(vf);
90 template<class ElementType>
91 ElementType parse_value(std::string str)
93 std::istringstream parser(str);
96 assert(!parser.fail());
100 template<class ElementType>
101 std::vector<ElementType> parse_string(std::string str,char delim)
103 std::istringstream ss(str);
105 std::vector<ElementType> result;
106 while (getline(ss,token,delim))
108 result.push_back(parse_value<ElementType>(token));
113 std::vector<std::pair<int,std::string> > vvMeshReader::GetROINames()
115 assert(filename!="");
117 reader.SetFileName(filename.c_str());
118 reader.SetMaxSizeLoadEntry(16384);
121 gdcm::SeqEntry * roi_info=reader.GetSeqEntry(0x3006,0x0020);
123 std::vector<std::pair<int, std::string> > roi_names;
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)));
132 std::vector<vvMesh::Pointer> vvMeshReader::readSelectedContours()
135 reader.SetFileName(filename.c_str());
136 reader.SetMaxSizeLoadEntry(16384);
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
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
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
161 std::string contour_type=j->GetEntryValue(0x3006,0x0042);
162 if (contour_type=="CLOSED_PLANAR ")
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
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());
177 for (unsigned int idx=0;idx<points.size();idx+=3)
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);
183 append->AddInput(contour);
185 else if (contour_type == "POINT ")
186 ; // silently ignore POINT type since we don't need them at the moment
188 std::cerr << "Warning: contour type " << contour_type << " not handled!" << std::endl;
191 current_roi->AddMesh(append->GetOutput());
192 result.push_back(current_roi);
196 //std::cerr << "Warning: ignoring ROI #" << roi_number << std::endl;
198 k=roi_info->GetNextSQItem(); //increment the second loop variable