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 CLITKIMAGEARITHMGENERICFILTER_CXX
19 #define CLITKIMAGEARITHMGENERICFILTER_CXX
21 #include "clitkImageArithmGenericFilter.h"
26 // class ImageArithmGenericFilter<args_info_clitkImageArithm>;
30 void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<float, 3u>, 3u > >()
32 typedef itk::Image< itk::Vector<float, 3u>, 3u > ImageType;
35 ImageType::Pointer input1 = this->GetInput<ImageType>(0);
37 // Set input image iterator
38 typedef itk::ImageRegionIterator<ImageType> IteratorType;
39 IteratorType it(input1, input1->GetLargestPossibleRegion());
42 ImageType::Pointer input2 = NULL;
46 // Special case for normalisation
47 if (mTypeOfOperation == 12) {
48 typedef itk::MinimumMaximumImageCalculator<ImageType> MinMaxFilterType;
49 typename MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
52 mScalar = ff->GetMaximum();
53 mTypeOfOperation = 11; // divide
57 if (mIsOperationUseASecondImage) {
59 input2 = this->GetInput<ImageType>(1);
60 // Set input image iterator
61 it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
63 if (!clitk::HaveSameSize<ImageType, ImageType>(input1, input2)) {
64 itkExceptionMacro(<< "The images (input and input2) must have the same size");
66 if(!clitk::HaveSameSpacing<ImageType, ImageType>(input1, input2)) {
67 itkWarningMacro(<< "The images (input and input2) do not have the same spacing. "
68 << "Using first input's information.");
73 // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
74 if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) {
75 // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
76 // << typeid(PixelType).name()
78 mOverwriteInputImage = false;
82 // ---------------- Overwrite input Image ---------------------
83 if (mOverwriteInputImage) {
84 // Set output iterator (to input1)
85 IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
86 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
87 else ComputeImage(it, ito);
88 this->SetNextOutput<ImageType>(input1);
90 // ---------------- Create new output Image ---------------------
92 /*if (mOutputIsFloat) {
93 // Create output image
94 typedef itk::Image<float,ImageType::ImageDimension> OutputImageType;
95 typename OutputImageType::Pointer output = OutputImageType::New();
96 output->SetRegions(input1->GetLargestPossibleRegion());
97 output->SetOrigin(input1->GetOrigin());
98 output->SetSpacing(input1->GetSpacing());
100 // Set output iterator
101 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
102 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
103 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
104 else ComputeImage(it, ito);
105 this->template SetNextOutput<OutputImageType>(output);
107 // Create output image
108 typedef ImageType OutputImageType;
109 OutputImageType::Pointer output = OutputImageType::New();
110 output->SetRegions(input1->GetLargestPossibleRegion());
111 output->SetOrigin(input1->GetOrigin());
112 output->SetSpacing(input1->GetSpacing());
114 // Set output iterator
115 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
116 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
117 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
118 else ComputeImage(it, ito);
119 this->SetNextOutput<OutputImageType>(output);
126 void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
127 itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
128 itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > >
129 (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it,
130 itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito)
132 typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter2;
137 typedef Iter2::PixelType PixelType;
139 PixelType scalar_vector;
140 scalar_vector.Fill(mScalar);
143 switch (mTypeOfOperation) {
145 while (!it.IsAtEnd()) {
146 ito.Set(it.Get() + scalar_vector);
152 while (!it.IsAtEnd()) {
153 ito.Set(it.Get() * mScalar);
160 while (!it.IsAtEnd()) {
162 ito.Set(mScalar / it.Get()));
163 else ito.Set(mDefaultPixelValue);
169 while (!it.IsAtEnd()) {
170 if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
171 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
177 while (!it.IsAtEnd()) {
178 if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
179 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
184 case 5: // Absolute value
185 while (!it.IsAtEnd()) {
186 if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
187 // <= zero to avoid warning for unsigned types
188 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
193 case 6: // Squared value
194 while (!it.IsAtEnd()) {
195 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
201 while (!it.IsAtEnd()) {
203 ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
204 else ito.Set(mDefaultPixelValue);
210 while (!it.IsAtEnd()) {
211 ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
217 while (!it.IsAtEnd()) {
219 ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
221 if (it.Get() ==0) ito.Set(0);
222 else ito.Set(mDefaultPixelValue);
229 while (!it.IsAtEnd()) {
230 ito.Set(PixelTypeDownCast<double, PixelType>((0x10000 - (double)it.Get())/mScalar));
237 while (!it.IsAtEnd()) {
238 ito.Set(it.Get() / mScalar);
244 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
252 void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
253 itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
254 itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
255 itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >
257 (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it1,
258 itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it2,
259 itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito)
261 typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter1;
262 typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter2;
263 typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter3;
268 typedef Iter3::PixelType PixelType;
270 switch (mTypeOfOperation) {
272 while (!ito.IsAtEnd()) {
273 ito.Set(it1.Get() + it2.Get());
281 while (!ito.IsAtEnd()) {
282 ito.Set(it1.Get() * it2.Get()) );
289 while (!ito.IsAtEnd()) {
291 ito.Set(it1.Get() / it2.Get()));
292 else ito.Set(mDefaultPixelValue);
299 while (!ito.IsAtEnd()) {
300 if (it1.Get() < it2.Get()) ito.Set(it2.Get());
307 while (!ito.IsAtEnd()) {
308 if (it1.Get() > it2.Get()) ito.Set(it2.Get());
315 case 5: // Absolute difference
316 while (!ito.IsAtEnd()) {
317 ito.Set(it2.Get()-it1.Get());
324 case 6: // Squared differences
325 while (!ito.IsAtEnd()) {
326 ito.Set(pow(it1.Get()-it2.Get(),2)));
333 case 7: // Difference
334 while (!ito.IsAtEnd()) {
335 ito.Set(it1.Get()-it2.Get());
342 case 8: // Relative Difference
343 while (!ito.IsAtEnd()) {
344 if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
353 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
360 #endif //define CLITKIMAGEARITHMGENERICFILTER_CXX