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 ===========================================================================**/
18 #ifndef CLITKSCINTIVOLSTATS_CXX
19 #define CLITKSCINTIVOLSTATS_CXX
22 #include "clitkScintivolStats_ggo.h"
24 #include "clitkImageCommon.h"
25 #include "clitkCommon.h"
27 #include "itkImageFileReader.h"
28 #include "itkBinaryImageToLabelMapFilter.h"
29 #include "itkLabelMapToLabelImageFilter.h"
30 #include "itkLabelStatisticsImageFilter.h"
31 #include "itkStatisticsImageFilter.h"
32 #include "itkExtractImageFilter.h"
36 //-------------------------------------------------------------------=
37 int main(int argc, char * argv[])
41 GGO(clitkScintivolStats, args_info);
44 typedef itk::Image<double, 3> Input3DType;
45 typedef itk::Image<double, 2> Input2DType;
46 typedef itk::Image<char, 3> Input3DMaskType;
47 typedef itk::Image<char, 2> Input2DMaskType;
48 typedef itk::ImageFileReader<Input3DType> InputReader3DType;
49 typedef itk::ImageFileReader<Input3DMaskType> InputReader3DMaskType;
50 typedef itk::ImageFileReader<Input2DMaskType> InputReader2DMaskType;
51 typedef itk::LabelStatisticsImageFilter<Input3DType, Input3DMaskType> LabelStatistics3DImageFilterType;
52 typedef itk::LabelStatisticsImageFilter<Input2DType, Input2DMaskType> LabelStatistics2DImageFilterType;
53 typedef itk::StatisticsImageFilter<Input2DType> Statistics2DImageFilterType;
54 typedef itk::ExtractImageFilter<Input3DType, Input2DType> ExtractImageFilter;
56 //Determine nbAcquisitionDynamic1 and nbAcquisitionDynamic2 (ie. the number of slice in dynamic1 and dynamic2)
57 //Start to open the dynamic1 and dynamic2 and get the 3rd dimension size
58 InputReader3DType::Pointer readerDynamic1 = InputReader3DType::New();
59 readerDynamic1->SetFileName(args_info.dynamic1_arg);
60 readerDynamic1->Update();
61 Input3DType::Pointer dynamic1 = readerDynamic1->GetOutput();
62 int nbAcquisitionDynamic1 = dynamic1->GetLargestPossibleRegion().GetSize()[2];
64 InputReader3DType::Pointer readerDynamic2 = InputReader3DType::New();
65 readerDynamic2->SetFileName(args_info.dynamic2_arg);
66 readerDynamic2->Update();
67 Input3DType::Pointer dynamic2 = readerDynamic2->GetOutput();
68 int nbAcquisitionDynamic2 = dynamic2->GetLargestPossibleRegion().GetSize()[2];
71 //If it's scatter or attenuation correction, just append results to the csv file
72 //If not open it normally and write the number of dynamic acquisition for the first and the second images (given by the ggo and from the dicom tag)
73 std::ofstream csvFile;
74 if (args_info.append_flag)
75 csvFile.open (args_info.output_arg, std::ios::app);
77 csvFile.open (args_info.output_arg);
78 csvFile << "1;" << nbAcquisitionDynamic1 << ";" << nbAcquisitionDynamic2 << "\n";
81 //Read Tomo image, total and remnant masks
82 InputReader3DType::Pointer readerTomo = InputReader3DType::New();
83 readerTomo->SetFileName(args_info.tomo_arg);
85 Input3DType::Pointer tomo = readerTomo->GetOutput();
87 InputReader3DMaskType::Pointer readerTotalLiver = InputReader3DMaskType::New();
88 readerTotalLiver->SetFileName(args_info.totalLiverMask_arg);
89 readerTotalLiver->Update();
91 InputReader3DMaskType::Pointer readerRemnantLiver = InputReader3DMaskType::New();
92 readerRemnantLiver->SetFileName(args_info.remnantLiverMask_arg);
93 readerRemnantLiver->Update();
95 //Find number of counts in tomo for total Liver and remnant Liver
96 LabelStatistics3DImageFilterType::Pointer labelStatisticsImageFilterTotalLiver = LabelStatistics3DImageFilterType::New();
97 labelStatisticsImageFilterTotalLiver->SetLabelInput(readerTotalLiver->GetOutput());
98 labelStatisticsImageFilterTotalLiver->SetInput(tomo);
99 labelStatisticsImageFilterTotalLiver->Update();
101 LabelStatistics3DImageFilterType::Pointer labelStatisticsImageFilterRemnantLiver = LabelStatistics3DImageFilterType::New();
102 labelStatisticsImageFilterRemnantLiver->SetLabelInput(readerRemnantLiver->GetOutput());
103 labelStatisticsImageFilterRemnantLiver->SetInput(tomo);
104 labelStatisticsImageFilterRemnantLiver->Update();
106 //Write them in the csv file
107 csvFile << args_info.acquisitionTimeTomo_arg << ";" << labelStatisticsImageFilterTotalLiver->GetSum(1) << ";" << labelStatisticsImageFilterRemnantLiver->GetSum(1) << "\n";
109 //Read liver and heart masks for dynamic1
110 InputReader2DMaskType::Pointer readerLiver = InputReader2DMaskType::New();
111 readerLiver->SetFileName(args_info.liverMask_arg);
112 readerLiver->Update();
114 InputReader2DMaskType::Pointer readerHeart = InputReader2DMaskType::New();
115 readerHeart->SetFileName(args_info.heartMask_arg);
116 readerHeart->Update();
118 for (unsigned int i=0; i<nbAcquisitionDynamic1; ++i) {
119 //Compute time frame in s
120 double timeFrame = args_info.frameDurationDynamic1_arg * i + args_info.frameDurationDynamic1_arg;
122 //Extract the Region of interest of the dynamic1
123 ExtractImageFilter::Pointer extractSlice = ExtractImageFilter::New();
124 Input3DType::IndexType start;
128 Input3DType::SizeType size;
129 size[0] = dynamic1->GetLargestPossibleRegion().GetSize()[0];
130 size[1] = dynamic1->GetLargestPossibleRegion().GetSize()[1];
132 Input3DType::RegionType desiredRegion;
133 desiredRegion.SetSize(size);
134 desiredRegion.SetIndex(start);
135 extractSlice->SetExtractionRegion(desiredRegion);
136 extractSlice->SetInput(dynamic1);
137 #if ITK_VERSION_MAJOR >= 4
138 extractSlice->SetDirectionCollapseToIdentity();
140 extractSlice->Update();
142 //Find number of counts in dynamic1 slice for Liver and heart
143 LabelStatistics2DImageFilterType::Pointer labelStatisticsImageFilterLiver = LabelStatistics2DImageFilterType::New();
144 labelStatisticsImageFilterLiver->SetLabelInput(readerLiver->GetOutput());
145 labelStatisticsImageFilterLiver->SetInput(extractSlice->GetOutput());
146 labelStatisticsImageFilterLiver->SetCoordinateTolerance(0.001);
147 labelStatisticsImageFilterLiver->SetDirectionTolerance(0.001);
148 labelStatisticsImageFilterLiver->Update();
150 LabelStatistics2DImageFilterType::Pointer labelStatisticsImageFilterHeart = LabelStatistics2DImageFilterType::New();
151 labelStatisticsImageFilterHeart->SetLabelInput(readerHeart->GetOutput());
152 labelStatisticsImageFilterHeart->SetInput(extractSlice->GetOutput());
153 labelStatisticsImageFilterHeart->SetCoordinateTolerance(0.001);
154 labelStatisticsImageFilterHeart->SetDirectionTolerance(0.001);
155 labelStatisticsImageFilterHeart->Update();
157 //Find number of counts in dynamic1 slice
158 Statistics2DImageFilterType::Pointer statisticsImageFilterDynamic1 = Statistics2DImageFilterType::New();
159 statisticsImageFilterDynamic1->SetInput(extractSlice->GetOutput());
160 statisticsImageFilterDynamic1->Update();
162 //Write them in the csv file
163 csvFile << timeFrame << ";" << labelStatisticsImageFilterLiver->GetSum(1) << ";" << labelStatisticsImageFilterHeart->GetSum(1) << ";" << statisticsImageFilterDynamic1->GetSum() << "\n";
166 //Read parenchyma mask for dynamic2
167 InputReader2DMaskType::Pointer readerParenchyma = InputReader2DMaskType::New();
168 readerParenchyma->SetFileName(args_info.parenchymaMask_arg);
169 readerParenchyma->Update();
171 for (unsigned int i=0; i<nbAcquisitionDynamic2; ++i) {
172 //Compute time frame in s
173 double timeFrame = args_info.frameDurationDynamic2_arg * i + args_info.acquisitionTimeDynamic2_arg + args_info.frameDurationDynamic2_arg;
175 //Extract the Region of interest of the dynamic1
176 ExtractImageFilter::Pointer extractSlice = ExtractImageFilter::New();
177 Input3DType::IndexType start;
181 Input3DType::SizeType size;
182 size[0] = dynamic2->GetLargestPossibleRegion().GetSize()[0];
183 size[1] = dynamic2->GetLargestPossibleRegion().GetSize()[1];
185 Input3DType::RegionType desiredRegion;
186 desiredRegion.SetSize(size);
187 desiredRegion.SetIndex(start);
188 extractSlice->SetExtractionRegion(desiredRegion);
189 extractSlice->SetInput(dynamic2);
190 #if ITK_VERSION_MAJOR >= 4
191 extractSlice->SetDirectionCollapseToIdentity();
193 extractSlice->Update();
195 //Find number of counts in dynamic2 slice for parenchyma
196 LabelStatistics2DImageFilterType::Pointer labelStatisticsImageFilterParenchyma = LabelStatistics2DImageFilterType::New();
197 labelStatisticsImageFilterParenchyma->SetLabelInput(readerParenchyma->GetOutput());
198 labelStatisticsImageFilterParenchyma->SetInput(extractSlice->GetOutput());
199 labelStatisticsImageFilterParenchyma->SetCoordinateTolerance(0.001);
200 labelStatisticsImageFilterParenchyma->SetDirectionTolerance(0.001);
201 labelStatisticsImageFilterParenchyma->Update();
203 //Write them in the csv file
204 csvFile << timeFrame << ";" << labelStatisticsImageFilterParenchyma->GetSum(1) << "\n";
209 // this is the end my friend
212 //-------------------------------------------------------------------=
214 #endif /* end #define CLITKSCINTIVOLSTATS_CXX */