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 ===========================================================================**/
20 #include "clitkDicom2Image_ggo.h"
21 #include "clitkCommon.h"
22 #include "clitkImageCommon.h"
23 #include "vvImageReader.h"
24 #include "vvImageWriter.h"
26 #include <vtkImageChangeInformation.h>
27 #if GDCM_MAJOR_VERSION == 2
28 #include <gdcmImageHelper.h>
29 #include <gdcmAttribute.h>
30 #include <gdcmReader.h>
33 //====================================================================
34 int main(int argc, char * argv[])
37 GGO(clitkDicom2Image, args_info);
38 std::vector<std::string> input_files;
39 ///if std_input is given, read the input files from stdin
40 if (args_info.std_input_given) {
45 input_files.push_back(tmp);
48 args_info.inputs_num=input_files.size();
49 } else for (unsigned int i=0; i<args_info.inputs_num; i++)
50 input_files.push_back(args_info.inputs[i]);
53 //===========================================
54 /// Get slices locations ...
55 std::vector<double> theorigin(3);
56 std::vector<double> sliceLocations;
57 for(unsigned int i=0; i<args_info.inputs_num; i++) {
58 //std::cout << "Reading <" << input_files[i] << std::endl;
59 #if GDCM_MAJOR_VERSION == 2
60 if (args_info.verbose_flag)
61 std::cout << "Using GDCM-2.x" << std::endl;
63 hreader.SetFileName(input_files[i].c_str());
65 theorigin = gdcm::ImageHelper::GetOriginValue(hreader.GetFile());
66 sliceLocations.push_back(theorigin[2]);
67 gdcm::Attribute<0x28, 0x100> pixel_size;
68 gdcm::DataSet& ds = hreader.GetFile().GetDataSet();
69 pixel_size.SetFromDataSet(ds);
70 if (pixel_size.GetValue() != 16)
72 std::cerr << "Pixel type 2 bytes ! " << std::endl;
73 std::cerr << "In file " << input_files[i] << std::endl;
77 if (args_info.verbose_flag)
78 std::cout << "Not using GDCM-2.x" << std::endl;
79 gdcm::File *header = new gdcm::File();
80 header->SetFileName(input_files[i]);
81 header->SetMaxSizeLoadEntry(16384); // required ?
83 theorigin[0] = header->GetXOrigin();
84 theorigin[1] = header->GetYOrigin();
85 theorigin[2] = header->GetZOrigin();
86 sliceLocations.push_back(theorigin[2]);
87 if (header->GetPixelSize() != 2) {
88 std::cerr << "Pixel type 2 bytes ! " << std::endl;
89 std::cerr << "In file " << input_files[i] << std::endl;
95 //===========================================
96 // Sort slices locations ...
97 std::vector<int> sliceIndex;
98 clitk::GetSortedIndex(sliceLocations, sliceIndex);
99 if (args_info.verboseSliceLocation_flag) {
100 std::cout << sliceLocations[sliceIndex[0]] << " -> "
101 << sliceIndex[0] << " / " << 0 << " => "
103 << input_files[sliceIndex[0]]
105 for(unsigned int i=1; i<sliceIndex.size(); i++) {
106 std::cout << sliceLocations[sliceIndex[i]] << " -> "
107 << sliceIndex[i] << " / " << i << " => "
108 << sliceLocations[sliceIndex[i]] - sliceLocations[sliceIndex[i-1]]
110 << input_files[sliceIndex[i]]
115 //===========================================
116 // Analyze slices locations ...
119 double tolerance = args_info.tolerance_arg;
120 double previous = sliceLocations[sliceIndex[0]];
121 for(unsigned int i=1; i<sliceIndex.size(); i++) {
122 currentDist = sliceLocations[sliceIndex[i]]-previous;
124 if (fabs(dist-currentDist) > tolerance) {
125 std::cout << "ERROR : " << std::endl
126 << "Current slice pos is = " << sliceLocations[sliceIndex[i]] << std::endl
127 << "Previous slice pos is = " << previous << std::endl
128 << "Current file is = " << input_files[sliceIndex[i]] << std::endl
129 << "Current index is = " << i << std::endl
130 << "Current sortindex is = " << sliceIndex[i] << std::endl
131 << "Current slice diff = " << dist << std::endl
132 << "Current error = " << fabs(dist-currentDist) << std::endl;
135 } else dist = currentDist;
136 previous = sliceLocations[sliceIndex[i]];
139 //===========================================
140 // Create ordered vector of filenames
141 std::vector<std::string> sorted_files;
142 sorted_files.resize(sliceIndex.size());
143 for(unsigned int i=0; i<sliceIndex.size(); i++)
144 sorted_files[i] = input_files[ sliceIndex[i] ];
146 //===========================================
148 vvImageReader::Pointer reader = vvImageReader::New();
149 reader->SetInputFilenames(sorted_files);
150 reader->Update(vvImageReader::DICOM);
151 if (reader->GetLastError().size() != 0) {
152 std::cerr << reader->GetLastError() << std::endl;
156 vvImage::Pointer image = reader->GetOutput();
157 vtkImageData* vtk_image = image->GetFirstVTKImageData();
158 vtkImageChangeInformation* modifier = vtkImageChangeInformation::New();
159 if (args_info.focal_origin_given) {
160 std::vector<double> spacing = image->GetSpacing();
161 std::vector<int> size = image->GetSize();
162 theorigin[0] = -spacing[0]*size[0]/2.0;
163 theorigin[1] = -spacing[1]*size[1]/2.0;
164 modifier->SetInput(vtk_image);
165 modifier->SetOutputOrigin(theorigin[0], theorigin[1], sliceLocations[sliceIndex[0]]);
167 vvImage::Pointer focal_image = vvImage::New();
168 focal_image->AddVtkImage(modifier->GetOutput());
172 vvImageWriter::Pointer writer = vvImageWriter::New();
173 writer->SetInput(image);
174 writer->SetOutputFileName(args_info.output_arg);