1 /*=========================================================================
3 Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
4 l'Image). All rights reserved. See Doc/License.txt or
5 http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
7 This software is distributed WITHOUT ANY WARRANTY; without even
8 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
9 PURPOSE. See the above copyright notices for more information.
11 =========================================================================*/
13 #ifndef CLITKIMAGEARITHMGENERICFILTER_CXX
14 #define CLITKIMAGEARITHMGENERICFILTER_CXX
17 -------------------------------------------------------------------
18 * @file clitkImageArithmGenericFilter.cxx
19 * @author David Sarrut <David.Sarrut@creatis.insa-lyon.fr>
23 -------------------------------------------------------------------*/
25 #include "clitkImageArithmGenericFilter.h"
27 //--------------------------------------------------------------------
28 clitk::ImageArithmGenericFilter::ImageArithmGenericFilter()
29 :ImageToImageGenericFilter<Self>("ImageArithmGenericFilter"),mTypeOfOperation(0) {
30 InitializeImageType<2>();
31 InitializeImageType<3>();
32 InitializeImageType<4>();
33 mIsOperationUseASecondImage = false;
35 //--------------------------------------------------------------------
37 //--------------------------------------------------------------------
38 template<unsigned int Dim>
39 void clitk::ImageArithmGenericFilter::InitializeImageType() {
40 ADD_IMAGE_TYPE(Dim, char);
41 ADD_IMAGE_TYPE(Dim, short);
42 /* ADD_IMAGE_TYPE(Dim, short);
43 ADD_IMAGE_TYPE(Dim, ushort;
44 ADD_IMAGE_TYPE(Dim, int);
45 ADD_IMAGE_TYPE(Dim, unsigned int);
46 ADD_IMAGE_TYPE(Dim, float);
47 ADD_IMAGE_TYPE(Dim, double);
50 //--------------------------------------------------------------------
53 //--------------------------------------------------------------------
54 template<class ImageType>
55 void clitk::ImageArithmGenericFilter::UpdateWithInputImageType() {
59 //DS TO BE CHANGED IN GetNumberOfInput() !!!!!!!!
61 if (mInputFilenames.size() == 0) {
62 std::cerr << "ERROR : please provide at least a input filename." << std::endl;
64 if (mInputFilenames.size() == 1) {
65 mIsOperationUseASecondImage = false;
67 if (mInputFilenames.size() == 2) {
68 mIsOperationUseASecondImage = true;
70 if (mInputFilenames.size() > 2) {
71 std::cerr << "ERROR : please provide only 1 or 2 input filenames." << std::endl;
75 // typename ImageType::Pointer input1 = clitk::readImage<ImageType>(mInputFilenames[0], mIOVerbose);
76 typename ImageType::Pointer input1 = this->template GetInput<ImageType>(0);
77 typename ImageType::Pointer outputImage;
79 // Read input2 (float is ok altough it could take too memory ...)
80 if (mIsOperationUseASecondImage) {
81 typedef itk::Image<float,ImageType::ImageDimension> ImageType2;
82 // typename ImageType2::Pointer input2 = clitk::readImage<ImageType2>(mInputFilenames[1], mIOVerbose);
83 typename ImageType2::Pointer input2 = this->template GetInput<ImageType2>(1);
84 outputImage = ComputeImage<ImageType, ImageType2>(input1, input2);
87 outputImage = ComputeImage<ImageType>(input1);
91 this->SetNextOutput<ImageType>(outputImage);
92 //clitk::writeImage<ImageType>(outputImage, mOutputFilename, mIOVerbose);
94 //--------------------------------------------------------------------
96 //--------------------------------------------------------------------
97 template<class ImageType>
98 typename ImageType::Pointer
99 clitk::ImageArithmGenericFilter::ComputeImage(typename ImageType::Pointer inputImage) {
101 typedef typename ImageType::PixelType PixelType;
102 typedef itk::ImageRegionIterator<ImageType> IteratorType;
103 IteratorType it(inputImage, inputImage->GetLargestPossibleRegion());
106 switch (mTypeOfOperation) {
108 while (!it.IsAtEnd()) {
109 it.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() + mScalar) );
114 while (!it.IsAtEnd()) {
115 it.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() * mScalar) );
120 while (!it.IsAtEnd()) {
122 it.Set(PixelTypeDownCast<double, PixelType>(mScalar / (double)it.Get()));
123 else it.Set(mDefaultPixelValue);
128 while (!it.IsAtEnd()) {
129 if (it.Get() < mScalar) it.Set(PixelTypeDownCast<double, PixelType>(mScalar));
134 while (!it.IsAtEnd()) {
135 if (it.Get() > mScalar) it.Set(PixelTypeDownCast<double, PixelType>(mScalar));
139 case 5: // Absolute value
140 while (!it.IsAtEnd()) {
141 if (it.Get() <= 0) it.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
142 // <= zero to avoid warning for unsigned types
146 case 6: // Squared value
147 while (!it.IsAtEnd()) {
148 it.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
153 while (!it.IsAtEnd()) {
155 it.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
156 else it.Set(mDefaultPixelValue);
161 while (!it.IsAtEnd()) {
162 it.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
167 while (!it.IsAtEnd()) {
169 it.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
171 if (it.Get() ==0) it.Set(0);
172 else it.Set(mDefaultPixelValue);
178 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
184 //--------------------------------------------------------------------
186 //--------------------------------------------------------------------
187 template<class ImageType1, class ImageType2>
188 typename ImageType1::Pointer
189 clitk::ImageArithmGenericFilter::ComputeImage(typename ImageType1::Pointer inputImage1,
190 typename ImageType2::Pointer inputImage2) {
192 typedef typename ImageType1::PixelType PixelType;
193 typedef itk::ImageRegionIterator<ImageType1> IteratorType;
194 IteratorType it1(inputImage1, inputImage1->GetLargestPossibleRegion());
197 typedef itk::ImageRegionConstIterator<ImageType2> ConstIteratorType;
198 ConstIteratorType it2(inputImage2, inputImage2->GetLargestPossibleRegion());
201 switch (mTypeOfOperation) {
203 while (!it1.IsAtEnd()) {
204 it1.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() + (double)it2.Get()) );
209 while (!it1.IsAtEnd()) {
210 it1.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() * (double)it2.Get()) );
215 while (!it1.IsAtEnd()) {
217 it1.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() / (double)it2.Get()));
218 else it1.Set(mDefaultPixelValue);
223 while (!it1.IsAtEnd()) {
224 if (it1.Get() < it2.Get()) it1.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
229 while (!it1.IsAtEnd()) {
230 if (it1.Get() > it2.Get()) it1.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
234 case 5: // Absolute difference
235 while (!it1.IsAtEnd()) {
236 it1.Set(PixelTypeDownCast<double, PixelType>(fabs(it2.Get()-it1.Get())));
240 case 6: // Squared differences
241 while (!it1.IsAtEnd()) {
242 it1.Set(PixelTypeDownCast<double, PixelType>(pow((double)it1.Get()-(double)it2.Get(),2)));
246 case 7: // Difference
247 while (!it1.IsAtEnd()) {
248 it1.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get()-(double)it2.Get()));
252 case 8: // Relative Difference
253 while (!it1.IsAtEnd()) {
254 if (it1.Get() != 0) it1.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
260 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
266 //--------------------------------------------------------------------
268 #endif //define CLITKIMAGEARITHMGENERICFILTER_CXX