From: Romulo Pinho Date: Thu, 31 Jan 2013 14:54:40 +0000 (+0100) Subject: Added clitkVectorArithm tool X-Git-Tag: v1.4.0~249^2~1^2~1 X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=2e890f7bcf7118cd1b49ff824fa8b5b6a6d9f044;p=clitk.git Added clitkVectorArithm tool - removed all vector support from clitkImageArithm --- diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 36fb8d9..7193772 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -11,6 +11,9 @@ ADD_LIBRARY(clitkBinarizeImageLib clitkBinarizeImageGenericFilter.cxx ${clitkBin WRAP_GGO(clitkImageArithm_GGO_C clitkImageArithm.ggo) ADD_LIBRARY(clitkImageArithmImageLib clitkImageArithmGenericFilter.cxx ${clitkImageArithm_GGO_C}) +WRAP_GGO(clitkVectorArithm_GGO_C clitkVectorArithm.ggo) +ADD_LIBRARY(clitkVectorArithmLib clitkVectorArithmGenericFilter.cxx ${clitkVectorArithm_GGO_C}) + WRAP_GGO(clitkResampleImage_GGO_C clitkResampleImage.ggo) ADD_LIBRARY(clitkResampleImageLib clitkResampleImageGenericFilter.cxx ${clitkResampleImage_GGO_C}) @@ -114,6 +117,10 @@ IF (CLITK_BUILD_TOOLS) TARGET_LINK_LIBRARIES(clitkImageArithm clitkImageArithmImageLib clitkCommon ${ITK_LIBRARIES} ) SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkImageArithm) + ADD_EXECUTABLE(clitkVectorArithm clitkVectorArithm.cxx) + TARGET_LINK_LIBRARIES(clitkVectorArithm clitkVectorArithmLib clitkCommon ${ITK_LIBRARIES} ) + SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkVectorArithm) + WRAP_GGO(clitkUnsharpMask_GGO_C clitkUnsharpMask.ggo) ADD_EXECUTABLE(clitkUnsharpMask clitkUnsharpMask.cxx ${clitkUnsharpMask_GGO_C}) TARGET_LINK_LIBRARIES(clitkUnsharpMask clitkCommon ${ITK_LIBRARIES} ) diff --git a/tools/clitkDicom2Image.cxx b/tools/clitkDicom2Image.cxx index 07a8565..2819fde 100644 --- a/tools/clitkDicom2Image.cxx +++ b/tools/clitkDicom2Image.cxx @@ -30,6 +30,8 @@ #include #endif +#include + //==================================================================== int main(int argc, char * argv[]) { @@ -52,8 +54,11 @@ int main(int argc, char * argv[]) //=========================================== /// Get slices locations ... - std::vector theorigin(3); - std::vector sliceLocations; + int series_number = -1; + std::set series_numbers; + std::map< int, std::vector > theorigin; + std::map< int, std::vector > sliceLocations; + std::map< int, std::vector > seriesFiles; for(unsigned int i=0; i pixel_size; gdcm::DataSet& ds = hreader.GetFile().GetDataSet(); + + if (args_info.extract_series_flag) { + gdcm::Attribute<0x20,0x11> series_number_att; + series_number_att.SetFromDataSet(hreader.GetFile().GetDataSet()); + series_number = series_number_att.GetValue(); + } + + series_numbers.insert(series_number); + theorigin[series_number] = gdcm::ImageHelper::GetOriginValue(hreader.GetFile()); + std::cout << "theorigin[series_number][0] " << theorigin[series_number][0] << std::endl; + std::cout << "theorigin[series_number][1] " << theorigin[series_number][1] << std::endl; + std::cout << "theorigin[series_number][2] " << theorigin[series_number][2] << std::endl; + sliceLocations[series_number].push_back(theorigin[series_number][2]); + seriesFiles[series_number].push_back(input_files[i]); + + gdcm::Attribute<0x28, 0x100> pixel_size; pixel_size.SetFromDataSet(ds); if (pixel_size.GetValue() != 16) { @@ -80,10 +98,18 @@ int main(int argc, char * argv[]) header->SetFileName(input_files[i]); header->SetMaxSizeLoadEntry(16384); // required ? header->Load(); - theorigin[0] = header->GetXOrigin(); - theorigin[1] = header->GetYOrigin(); - theorigin[2] = header->GetZOrigin(); - sliceLocations.push_back(theorigin[2]); + + if (args_info.extract_series_flag) { + series_number = atoi(header->GetEntryValue(0x20,0x11).c_str()); + } + + series_numbers.insert(series_number); + theorigin[series_number].resize(3); + theorigin[series_number][0] = header->GetXOrigin(); + theorigin[series_number][1] = header->GetYOrigin(); + theorigin[series_number][2] = header->GetZOrigin(); + sliceLocations[series_number].push_back(theorigin[series_number][2]); + seriesFiles[series_number].push_back(input_files[i]); if (header->GetPixelSize() != 2) { std::cerr << "Pixel type 2 bytes ! " << std::endl; std::cerr << "In file " << input_files[i] << std::endl; @@ -94,87 +120,103 @@ int main(int argc, char * argv[]) //=========================================== // Sort slices locations ... - std::vector sliceIndex; - clitk::GetSortedIndex(sliceLocations, sliceIndex); - if (args_info.verboseSliceLocation_flag) { - std::cout << sliceLocations[sliceIndex[0]] << " -> " - << sliceIndex[0] << " / " << 0 << " => " - << "0 mm " - << input_files[sliceIndex[0]] - << std::endl; - for(unsigned int i=1; i " - << sliceIndex[i] << " / " << i << " => " - << sliceLocations[sliceIndex[i]] - sliceLocations[sliceIndex[i-1]] - << "mm " - << input_files[sliceIndex[i]] + std::set::iterator sn = series_numbers.begin(); + while ( sn != series_numbers.end() ) { + std::vector locs = sliceLocations[*sn]; + std::vector origin = theorigin[*sn]; + std::vector files = seriesFiles[*sn]; + std::vector sliceIndex; + clitk::GetSortedIndex(locs, sliceIndex); + if (args_info.verboseSliceLocation_flag) { + std::cout << locs[sliceIndex[0]] << " -> " + << sliceIndex[0] << " / " << 0 << " => " + << "0 mm " + << files[sliceIndex[0]] << std::endl; - } - } - - //=========================================== - // Analyze slices locations ... - double currentDist; - double dist=0; - double tolerance = args_info.tolerance_arg; - double previous = sliceLocations[sliceIndex[0]]; - for(unsigned int i=1; i tolerance) { - std::cout << "ERROR : " << std::endl - << "Current slice pos is = " << sliceLocations[sliceIndex[i]] << std::endl - << "Previous slice pos is = " << previous << std::endl - << "Current file is = " << input_files[sliceIndex[i]] << std::endl - << "Current index is = " << i << std::endl - << "Current sortindex is = " << sliceIndex[i] << std::endl - << "Current slice diff = " << dist << std::endl - << "Current error = " << fabs(dist-currentDist) << std::endl; - exit(1); + for(unsigned int i=1; i " + << sliceIndex[i] << " / " << i << " => " + << locs[sliceIndex[i]] - locs[sliceIndex[i-1]] + << "mm " + << files[sliceIndex[i]] + << std::endl; } - } else dist = currentDist; - previous = sliceLocations[sliceIndex[i]]; - } + } - //=========================================== - // Create ordered vector of filenames - std::vector sorted_files; - sorted_files.resize(sliceIndex.size()); - for(unsigned int i=0; i tolerance) { + std::cout << "ERROR : " << std::endl + << "Current slice pos is = " << locs[sliceIndex[i]] << std::endl + << "Previous slice pos is = " << previous << std::endl + << "Current file is = " << files[sliceIndex[i]] << std::endl + << "Current index is = " << i << std::endl + << "Current sortindex is = " << sliceIndex[i] << std::endl + << "Current slice diff = " << dist << std::endl + << "Current error = " << fabs(dist-currentDist) << std::endl; + exit(1); + } + } else dist = currentDist; + previous = locs[sliceIndex[i]]; + } - //=========================================== - // Read write serie - vvImageReader::Pointer reader = vvImageReader::New(); - reader->SetInputFilenames(sorted_files); - reader->Update(vvImageReader::DICOM); - if (reader->GetLastError().size() != 0) { - std::cerr << reader->GetLastError() << std::endl; - return 1; - } - - vvImage::Pointer image = reader->GetOutput(); - vtkImageData* vtk_image = image->GetFirstVTKImageData(); - vtkImageChangeInformation* modifier = vtkImageChangeInformation::New(); - if (args_info.focal_origin_given) { - std::vector spacing = image->GetSpacing(); - std::vector size = image->GetSize(); - theorigin[0] = -spacing[0]*size[0]/2.0; - theorigin[1] = -spacing[1]*size[1]/2.0; - modifier->SetInput(vtk_image); - modifier->SetOutputOrigin(theorigin[0], theorigin[1], sliceLocations[sliceIndex[0]]); - modifier->Update(); - vvImage::Pointer focal_image = vvImage::New(); - focal_image->AddVtkImage(modifier->GetOutput()); - image = focal_image; - } + //=========================================== + // Create ordered vector of filenames + std::vector sorted_files; + sorted_files.resize(sliceIndex.size()); + for(unsigned int i=0; iSetInput(image); - writer->SetOutputFileName(args_info.output_arg); - writer->Update(); + //=========================================== + // Read write serie + vvImageReader::Pointer reader = vvImageReader::New(); + reader->SetInputFilenames(sorted_files); + reader->Update(vvImageReader::DICOM); + if (reader->GetLastError().size() != 0) { + std::cerr << reader->GetLastError() << std::endl; + return 1; + } + + vvImage::Pointer image = reader->GetOutput(); + vtkImageData* vtk_image = image->GetFirstVTKImageData(); + vtkImageChangeInformation* modifier = vtkImageChangeInformation::New(); + if (args_info.focal_origin_given) { + std::vector spacing = image->GetSpacing(); + std::vector size = image->GetSize(); + origin[0] = -spacing[0]*size[0]/2.0; + origin[1] = -spacing[1]*size[1]/2.0; + modifier->SetInput(vtk_image); + modifier->SetOutputOrigin(origin[0], origin[1], locs[sliceIndex[0]]); + modifier->Update(); + vvImage::Pointer focal_image = vvImage::New(); + focal_image->AddVtkImage(modifier->GetOutput()); + image = focal_image; + } - modifier->Delete(); + std::string outfile; + if (series_numbers.size() == 1) + outfile = args_info.output_arg; + else { + std::ostringstream name; + name << *sn << "_" << args_info.output_arg; + outfile = name.str(); + } + vvImageWriter::Pointer writer = vvImageWriter::New(); + writer->SetInput(image); + writer->SetOutputFileName(outfile); + writer->Update(); + modifier->Delete(); + + sn++; + } + return 0; } diff --git a/tools/clitkDicom2Image.ggo b/tools/clitkDicom2Image.ggo index fb51135..a2d3a89 100644 --- a/tools/clitkDicom2Image.ggo +++ b/tools/clitkDicom2Image.ggo @@ -10,3 +10,4 @@ option "tolerance" t "Tolerance for slice position" double default="0" no option "output" o "Output image filename" string yes option "std_input" - "Take the like of input file from stdin, to support huge lists of filenames" flag off option "focal_origin" - "Output files with FOCAL-like origin, instead of the origin present in the dicom files" flag off +option "extract_series" s "Identify different series in the file list and create the MHDs accordinly" flag off diff --git a/tools/clitkImageArithmGenericFilter.cxx b/tools/clitkImageArithmGenericFilter.cxx index f0d55c2..2509b1a 100644 --- a/tools/clitkImageArithmGenericFilter.cxx +++ b/tools/clitkImageArithmGenericFilter.cxx @@ -25,669 +25,6 @@ namespace clitk { // template<> // class ImageArithmGenericFilter; -//////////////////////////////////////////////////////////////////// - // specializations for itk::Vector, 3u - template<> - template<> - void ImageArithmGenericFilter::UpdateWithInputImageType< itk::Image< itk::Vector, 3u > >() - { - typedef itk::Image< itk::Vector, 3u > ImageType; - - // Read input1 - ImageType::Pointer input1 = this->GetInput(0); - - // Set input image iterator - typedef itk::ImageRegionIterator IteratorType; - IteratorType it(input1, input1->GetLargestPossibleRegion()); - - // typedef input2 - ImageType::Pointer input2 = NULL; - IteratorType it2; - - /* - // Special case for normalisation - if (mTypeOfOperation == 12) { - typedef itk::MinimumMaximumImageCalculator MinMaxFilterType; - typename MinMaxFilterType::Pointer ff = MinMaxFilterType::New(); - ff->SetImage(input1); - ff->ComputeMaximum(); - mScalar = ff->GetMaximum(); - mTypeOfOperation = 11; // divide - } - */ - - if (mIsOperationUseASecondImage) { - // Read input2 - input2 = this->GetInput(1); - // Set input image iterator - it2 = IteratorType(input2, input2->GetLargestPossibleRegion()); - // Check dimension - if (!clitk::HaveSameSize(input1, input2)) { - itkExceptionMacro(<< "The images (input and input2) must have the same size"); - } - if(!clitk::HaveSameSpacing(input1, input2)) { - itkWarningMacro(<< "The images (input and input2) do not have the same spacing. " - << "Using first input's information."); - } - } - - /* - // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite - if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) { - // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is " - // << typeid(PixelType).name() - // << std::endl; - mOverwriteInputImage = false; - } - */ - - // ---------------- Overwrite input Image --------------------- - if (mOverwriteInputImage) { - // Set output iterator (to input1) - IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion()); - if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito); - else ComputeImage(it, ito); - this->SetNextOutput(input1); - } - // ---------------- Create new output Image --------------------- - else { - /*if (mOutputIsFloat) { - // Create output image - typedef itk::Image OutputImageType; - typename OutputImageType::Pointer output = OutputImageType::New(); - output->SetRegions(input1->GetLargestPossibleRegion()); - output->SetOrigin(input1->GetOrigin()); - output->SetSpacing(input1->GetSpacing()); - output->Allocate(); - // Set output iterator - typedef itk::ImageRegionIterator IteratorOutputType; - IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion()); - if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito); - else ComputeImage(it, ito); - this->template SetNextOutput(output); - } else*/ { - // Create output image - typedef ImageType OutputImageType; - OutputImageType::Pointer output = OutputImageType::New(); - output->SetRegions(input1->GetLargestPossibleRegion()); - output->SetOrigin(input1->GetOrigin()); - output->SetSpacing(input1->GetSpacing()); - output->Allocate(); - // Set output iterator - typedef itk::ImageRegionIterator IteratorOutputType; - IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion()); - if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito); - else ComputeImage(it, ito); - this->SetNextOutput(output); - } - } - } - - template<> - template<> - void ImageArithmGenericFilter::ComputeImage< - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > > - (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito) - { - typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter2; - - ito.GoToBegin(); - it.GoToBegin(); - - typedef Iter2::PixelType PixelType; - - PixelType scalar_vector; - scalar_vector.Fill(mScalar); - - // Perform operation - switch (mTypeOfOperation) { - case 0: // Addition - while (!it.IsAtEnd()) { - ito.Set(it.Get() + scalar_vector); - ++it; - ++ito; - } - break; - case 1: // Multiply - while (!it.IsAtEnd()) { - ito.Set(it.Get() * mScalar); - ++it; - ++ito; - } - break; - /* - case 2: // Inverse - while (!it.IsAtEnd()) { - if (it.Get() != 0) - ito.Set(mScalar / it.Get())); - else ito.Set(mDefaultPixelValue); - ++it; - ++ito; - } - break; - case 3: // Max - while (!it.IsAtEnd()) { - if (it.Get() < mScalar) ito.Set(PixelTypeDownCast(mScalar)); - else ito.Set(PixelTypeDownCast(it.Get())); - ++it; - ++ito; - } - break; - case 4: // Min - while (!it.IsAtEnd()) { - if (it.Get() > mScalar) ito.Set(PixelTypeDownCast(mScalar)); - else ito.Set(PixelTypeDownCast(it.Get())); - ++it; - ++ito; - } - break; - case 5: // Absolute value - while (!it.IsAtEnd()) { - if (it.Get() <= 0) ito.Set(PixelTypeDownCast(-it.Get())); - // <= zero to avoid warning for unsigned types - else ito.Set(PixelTypeDownCast(it.Get())); - ++it; - ++ito; - } - break; - case 6: // Squared value - while (!it.IsAtEnd()) { - ito.Set(PixelTypeDownCast((double)it.Get()*(double)it.Get())); - ++it; - ++ito; - } - break; - case 7: // Log - while (!it.IsAtEnd()) { - if (it.Get() > 0) - ito.Set(PixelTypeDownCast(log((double)it.Get()))); - else ito.Set(mDefaultPixelValue); - ++it; - ++ito; - } - break; - case 8: // exp - while (!it.IsAtEnd()) { - ito.Set(PixelTypeDownCast(exp((double)it.Get()))); - ++it; - ++ito; - } - break; - case 9: // sqrt - while (!it.IsAtEnd()) { - if (it.Get() > 0) - ito.Set(PixelTypeDownCast(sqrt((double)it.Get()))); - else { - if (it.Get() ==0) ito.Set(0); - else ito.Set(mDefaultPixelValue); - } - ++it; - ++ito; - } - break; - case 10: // exp - while (!it.IsAtEnd()) { - ito.Set(PixelTypeDownCast((0x10000 - (double)it.Get())/mScalar)); - ++it; - ++ito; - } - break; - */ - case 11: // divide - while (!it.IsAtEnd()) { - ito.Set(it.Get() / mScalar); - ++it; - ++ito; - } - break; - default: // error ? - std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl; - exit(-1); - } - - } - - template<> - template<> - void ImageArithmGenericFilter::ComputeImage< - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > - > - (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it1, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it2, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito) - { - typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter1; - typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter2; - typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter3; - - it1.GoToBegin(); - it2.GoToBegin(); - ito.GoToBegin(); - typedef Iter3::PixelType PixelType; - - switch (mTypeOfOperation) { - case 0: // Addition - while (!ito.IsAtEnd()) { - ito.Set(it1.Get() + it2.Get()); - ++it1; - ++it2; - ++ito; - } - break; - /* - case 1: // Multiply - while (!ito.IsAtEnd()) { - ito.Set(it1.Get() * it2.Get()) ); - ++it1; - ++it2; - ++ito; - } - break; - case 2: // Divide - while (!ito.IsAtEnd()) { - if (it1.Get() != 0) - ito.Set(it1.Get() / it2.Get())); - else ito.Set(mDefaultPixelValue); - ++it1; - ++it2; - ++ito; - } - break; - case 3: // Max - while (!ito.IsAtEnd()) { - if (it1.Get() < it2.Get()) ito.Set(it2.Get()); - ++it1; - ++it2; - ++ito; - } - break; - case 4: // Min - while (!ito.IsAtEnd()) { - if (it1.Get() > it2.Get()) ito.Set(it2.Get()); - ++it1; - ++it2; - ++ito; - } - break; - */ - case 5: // Absolute difference - while (!ito.IsAtEnd()) { - ito.Set(it2.Get()-it1.Get()); - ++it1; - ++it2; - ++ito; - } - break; - /* - case 6: // Squared differences - while (!ito.IsAtEnd()) { - ito.Set(pow(it1.Get()-it2.Get(),2))); - ++it1; - ++it2; - ++ito; - } - break; - */ - case 7: // Difference - while (!ito.IsAtEnd()) { - ito.Set(it1.Get()-it2.Get()); - ++it1; - ++it2; - ++ito; - } - break; - /* - case 8: // Relative Difference - while (!ito.IsAtEnd()) { - if (it1.Get() != 0) ito.Set(PixelTypeDownCast(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get()); - else ito.Set(0.0); - ++it1; - ++it2; - ++ito; - } - break; - */ - default: // error ? - std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl; - exit(-1); - } - } - -//////////////////////////////////////////////////////////////////// - // specializations for itk::Vector, 3u - template<> - template<> - void ImageArithmGenericFilter::UpdateWithInputImageType< itk::Image< itk::Vector, 3u > >() - { - typedef itk::Image< itk::Vector, 3u > ImageType; - - // Read input1 - ImageType::Pointer input1 = this->GetInput(0); - - // Set input image iterator - typedef itk::ImageRegionIterator IteratorType; - IteratorType it(input1, input1->GetLargestPossibleRegion()); - - // typedef input2 - ImageType::Pointer input2 = NULL; - IteratorType it2; - - /* - // Special case for normalisation - if (mTypeOfOperation == 12) { - typedef itk::MinimumMaximumImageCalculator MinMaxFilterType; - typename MinMaxFilterType::Pointer ff = MinMaxFilterType::New(); - ff->SetImage(input1); - ff->ComputeMaximum(); - mScalar = ff->GetMaximum(); - mTypeOfOperation = 11; // divide - } - */ - - if (mIsOperationUseASecondImage) { - // Read input2 - input2 = this->GetInput(1); - // Set input image iterator - it2 = IteratorType(input2, input2->GetLargestPossibleRegion()); - // Check dimension - if (!clitk::HaveSameSize(input1, input2)) { - itkExceptionMacro(<< "The images (input and input2) must have the same size"); - } - if(!clitk::HaveSameSpacing(input1, input2)) { - itkWarningMacro(<< "The images (input and input2) do not have the same spacing. " - << "Using first input's information."); - } - } - - /* - // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite - if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) { - // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is " - // << typeid(PixelType).name() - // << std::endl; - mOverwriteInputImage = false; - } - */ - - // ---------------- Overwrite input Image --------------------- - if (mOverwriteInputImage) { - // Set output iterator (to input1) - IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion()); - if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito); - else ComputeImage(it, ito); - this->SetNextOutput(input1); - } - // ---------------- Create new output Image --------------------- - else { - /*if (mOutputIsFloat) { - // Create output image - typedef itk::Image OutputImageType; - typename OutputImageType::Pointer output = OutputImageType::New(); - output->SetRegions(input1->GetLargestPossibleRegion()); - output->SetOrigin(input1->GetOrigin()); - output->SetSpacing(input1->GetSpacing()); - output->Allocate(); - // Set output iterator - typedef itk::ImageRegionIterator IteratorOutputType; - IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion()); - if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito); - else ComputeImage(it, ito); - this->template SetNextOutput(output); - } else*/ { - // Create output image - typedef ImageType OutputImageType; - OutputImageType::Pointer output = OutputImageType::New(); - output->SetRegions(input1->GetLargestPossibleRegion()); - output->SetOrigin(input1->GetOrigin()); - output->SetSpacing(input1->GetSpacing()); - output->Allocate(); - // Set output iterator - typedef itk::ImageRegionIterator IteratorOutputType; - IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion()); - if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito); - else ComputeImage(it, ito); - this->SetNextOutput(output); - } - } - } - - template<> - template<> - void ImageArithmGenericFilter::ComputeImage< - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > > - (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito) - { - typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter2; - - ito.GoToBegin(); - it.GoToBegin(); - - typedef Iter2::PixelType PixelType; - - PixelType scalar_vector; - scalar_vector.Fill(mScalar); - - // Perform operation - switch (mTypeOfOperation) { - case 0: // Addition - while (!it.IsAtEnd()) { - ito.Set(it.Get() + scalar_vector); - ++it; - ++ito; - } - break; - case 1: // Multiply - while (!it.IsAtEnd()) { - ito.Set(it.Get() * mScalar); - ++it; - ++ito; - } - break; - /* - case 2: // Inverse - while (!it.IsAtEnd()) { - if (it.Get() != 0) - ito.Set(mScalar / it.Get())); - else ito.Set(mDefaultPixelValue); - ++it; - ++ito; - } - break; - case 3: // Max - while (!it.IsAtEnd()) { - if (it.Get() < mScalar) ito.Set(PixelTypeDownCast(mScalar)); - else ito.Set(PixelTypeDownCast(it.Get())); - ++it; - ++ito; - } - break; - case 4: // Min - while (!it.IsAtEnd()) { - if (it.Get() > mScalar) ito.Set(PixelTypeDownCast(mScalar)); - else ito.Set(PixelTypeDownCast(it.Get())); - ++it; - ++ito; - } - break; - case 5: // Absolute value - while (!it.IsAtEnd()) { - if (it.Get() <= 0) ito.Set(PixelTypeDownCast(-it.Get())); - // <= zero to avoid warning for unsigned types - else ito.Set(PixelTypeDownCast(it.Get())); - ++it; - ++ito; - } - break; - case 6: // Squared value - while (!it.IsAtEnd()) { - ito.Set(PixelTypeDownCast((double)it.Get()*(double)it.Get())); - ++it; - ++ito; - } - break; - case 7: // Log - while (!it.IsAtEnd()) { - if (it.Get() > 0) - ito.Set(PixelTypeDownCast(log((double)it.Get()))); - else ito.Set(mDefaultPixelValue); - ++it; - ++ito; - } - break; - case 8: // exp - while (!it.IsAtEnd()) { - ito.Set(PixelTypeDownCast(exp((double)it.Get()))); - ++it; - ++ito; - } - break; - case 9: // sqrt - while (!it.IsAtEnd()) { - if (it.Get() > 0) - ito.Set(PixelTypeDownCast(sqrt((double)it.Get()))); - else { - if (it.Get() ==0) ito.Set(0); - else ito.Set(mDefaultPixelValue); - } - ++it; - ++ito; - } - break; - case 10: // exp - while (!it.IsAtEnd()) { - ito.Set(PixelTypeDownCast((0x10000 - (double)it.Get())/mScalar)); - ++it; - ++ito; - } - break; - */ - case 11: // divide - while (!it.IsAtEnd()) { - ito.Set(it.Get() / mScalar); - ++it; - ++ito; - } - break; - default: // error ? - std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl; - exit(-1); - } - - } - - template<> - template<> - void ImageArithmGenericFilter::ComputeImage< - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > - > - (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it1, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it2, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito) - { - typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter1; - typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter2; - typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter3; - - it1.GoToBegin(); - it2.GoToBegin(); - ito.GoToBegin(); - typedef Iter3::PixelType PixelType; - - switch (mTypeOfOperation) { - case 0: // Addition - while (!ito.IsAtEnd()) { - ito.Set(it1.Get() + it2.Get()); - ++it1; - ++it2; - ++ito; - } - break; - /* - case 1: // Multiply - while (!ito.IsAtEnd()) { - ito.Set(it1.Get() * it2.Get()) ); - ++it1; - ++it2; - ++ito; - } - break; - case 2: // Divide - while (!ito.IsAtEnd()) { - if (it1.Get() != 0) - ito.Set(it1.Get() / it2.Get())); - else ito.Set(mDefaultPixelValue); - ++it1; - ++it2; - ++ito; - } - break; - case 3: // Max - while (!ito.IsAtEnd()) { - if (it1.Get() < it2.Get()) ito.Set(it2.Get()); - ++it1; - ++it2; - ++ito; - } - break; - case 4: // Min - while (!ito.IsAtEnd()) { - if (it1.Get() > it2.Get()) ito.Set(it2.Get()); - ++it1; - ++it2; - ++ito; - } - break; - */ - case 5: // Absolute difference - while (!ito.IsAtEnd()) { - ito.Set(it2.Get()-it1.Get()); - ++it1; - ++it2; - ++ito; - } - break; - /* - case 6: // Squared differences - while (!ito.IsAtEnd()) { - ito.Set(pow(it1.Get()-it2.Get(),2))); - ++it1; - ++it2; - ++ito; - } - break; - */ - case 7: // Difference - while (!ito.IsAtEnd()) { - ito.Set(it1.Get()-it2.Get()); - ++it1; - ++it2; - ++ito; - } - break; - /* - case 8: // Relative Difference - while (!ito.IsAtEnd()) { - if (it1.Get() != 0) ito.Set(PixelTypeDownCast(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get()); - else ito.Set(0.0); - ++it1; - ++it2; - ++ito; - } - break; - */ - default: // error ? - std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl; - exit(-1); - } - } } diff --git a/tools/clitkImageArithmGenericFilter.h b/tools/clitkImageArithmGenericFilter.h index 39bfeb9..3c3579f 100644 --- a/tools/clitkImageArithmGenericFilter.h +++ b/tools/clitkImageArithmGenericFilter.h @@ -97,49 +97,6 @@ namespace clitk { }; // end class ImageArithmGenericFilter - // specializations for itk::Vector, 3u - template<> template<> - void ImageArithmGenericFilter::UpdateWithInputImageType< itk::Image< itk::Vector, 3u > >(); - - template<> template<> - void ImageArithmGenericFilter::ComputeImage< - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > - > - (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito); - - template<> template<> - void ImageArithmGenericFilter::ComputeImage< - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > - > - (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it1, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it2, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito); - - // specializations for itk::Vector, 3u - template<> template<> - void ImageArithmGenericFilter::UpdateWithInputImageType< itk::Image< itk::Vector, 3u > >(); - - template<> template<> - void ImageArithmGenericFilter::ComputeImage< - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > - > - (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito); - - template<> template<> - void ImageArithmGenericFilter::ComputeImage< - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > - > - (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it1, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it2, - itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito); } // end namespace //-------------------------------------------------------------------- diff --git a/tools/clitkImageArithmGenericFilter.txx b/tools/clitkImageArithmGenericFilter.txx index 1316fc1..4c2f0a5 100644 --- a/tools/clitkImageArithmGenericFilter.txx +++ b/tools/clitkImageArithmGenericFilter.txx @@ -45,8 +45,6 @@ template void ImageArithmGenericFilter::InitializeImageType() { ADD_DEFAULT_IMAGE_TYPES(Dim); - ADD_VEC_IMAGE_TYPE(3u,3u,float); - ADD_VEC_IMAGE_TYPE(3u,3u,double); } //-------------------------------------------------------------------- diff --git a/tools/clitkVectorArithm.cxx b/tools/clitkVectorArithm.cxx new file mode 100644 index 0000000..bc3ac3a --- /dev/null +++ b/tools/clitkVectorArithm.cxx @@ -0,0 +1,52 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://www.centreleonberard.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html +===========================================================================**/ +#ifndef CLITKVECTORARITHM_CXX +#define CLITKVECTORARITHM_CXX +/** + ------------------------------------------------- + * @file clitkVectorArithm.cxx + * @author David Sarrut + * @date 23 Feb 2008 08:37:53 + -------------------------------------------------*/ + +// clitk include +#include "clitkVectorArithm_ggo.h" +#include "clitkVectorArithmGenericFilter.h" +#include "clitkIO.h" + +//-------------------------------------------------------------------- +int main(int argc, char * argv[]) +{ + + // Init command line + GGO(clitkVectorArithm, args_info); + CLITK_INIT; + + // Creation of a generic filter + typedef clitk::VectorArithmGenericFilter FilterType; + FilterType::Pointer filter = FilterType::New(); + + // Go ! + filter->SetArgsInfo(args_info); + CLITK_TRY_CATCH_EXIT(filter->Update()); + + // this is the end my friend + return EXIT_SUCCESS; +} // end main + +#endif //define CLITKVECTORARITHM_CXX diff --git a/tools/clitkVectorArithm.ggo b/tools/clitkVectorArithm.ggo new file mode 100644 index 0000000..5c7d65e --- /dev/null +++ b/tools/clitkVectorArithm.ggo @@ -0,0 +1,21 @@ +#File clitkVectorArithm.ggo +package "clitkVectorArithm" +version "1.0" +purpose "Perform an arithmetic operation (+-*/ ...) using two images or using an image and a scalar value." + + +option "config" - "Config file" string no +option "verbose" v "Verbose" flag off +option "imagetypes" - "Display allowed image types" flag off + +option "input1" i "Input first image filename" string yes +option "input2" j "Input second image filename" string no +option "output" o "Output image filename" string yes + +option "scalar" s "Scalar value" double no +option "operation" t "Type of operation : \n With another image : 0=add, 1=multiply (dotproduct), 7=difference, 9=crossproduct\n; For 'scalar' : 0=add, 1=multiply, 5=absval (magnitude), 6=squared magnitude, 11=divide, 12=normalize" int default="0" no +option "pixelValue" - "Default value for NaN/Inf" double default="0.0" no + +option "setFloatOutput" f "Set output to float pixel type" flag off + +#option "binary" b "Mask for binary operation" string no diff --git a/tools/clitkVectorArithmGenericFilter.cxx b/tools/clitkVectorArithmGenericFilter.cxx new file mode 100644 index 0000000..ff7b2d7 --- /dev/null +++ b/tools/clitkVectorArithmGenericFilter.cxx @@ -0,0 +1,31 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://www.centreleonberard.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html +===========================================================================**/ +#ifndef CLITKVECTORARITHMGENERICFILTER_CXX +#define CLITKVECTORARITHMGENERICFILTER_CXX + +#include "clitkVectorArithmGenericFilter.h" + +namespace clitk { + // Specialisation +// template<> +// class VectorArithmGenericFilter; + + +} + +#endif //define CLITKVECTORARITHMGENERICFILTER_CXX diff --git a/tools/clitkVectorArithmGenericFilter.h b/tools/clitkVectorArithmGenericFilter.h new file mode 100644 index 0000000..1196b71 --- /dev/null +++ b/tools/clitkVectorArithmGenericFilter.h @@ -0,0 +1,117 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://www.centreleonberard.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html +===========================================================================**/ +#ifndef CLITKVectorArithmGENERICFILTER_H +#define CLITKVectorArithmGENERICFILTER_H +/** + ------------------------------------------------------------------- + * @file clitkVectorArithmGenericFilter.h + * @author David Sarrut + * @date 23 Feb 2008 08:37:53 + + * @brief + -------------------------------------------------------------------*/ + +// clitk include +#include "clitkCommon.h" +#include "clitkImageToImageGenericFilter.h" +#include "clitkVectorArithm_ggo.h" + +// itk include +#include "itkImage.h" +#include "itkImageIOBase.h" +#include "itkImageRegionIterator.h" +#include "itkImageRegionConstIterator.h" + +//-------------------------------------------------------------------- +namespace clitk { + + template + class ITK_EXPORT VectorArithmGenericFilter: + public clitk::ImageToImageGenericFilter > { + + public: + + // Constructor + VectorArithmGenericFilter (); + + // Types + typedef VectorArithmGenericFilter Self; + typedef ImageToImageGenericFilterBase Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + // New + itkNewMacro(Self); + + //-------------------------------------------------------------------- + void SetArgsInfo(const args_info_type & a); + + // Set methods + void SetDefaultPixelValue (double value) { mDefaultPixelValue = value ;} + void SetTypeOfOperation (int value) { mTypeOfOperation = value ;} + void SetScalar (double value) { mScalar = value ;} + void EnableOverwriteInputImage(bool b); + + // Get methods + double GetDefaultPixelValue () { return mDefaultPixelValue ;} + int GetTypeOfOperation () { return mTypeOfOperation ;} + double GetScalar () { return mScalar ;} + + //-------------------------------------------------------------------- + // Main function called each time the filter is updated + template + void UpdateWithInputImageType(); + + protected: + template void InitializeImageType(); + bool mIsOperationUseASecondImage; + double mScalar; + double mDefaultPixelValue; + int mTypeOfOperation; + args_info_type mArgsInfo; + bool mOverwriteInputImage; + bool mOutputIsFloat; + bool mIsOutputScalar; + + template + void ComputeImage(Iter1 it, Iter2 ito); + + template + void ComputeImage(Iter1 it1, Iter2 it2, Iter3 ito); + + template + void ComputeScalarImage(Iter1 it, Iter2 ito); + + template + void ComputeScalarImage(Iter1 it1, Iter2 it2, Iter3 ito); + + //-------------------------------------------------------------------- + + }; // end class VectorArithmGenericFilter + + // specializations for itk::Vector, 3u +} // end namespace +//-------------------------------------------------------------------- + + +#ifndef ITK_MANUAL_INSTANTIATION +#include "clitkVectorArithmGenericFilter.txx" +#endif + +#endif //#define CLITKVectorArithmGENERICFILTER_H + diff --git a/tools/clitkVectorArithmGenericFilter.txx b/tools/clitkVectorArithmGenericFilter.txx new file mode 100644 index 0000000..38c3b87 --- /dev/null +++ b/tools/clitkVectorArithmGenericFilter.txx @@ -0,0 +1,510 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://www.centreleonberard.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html +===========================================================================**/ +#ifndef CLITKVectorArithmGENERICFILTER_TXX +#define CLITKVectorArithmGENERICFILTER_TXX + +#include "clitkImageCommon.h" + +#include "itkMinimumMaximumImageCalculator.h" + +namespace clitk +{ + +//-------------------------------------------------------------------- +template +VectorArithmGenericFilter::VectorArithmGenericFilter() + :ImageToImageGenericFilter("VectorArithmGenericFilter"),mTypeOfOperation(0) +{ + InitializeImageType<3>(); + mIsOperationUseASecondImage = false; + mIsOutputScalar = false; + mOverwriteInputImage = true; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +template +void VectorArithmGenericFilter::InitializeImageType() +{ + ADD_VEC_IMAGE_TYPE(Dim,3u,float); + ADD_VEC_IMAGE_TYPE(Dim,3u,double); + ADD_VEC_IMAGE_TYPE(Dim,2u,float); + ADD_VEC_IMAGE_TYPE(Dim,2u,double); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void VectorArithmGenericFilter::EnableOverwriteInputImage(bool b) +{ + mOverwriteInputImage = b; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void VectorArithmGenericFilter::SetArgsInfo(const args_info_type & a) +{ + mArgsInfo=a; + + // Set value + this->SetIOVerbose(mArgsInfo.verbose_flag); + mTypeOfOperation = mArgsInfo.operation_arg; + mDefaultPixelValue = mArgsInfo.pixelValue_arg; + mScalar = mArgsInfo.scalar_arg; + mOutputIsFloat = mArgsInfo.setFloatOutput_flag; + + if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes(); + + if (mArgsInfo.input1_given) this->AddInputFilename(mArgsInfo.input1_arg); + if (mArgsInfo.input2_given) { + mIsOperationUseASecondImage = true; + this->AddInputFilename(mArgsInfo.input2_arg); + if (mArgsInfo.operation_arg == 1) + mIsOutputScalar = true; + } + else if (mArgsInfo.operation_arg == 5 || mArgsInfo.operation_arg == 6) + mIsOutputScalar = true; + + if (mArgsInfo.output_given) this->SetOutputFilename(mArgsInfo.output_arg); + + // Check type of operation (with scalar or with other image) + if ((mArgsInfo.input2_given) && (mArgsInfo.scalar_given)) { + std::cerr << "ERROR : you cannot provide both --scalar and --input2 option" << std::endl; + exit(-1); + } + if ((!mArgsInfo.input2_given) && (!mArgsInfo.scalar_given)) { + if (mArgsInfo.operation_arg < 5) { + std::cerr << "Such operation need the --scalar option." << std::endl; + exit(-1); + } + } +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +template +void VectorArithmGenericFilter::UpdateWithInputImageType() +{ + // Read input1 + typename ImageType::Pointer input1 = this->template GetInput(0); + + // Set input image iterator + typedef itk::ImageRegionIterator IteratorType; + IteratorType it(input1, input1->GetLargestPossibleRegion()); + + // typedef input2 + typename ImageType::Pointer input2 = NULL; + IteratorType it2; + + // Special case for normalisation + /* + if (mTypeOfOperation == 12) { + typedef itk::MinimumMaximumImageCalculator MinMaxFilterType; + typename MinMaxFilterType::Pointer ff = MinMaxFilterType::New(); + ff->SetImage(input1); + ff->ComputeMaximum(); + mScalar = ff->GetMaximum(); + mTypeOfOperation = 11; // divide + } + */ + + if (mIsOperationUseASecondImage) { + // Read input2 + input2 = this->template GetInput(1); + // Set input image iterator + it2 = IteratorType(input2, input2->GetLargestPossibleRegion()); + // Check dimension + if (!clitk::HaveSameSize(input1, input2)) { + itkExceptionMacro(<< "The images (input and input2) must have the same size"); + } + if(!clitk::HaveSameSpacing(input1, input2)) { + itkWarningMacro(<< "The images (input and input2) do not have the same spacing. " + << "Using first input's information."); + } + } + + // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite + if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) { + // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is " + // << typeid(PixelType).name() + // << std::endl; + mOverwriteInputImage = false; + } + + // ---------------- Overwrite input Image --------------------- + if (mOverwriteInputImage && !mIsOutputScalar) { + // Set output iterator (to input1) + IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion()); + if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito); + else ComputeImage(it, ito); + this->template SetNextOutput(input1); + } + + // ---------------- Create new output Image --------------------- + else { + // Create output image + if (!mIsOutputScalar) { + typedef ImageType OutputImageType; + typename OutputImageType::Pointer output = OutputImageType::New(); + output->SetRegions(input1->GetLargestPossibleRegion()); + output->SetOrigin(input1->GetOrigin()); + output->SetSpacing(input1->GetSpacing()); + output->Allocate(); + // Set output iterator + typedef itk::ImageRegionIterator IteratorOutputType; + IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion()); + if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito); + else ComputeImage(it, ito); + this->template SetNextOutput(output); + } + else { + // Create scalar output image + typedef itk::Image OutputImageType; + typename OutputImageType::Pointer output = OutputImageType::New(); + output->SetRegions(input1->GetLargestPossibleRegion()); + output->SetOrigin(input1->GetOrigin()); + output->SetSpacing(input1->GetSpacing()); + output->Allocate(); + // Set output iterator + typedef itk::ImageRegionIterator IteratorOutputType; + IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion()); + if (mIsOperationUseASecondImage) ComputeScalarImage(it, it2, ito); + else ComputeScalarImage(it, ito); + this->template SetNextOutput(output); + } + } +} +//-------------------------------------------------------------------- + +//-------------------------------------------------------------------- +template +template +void VectorArithmGenericFilter::ComputeImage(Iter1 it1, Iter2 it2, Iter3 ito) +{ + it1.GoToBegin(); + it2.GoToBegin(); + ito.GoToBegin(); + typedef typename Iter3::PixelType PixelType; + + switch (mTypeOfOperation) { + case 0: // Addition + while (!ito.IsAtEnd()) { + ito.Set(it1.Get() + it2.Get()); + ++it1; + ++it2; + ++ito; + } + break; + /* + case 1: // Multiply + while (!ito.IsAtEnd()) { + ito.Set(it1.Get() * it2.Get()) ); + ++it1; + ++it2; + ++ito; + } + break; + */ + /* + case 2: // Divide + while (!ito.IsAtEnd()) { + if (it1.Get() != 0) + ito.Set(it1.Get() / it2.Get())); + else ito.Set(mDefaultPixelValue); + ++it1; + ++it2; + ++ito; + } + break; + case 3: // Max + while (!ito.IsAtEnd()) { + if (it1.Get() < it2.Get()) ito.Set(it2.Get()); + ++it1; + ++it2; + ++ito; + } + break; + case 4: // Min + while (!ito.IsAtEnd()) { + if (it1.Get() > it2.Get()) ito.Set(it2.Get()); + ++it1; + ++it2; + ++ito; + } + break; + */ + /* + case 5: // Absolute difference + while (!ito.IsAtEnd()) { + ito.Set(it2.Get()-it1.Get()); + ++it1; + ++it2; + ++ito; + } + break; + */ + /* + case 6: // Squared differences + while (!ito.IsAtEnd()) { + ito.Set(pow(it1.Get()-it2.Get(),2))); + ++it1; + ++it2; + ++ito; + } + break; + */ + case 7: // Difference + while (!ito.IsAtEnd()) { + ito.Set(it1.Get()-it2.Get()); + ++it1; + ++it2; + ++ito; + } + break; + /* + case 8: // Relative Difference + while (!ito.IsAtEnd()) { + if (it1.Get() != 0) ito.Set(PixelTypeDownCast(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get()); + else ito.Set(0.0); + ++it1; + ++it2; + ++ito; + } + break; + */ + /* + case 9: // CrossProduct + while (!ito.IsAtEnd()) { + ito.Set(it1.Get()^it2.Get()); + ++it1; + ++it2; + ++ito; + } + break; + */ + default: // error ? + std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl; + exit(-1); + } +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +template +void clitk::VectorArithmGenericFilter::ComputeImage(Iter1 it, Iter2 ito) +{ + ito.GoToBegin(); + it.GoToBegin(); + + typedef typename Iter1::PixelType PixelType; + + PixelType scalar_vector; + scalar_vector.Fill(mScalar); + + // Perform operation + switch (mTypeOfOperation) { + case 0: // Addition + while (!it.IsAtEnd()) { + ito.Set(it.Get() + scalar_vector); + ++it; + ++ito; + } + break; + case 1: // Multiply + while (!it.IsAtEnd()) { + ito.Set(it.Get() * mScalar); + ++it; + ++ito; + } + break; + /* + case 2: // Inverse + while (!it.IsAtEnd()) { + if (it.Get() != 0) + ito.Set(mScalar / it.Get())); + else ito.Set(mDefaultPixelValue); + ++it; + ++ito; + } + break; + case 3: // Max + while (!it.IsAtEnd()) { + if (it.Get() < mScalar) ito.Set(PixelTypeDownCast(mScalar)); + else ito.Set(PixelTypeDownCast(it.Get())); + ++it; + ++ito; + } + break; + case 4: // Min + while (!it.IsAtEnd()) { + if (it.Get() > mScalar) ito.Set(PixelTypeDownCast(mScalar)); + else ito.Set(PixelTypeDownCast(it.Get())); + ++it; + ++ito; + } + break; + */ + /* + case 5: // Absolute value + while (!it.IsAtEnd()) { + ito.Set(PixelTypeDownCast(it.GetNorm())); + ++it; + ++ito; + } + break; + case 6: // Squared value + while (!it.IsAtEnd()) { + ito.Set(PixelTypeDownCast(it.GetSquaredNorm()); + ++it; + ++ito; + } + break; + */ + /* + case 7: // Log + while (!it.IsAtEnd()) { + if (it.Get() > 0) + ito.Set(PixelTypeDownCast(log((double)it.Get()))); + else ito.Set(mDefaultPixelValue); + ++it; + ++ito; + } + break; + case 8: // exp + while (!it.IsAtEnd()) { + ito.Set(PixelTypeDownCast(exp((double)it.Get()))); + ++it; + ++ito; + } + break; + case 9: // sqrt + while (!it.IsAtEnd()) { + if (it.Get() > 0) + ito.Set(PixelTypeDownCast(sqrt((double)it.Get()))); + else { + if (it.Get() ==0) ito.Set(0); + else ito.Set(mDefaultPixelValue); + } + ++it; + ++ito; + } + break; + case 10: // exp + while (!it.IsAtEnd()) { + ito.Set(PixelTypeDownCast((0x10000 - (double)it.Get())/mScalar)); + ++it; + ++ito; + } + break; + */ + case 11: // divide + while (!it.IsAtEnd()) { + ito.Set(it.Get() / mScalar); + ++it; + ++ito; + } + break; + case 12: // normalize + while (!it.IsAtEnd()) { + PixelType n = it.Get(); + n.Normalize(); + ito.Set(n); + ++it; + ++ito; + } + break; + default: // error ? + std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl; + exit(-1); + } +} +//-------------------------------------------------------------------- + +//-------------------------------------------------------------------- +template +template +void VectorArithmGenericFilter::ComputeScalarImage(Iter1 it1, Iter2 it2, Iter3 ito) +{ + it1.GoToBegin(); + it2.GoToBegin(); + ito.GoToBegin(); + typedef typename Iter3::PixelType PixelType; + + switch (mTypeOfOperation) { + case 1: // Multiply + while (!ito.IsAtEnd()) { + ito.Set(it1.Get() * it2.Get()); + ++it1; + ++it2; + ++ito; + } + break; + default: // error ? + std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl; + exit(-1); + } +} +//-------------------------------------------------------------------- + +//-------------------------------------------------------------------- +template +template +void clitk::VectorArithmGenericFilter::ComputeScalarImage(Iter1 it, Iter2 ito) +{ + ito.GoToBegin(); + it.GoToBegin(); + + typedef typename Iter2::PixelType PixelType; + + + // Perform operation + switch (mTypeOfOperation) { + case 5: // Absolute value + while (!it.IsAtEnd()) { + ito.Set(PixelTypeDownCast(it.Get().GetNorm())); + ++it; + ++ito; + } + break; + case 6: // Squared value + while (!it.IsAtEnd()) { + ito.Set(PixelTypeDownCast(it.Get().GetSquaredNorm())); + ++it; + ++ito; + } + break; + default: // error ? + std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl; + exit(-1); + } +} +//-------------------------------------------------------------------- + + +} // end namespace + +#endif //#define CLITKVectorArithmGENERICFILTER_TXX