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})
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} )
#include <gdcmReader.h>
#endif
+#include <set>
+
//====================================================================
int main(int argc, char * argv[])
{
//===========================================
/// Get slices locations ...
- std::vector<double> theorigin(3);
- std::vector<double> sliceLocations;
+ int series_number = -1;
+ std::set<int> series_numbers;
+ std::map< int, std::vector<double> > theorigin;
+ std::map< int, std::vector<double> > sliceLocations;
+ std::map< int, std::vector<std::string> > seriesFiles;
for(unsigned int i=0; i<args_info.inputs_num; i++) {
//std::cout << "Reading <" << input_files[i] << std::endl;
#if GDCM_MAJOR_VERSION == 2
gdcm::Reader hreader;
hreader.SetFileName(input_files[i].c_str());
hreader.Read();
- theorigin = gdcm::ImageHelper::GetOriginValue(hreader.GetFile());
- sliceLocations.push_back(theorigin[2]);
- gdcm::Attribute<0x28, 0x100> 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)
{
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;
//===========================================
// Sort slices locations ...
- std::vector<int> 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.size(); i++) {
- std::cout << sliceLocations[sliceIndex[i]] << " -> "
- << sliceIndex[i] << " / " << i << " => "
- << sliceLocations[sliceIndex[i]] - sliceLocations[sliceIndex[i-1]]
- << "mm "
- << input_files[sliceIndex[i]]
+ std::set<int>::iterator sn = series_numbers.begin();
+ while ( sn != series_numbers.end() ) {
+ std::vector<double> locs = sliceLocations[*sn];
+ std::vector<double> origin = theorigin[*sn];
+ std::vector<std::string> files = seriesFiles[*sn];
+ std::vector<int> 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<sliceIndex.size(); i++) {
- currentDist = sliceLocations[sliceIndex[i]]-previous;
- if (i!=1) {
- if (fabs(dist-currentDist) > 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.size(); i++) {
+ std::cout << locs[sliceIndex[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<std::string> sorted_files;
- sorted_files.resize(sliceIndex.size());
- for(unsigned int i=0; i<sliceIndex.size(); i++)
- sorted_files[i] = input_files[ sliceIndex[i] ];
+ //===========================================
+ // Analyze slices locations ...
+ double currentDist;
+ double dist=0;
+ double tolerance = args_info.tolerance_arg;
+ double previous = locs[sliceIndex[0]];
+ for(unsigned int i=1; i<sliceIndex.size(); i++) {
+ currentDist = locs[sliceIndex[i]]-previous;
+ if (i!=1) {
+ if (fabs(dist-currentDist) > 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<double> spacing = image->GetSpacing();
- std::vector<int> 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<std::string> sorted_files;
+ sorted_files.resize(sliceIndex.size());
+ for(unsigned int i=0; i<sliceIndex.size(); i++)
+ sorted_files[i] = files[ sliceIndex[i] ];
- vvImageWriter::Pointer writer = vvImageWriter::New();
- writer->SetInput(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<double> spacing = image->GetSpacing();
+ std::vector<int> 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;
}
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
// template<>
// class ImageArithmGenericFilter<args_info_clitkImageArithm>;
-////////////////////////////////////////////////////////////////////
- // specializations for itk::Vector<float, 3u>, 3u
- template<>
- template<>
- void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<float, 3u>, 3u > >()
- {
- typedef itk::Image< itk::Vector<float, 3u>, 3u > ImageType;
-
- // Read input1
- ImageType::Pointer input1 = this->GetInput<ImageType>(0);
-
- // Set input image iterator
- typedef itk::ImageRegionIterator<ImageType> IteratorType;
- IteratorType it(input1, input1->GetLargestPossibleRegion());
-
- // typedef input2
- ImageType::Pointer input2 = NULL;
- IteratorType it2;
-
- /*
- // Special case for normalisation
- if (mTypeOfOperation == 12) {
- typedef itk::MinimumMaximumImageCalculator<ImageType> 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<ImageType>(1);
- // Set input image iterator
- it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
- // Check dimension
- if (!clitk::HaveSameSize<ImageType, ImageType>(input1, input2)) {
- itkExceptionMacro(<< "The images (input and input2) must have the same size");
- }
- if(!clitk::HaveSameSpacing<ImageType, ImageType>(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<ImageType>(input1);
- }
- // ---------------- Create new output Image ---------------------
- else {
- /*if (mOutputIsFloat) {
- // Create output image
- typedef itk::Image<float,ImageType::ImageDimension> 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<OutputImageType> IteratorOutputType;
- IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
- if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
- else ComputeImage(it, ito);
- this->template SetNextOutput<OutputImageType>(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<OutputImageType> IteratorOutputType;
- IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
- if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
- else ComputeImage(it, ito);
- this->SetNextOutput<OutputImageType>(output);
- }
- }
- }
-
- template<>
- template<>
- void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
- itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
- itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > >
- (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it,
- itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito)
- {
- typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 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<double, PixelType>(mScalar));
- else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
- ++it;
- ++ito;
- }
- break;
- case 4: // Min
- while (!it.IsAtEnd()) {
- if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
- else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
- ++it;
- ++ito;
- }
- break;
- case 5: // Absolute value
- while (!it.IsAtEnd()) {
- if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
- // <= zero to avoid warning for unsigned types
- else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
- ++it;
- ++ito;
- }
- break;
- case 6: // Squared value
- while (!it.IsAtEnd()) {
- ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
- ++it;
- ++ito;
- }
- break;
- case 7: // Log
- while (!it.IsAtEnd()) {
- if (it.Get() > 0)
- ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
- else ito.Set(mDefaultPixelValue);
- ++it;
- ++ito;
- }
- break;
- case 8: // exp
- while (!it.IsAtEnd()) {
- ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
- ++it;
- ++ito;
- }
- break;
- case 9: // sqrt
- while (!it.IsAtEnd()) {
- if (it.Get() > 0)
- ito.Set(PixelTypeDownCast<double, PixelType>(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<double, PixelType>((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<args_info_clitkImageArithm>::ComputeImage<
- itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
- itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
- itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >
- >
- (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it1,
- itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it2,
- itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito)
- {
- typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter1;
- typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter2;
- typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 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, PixelType>(((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<double, 3u>, 3u
- template<>
- template<>
- void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<double, 3u>, 3u > >()
- {
- typedef itk::Image< itk::Vector<double, 3u>, 3u > ImageType;
-
- // Read input1
- ImageType::Pointer input1 = this->GetInput<ImageType>(0);
-
- // Set input image iterator
- typedef itk::ImageRegionIterator<ImageType> IteratorType;
- IteratorType it(input1, input1->GetLargestPossibleRegion());
-
- // typedef input2
- ImageType::Pointer input2 = NULL;
- IteratorType it2;
-
- /*
- // Special case for normalisation
- if (mTypeOfOperation == 12) {
- typedef itk::MinimumMaximumImageCalculator<ImageType> 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<ImageType>(1);
- // Set input image iterator
- it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
- // Check dimension
- if (!clitk::HaveSameSize<ImageType, ImageType>(input1, input2)) {
- itkExceptionMacro(<< "The images (input and input2) must have the same size");
- }
- if(!clitk::HaveSameSpacing<ImageType, ImageType>(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<ImageType>(input1);
- }
- // ---------------- Create new output Image ---------------------
- else {
- /*if (mOutputIsFloat) {
- // Create output image
- typedef itk::Image<float,ImageType::ImageDimension> 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<OutputImageType> IteratorOutputType;
- IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
- if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
- else ComputeImage(it, ito);
- this->template SetNextOutput<OutputImageType>(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<OutputImageType> IteratorOutputType;
- IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
- if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
- else ComputeImage(it, ito);
- this->SetNextOutput<OutputImageType>(output);
- }
- }
- }
-
- template<>
- template<>
- void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
- itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >,
- itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > >
- (itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it,
- itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > ito)
- {
- typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 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<double, PixelType>(mScalar));
- else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
- ++it;
- ++ito;
- }
- break;
- case 4: // Min
- while (!it.IsAtEnd()) {
- if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
- else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
- ++it;
- ++ito;
- }
- break;
- case 5: // Absolute value
- while (!it.IsAtEnd()) {
- if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
- // <= zero to avoid warning for unsigned types
- else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
- ++it;
- ++ito;
- }
- break;
- case 6: // Squared value
- while (!it.IsAtEnd()) {
- ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
- ++it;
- ++ito;
- }
- break;
- case 7: // Log
- while (!it.IsAtEnd()) {
- if (it.Get() > 0)
- ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
- else ito.Set(mDefaultPixelValue);
- ++it;
- ++ito;
- }
- break;
- case 8: // exp
- while (!it.IsAtEnd()) {
- ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
- ++it;
- ++ito;
- }
- break;
- case 9: // sqrt
- while (!it.IsAtEnd()) {
- if (it.Get() > 0)
- ito.Set(PixelTypeDownCast<double, PixelType>(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<double, PixelType>((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<args_info_clitkImageArithm>::ComputeImage<
- itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >,
- itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >,
- itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >
- >
- (itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it1,
- itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it2,
- itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > ito)
- {
- typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > Iter1;
- typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > Iter2;
- typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 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, PixelType>(((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);
- }
- }
}
}; // end class ImageArithmGenericFilter
- // specializations for itk::Vector<float, 3u>, 3u
- template<> template<>
- void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<float, 3u>, 3u > >();
-
- template<> template<>
- void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
- itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
- itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >
- >
- (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it,
- itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito);
-
- template<> template<>
- void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
- itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
- itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
- itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >
- >
- (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it1,
- itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it2,
- itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito);
-
- // specializations for itk::Vector<double, 3u>, 3u
- template<> template<>
- void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<double, 3u>, 3u > >();
-
- template<> template<>
- void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
- itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >,
- itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >
- >
- (itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it,
- itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > ito);
-
- template<> template<>
- void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
- itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >,
- itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >,
- itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >
- >
- (itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it1,
- itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it2,
- itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > ito);
} // end namespace
//--------------------------------------------------------------------
void ImageArithmGenericFilter<args_info_type>::InitializeImageType()
{
ADD_DEFAULT_IMAGE_TYPES(Dim);
- ADD_VEC_IMAGE_TYPE(3u,3u,float);
- ADD_VEC_IMAGE_TYPE(3u,3u,double);
}
//--------------------------------------------------------------------
--- /dev/null
+/*=========================================================================
+ 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 <David.Sarrut@creatis.insa-lyon.fr>
+ * @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<args_info_clitkVectorArithm> 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
--- /dev/null
+#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
--- /dev/null
+/*=========================================================================
+ 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<args_info_clitkVECTORARITHM>;
+
+
+}
+
+#endif //define CLITKVECTORARITHMGENERICFILTER_CXX
--- /dev/null
+/*=========================================================================
+ 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 <David.Sarrut@creatis.insa-lyon.fr>
+ * @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 args_info_type>
+ class ITK_EXPORT VectorArithmGenericFilter:
+ public clitk::ImageToImageGenericFilter<VectorArithmGenericFilter<args_info_type> > {
+
+ public:
+
+ // Constructor
+ VectorArithmGenericFilter ();
+
+ // Types
+ typedef VectorArithmGenericFilter Self;
+ typedef ImageToImageGenericFilterBase Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> 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<class InputImageType>
+ void UpdateWithInputImageType();
+
+ protected:
+ template<unsigned int Dim> void InitializeImageType();
+ bool mIsOperationUseASecondImage;
+ double mScalar;
+ double mDefaultPixelValue;
+ int mTypeOfOperation;
+ args_info_type mArgsInfo;
+ bool mOverwriteInputImage;
+ bool mOutputIsFloat;
+ bool mIsOutputScalar;
+
+ template<class Iter1, class Iter2>
+ void ComputeImage(Iter1 it, Iter2 ito);
+
+ template<class Iter1, class Iter2, class Iter3>
+ void ComputeImage(Iter1 it1, Iter2 it2, Iter3 ito);
+
+ template<class Iter1, class Iter2>
+ void ComputeScalarImage(Iter1 it, Iter2 ito);
+
+ template<class Iter1, class Iter2, class Iter3>
+ void ComputeScalarImage(Iter1 it1, Iter2 it2, Iter3 ito);
+
+ //--------------------------------------------------------------------
+
+ }; // end class VectorArithmGenericFilter
+
+ // specializations for itk::Vector<float, 3u>, 3u
+} // end namespace
+//--------------------------------------------------------------------
+
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "clitkVectorArithmGenericFilter.txx"
+#endif
+
+#endif //#define CLITKVectorArithmGENERICFILTER_H
+
--- /dev/null
+/*=========================================================================
+ 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<class args_info_type>
+VectorArithmGenericFilter<args_info_type>::VectorArithmGenericFilter()
+ :ImageToImageGenericFilter<Self>("VectorArithmGenericFilter"),mTypeOfOperation(0)
+{
+ InitializeImageType<3>();
+ mIsOperationUseASecondImage = false;
+ mIsOutputScalar = false;
+ mOverwriteInputImage = true;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<class args_info_type>
+template<unsigned int Dim>
+void VectorArithmGenericFilter<args_info_type>::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<class args_info_type>
+void VectorArithmGenericFilter<args_info_type>::EnableOverwriteInputImage(bool b)
+{
+ mOverwriteInputImage = b;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<class args_info_type>
+void VectorArithmGenericFilter<args_info_type>::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<class args_info_type>
+template<class ImageType>
+void VectorArithmGenericFilter<args_info_type>::UpdateWithInputImageType()
+{
+ // Read input1
+ typename ImageType::Pointer input1 = this->template GetInput<ImageType>(0);
+
+ // Set input image iterator
+ typedef itk::ImageRegionIterator<ImageType> 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<ImageType> 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<ImageType>(1);
+ // Set input image iterator
+ it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
+ // Check dimension
+ if (!clitk::HaveSameSize<ImageType, ImageType>(input1, input2)) {
+ itkExceptionMacro(<< "The images (input and input2) must have the same size");
+ }
+ if(!clitk::HaveSameSpacing<ImageType, ImageType>(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<ImageType>(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<OutputImageType> IteratorOutputType;
+ IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
+ if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
+ else ComputeImage(it, ito);
+ this->template SetNextOutput<OutputImageType>(output);
+ }
+ else {
+ // Create scalar output image
+ typedef itk::Image<typename ImageType::PixelType::ValueType, ImageType::ImageDimension> 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<OutputImageType> IteratorOutputType;
+ IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
+ if (mIsOperationUseASecondImage) ComputeScalarImage(it, it2, ito);
+ else ComputeScalarImage(it, ito);
+ this->template SetNextOutput<OutputImageType>(output);
+ }
+ }
+}
+//--------------------------------------------------------------------
+
+//--------------------------------------------------------------------
+template<class args_info_type>
+template<class Iter1, class Iter2, class Iter3>
+void VectorArithmGenericFilter<args_info_type>::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, PixelType>(((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<class args_info_type>
+template<class Iter1, class Iter2>
+void clitk::VectorArithmGenericFilter<args_info_type>::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<double, PixelType>(mScalar));
+ else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
+ ++it;
+ ++ito;
+ }
+ break;
+ case 4: // Min
+ while (!it.IsAtEnd()) {
+ if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
+ else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
+ ++it;
+ ++ito;
+ }
+ break;
+ */
+ /*
+ case 5: // Absolute value
+ while (!it.IsAtEnd()) {
+ ito.Set(PixelTypeDownCast<double, PixelType>(it.GetNorm()));
+ ++it;
+ ++ito;
+ }
+ break;
+ case 6: // Squared value
+ while (!it.IsAtEnd()) {
+ ito.Set(PixelTypeDownCast<double, PixelType>(it.GetSquaredNorm());
+ ++it;
+ ++ito;
+ }
+ break;
+ */
+ /*
+ case 7: // Log
+ while (!it.IsAtEnd()) {
+ if (it.Get() > 0)
+ ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
+ else ito.Set(mDefaultPixelValue);
+ ++it;
+ ++ito;
+ }
+ break;
+ case 8: // exp
+ while (!it.IsAtEnd()) {
+ ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
+ ++it;
+ ++ito;
+ }
+ break;
+ case 9: // sqrt
+ while (!it.IsAtEnd()) {
+ if (it.Get() > 0)
+ ito.Set(PixelTypeDownCast<double, PixelType>(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<double, PixelType>((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<class args_info_type>
+template<class Iter1, class Iter2, class Iter3>
+void VectorArithmGenericFilter<args_info_type>::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<class args_info_type>
+template<class Iter1, class Iter2>
+void clitk::VectorArithmGenericFilter<args_info_type>::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<double, PixelType>(it.Get().GetNorm()));
+ ++it;
+ ++ito;
+ }
+ break;
+ case 6: // Squared value
+ while (!it.IsAtEnd()) {
+ ito.Set(PixelTypeDownCast<double, PixelType>(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
git://git.creatis.insa-lyon.fr / clitk.git / commitdiff