ITK/Examples/EdgesAndGradients/SobelEdgeDetectionImageFilter: Difference between revisions

From KitwarePublic
< ITK‎ | Examples
Jump to navigationJump to search
(Added Quickview and user-specified file.)
(Added Quickview and user-specified file.)
Line 1: Line 1:
==SobelEdgeDetectionImageFilter.cxx==
==SobelEdgeDetectionImageFilter.cxx==
<source lang="cpp">
<source lang="cpp">
/*=========================================================================
#include "itkImage.h"
*
#include "itkSobelEdgeDetectionImageFilter.h"
*  Copyright Insight Software Consortium
#include "itkImageFileReader.h"
*
*  Licensed under the Apache License, Version 2.0 (the "License");
*  you may not use this file except in compliance with the License.
*  You may obtain a copy of the License at
*
*        http://www.apache.org/licenses/LICENSE-2.0.txt
*
*  Unless required by applicable law or agreed to in writing, software
*  distributed under the License is distributed on an "AS IS" BASIS,
*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
*  See the License for the specific language governing permissions and
*  limitations under the License.
*
*=========================================================================*/


#include "itkImageRegistrationMethod.h"
#include "itksys/SystemTools.hxx"
#include "itkAffineTransform.h"
#include "itkMutualInformationImageToImageMetric.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkGradientDescentOptimizer.h"


#include "itkTextOutput.h"
#include "QuickView.h"
#include "itkImageRegionIterator.h"


namespace
typedef itk::Image<unsigned char, 2>  UnsignedCharImageType;
{
typedef itk::Image<float, 2>         FloatImageType;
 
double F( itk::Vector<double,3> & v );
}


#if ITK_VERSION_MAJOR >= 4
typedef UnsignedCharImageType ImageType;
#else
typedef FloatImageType ImageType;
#endif
static void CreateImage(ImageType::Pointer);


/**
int main(int argc, char *argv[])
*  This program test one instantiation of the itk::ImageRegistrationMethod class
*
*  This file tests the combination of:
*  - MutualInformation
*  - AffineTransform
*  - GradientDescentOptimizer
*  - LinearInterpolateImageFunction
*
*  The test image pattern consists of a 3D gaussian in the middle
*  with some directional pattern on the outside.
*  One image is scaled and shifted relative to the other.
*
* Notes:
* =======
* This example performs an affine registration
* between a moving (source) and fixed (target) image using mutual information.
* It uses the optimization method of Viola and Wells to find the
* best affine transform to register the moving image onto the fixed
* image.
*
* The mutual information value and its derivatives are estimated
* using spatial sampling. The performance
* of the registration depends on good choices of the parameters
* used to estimate the mutual information. Refer to the documentation
* for MutualInformationImageToImageMetric for details on these
* parameters and how to set them.
*
* The registration uses a simple stochastic gradient ascent scheme. Steps
* are repeatedly taken that are proportional to the approximate
* deriviative of the mutual information with respect to the affine
* transform parameters. The stepsize is governed by the LearningRate
* parameter.
*
* Since the parameters of the linear part is different in magnitude
* to the parameters in the offset part, scaling is required
* to improve convergence. The scaling can set via the optimizer.
*
* NB: In the Viola and Wells paper, the scaling is specified by
* using different learning rates for the linear and offset part.
* The following formula translate their scaling parameters to
* those used in this framework:
*
* LearningRate = lambda_R
* TranslationScale = sqrt( lambda_T / lambda_R );
*
* In the optimizer's scale transform set the scaling for
* all the translation parameters to TranslationScale^{-2}.
* Set the scale for all other parameters to 1.0.
*
* Note: the optimization performance can be improved by
* setting the image origin to center of mass of the image.
*
* Implementaton of this example and related components are based on:
* Viola, P. and Wells III, W. (1997).
* "Alignment by Maximization of Mutual Information"
* International Journal of Computer Vision, 24(2):137-154
*
*/
int itkImageRegistrationMethodTest_13(int, char* [] )
{
{
 
   ImageType::Pointer image = ImageType::New();
  itk::OutputWindow::SetInstance(itk::TextOutput::New().GetPointer());
   if (argc < 2)
 
  bool pass = true;
 
  const unsigned int dimension = 3;
  unsigned int j;
 
  typedef float  PixelType;
 
  // Fixed Image Type
   typedef itk::Image<PixelType,dimension>              FixedImageType;
 
  // Moving Image Type
  typedef itk::Image<PixelType,dimension>              MovingImageType;
 
  // Transform Type
  typedef itk::AffineTransform< double,dimension >  TransformType;
 
  // Optimizer Type
  typedef itk::GradientDescentOptimizer            OptimizerType;
 
  // Metric Type
  typedef itk::MutualInformationImageToImageMetric<
                                    FixedImageType,
                                    MovingImageType >    MetricType;
 
  // Interpolation technique
  typedef itk:: LinearInterpolateImageFunction<
                                    MovingImageType,
                                    double          >    InterpolatorType;
 
  // Registration Method
  typedef itk::ImageRegistrationMethod<
                                    FixedImageType,
                                    MovingImageType >    RegistrationType;
 
 
  MetricType::Pointer         metric        = MetricType::New();
   TransformType::Pointer      transform    = TransformType::New();
  OptimizerType::Pointer      optimizer    = OptimizerType::New();
  FixedImageType::Pointer    fixedImage    = FixedImageType::New();
  MovingImageType::Pointer    movingImage  = MovingImageType::New();
  InterpolatorType::Pointer  interpolator  = InterpolatorType::New();
  RegistrationType::Pointer  registration  = RegistrationType::New();
 
  /*********************************************************
  * Set up the two input images.
  * One image scaled and shifted with respect to the other.
  **********************************************************/
  double displacement[dimension] = {7,3,2};
  double scale[dimension] = { 0.80, 1.0, 1.0 };
 
  FixedImageType::SizeType size = {{100,100,40}};
  FixedImageType::IndexType index = {{0,0,0}};
  FixedImageType::RegionType region;
  region.SetSize( size );
  region.SetIndex( index );
 
  fixedImage->SetLargestPossibleRegion( region );
  fixedImage->SetBufferedRegion( region );
  fixedImage->SetRequestedRegion( region );
  fixedImage->Allocate();
 
  movingImage->SetLargestPossibleRegion( region );
  movingImage->SetBufferedRegion( region );
  movingImage->SetRequestedRegion( region );
  movingImage->Allocate();
 
 
  typedef itk::ImageRegionIterator<MovingImageType> MovingImageIterator;
  typedef itk::ImageRegionIterator<FixedImageType>  FixedImageIterator;
 
  itk::Point<double,dimension> center;
  for ( j = 0; j < dimension; j++ )
     {
     {
     center[j] = 0.5 *  (double)region.GetSize()[j];
     CreateImage(image);
     }
     }
 
   else
   itk::Point<double,dimension> p;
  itk::Vector<double,dimension> d;
 
  MovingImageIterator mIter( movingImage, region );
  FixedImageIterator  fIter( fixedImage, region );
 
  while( !mIter.IsAtEnd() )
     {
     {
     for ( j = 0; j < dimension; j++ )
     typedef itk::ImageFileReader<ImageType> ReaderType;
      {
    ReaderType::Pointer reader = ReaderType::New();
      p[j] = mIter.GetIndex()[j];
     reader->SetFileName (argv[1]);
      }
     reader->Update();
 
     image = reader->GetOutput();
     d = p - center;
     }
 
    fIter.Set( (PixelType) F(d) );
 
     for ( j = 0; j < dimension; j++ )
      {
      d[j] = d[j] * scale[j] + displacement[j];
      }
 
     mIter.Set( (PixelType) F(d) );
 
     ++fIter;
    ++mIter;


    }
  typedef itk::SobelEdgeDetectionImageFilter <ImageType, FloatImageType>
          SobelEdgeDetectionImageFilterType;
  SobelEdgeDetectionImageFilterType::Pointer sobelFilter
          = SobelEdgeDetectionImageFilterType::New();
  sobelFilter->SetInput(image);


   // set the image origin to be center of the image
   QuickView viewer;
  double transCenter[dimension];
   if (argc > 1)
   for ( j = 0; j < dimension; j++ )
     {
     {
     transCenter[j] = -0.5 * double(size[j]);
     viewer.AddImage<ImageType>(
      image,
      true,
      itksys::SystemTools::GetFilenameName(argv[1]));
     }
     }
 
   else
   movingImage->SetOrigin( transCenter );
  fixedImage->SetOrigin( transCenter );
 
 
  /******************************************************************
  * Set up the optimizer.
  ******************************************************************/
 
  // set the translation scale
  typedef OptimizerType::ScalesType ScalesType;
  ScalesType parametersScales( transform->GetNumberOfParameters() );
 
  parametersScales.Fill( 1.0 );
 
  for ( j = 9; j < 12; j++ )
     {
     {
     parametersScales[j] = 0.001;
     viewer.AddImage<ImageType>(image);
     }
     }
  viewer.AddImage<FloatImageType>(
    sobelFilter->GetOutput(),
    true,
    "SobelEdgeDetectionImageFilter");


   optimizer->SetScales( parametersScales );
   viewer.Visualize();


   // need to maximize for mutual information
   return EXIT_SUCCESS;
  optimizer->MaximizeOn();
}


  /******************************************************************
void CreateImage(ImageType::Pointer image)
  * Set up the metric.
{
  ******************************************************************/
   itk::Index<2> start;
  metric->SetMovingImageStandardDeviation( 5.0 );
   start.Fill(0);
  metric->SetFixedImageStandardDeviation( 5.0 );
  metric->SetNumberOfSpatialSamples( 100 );
  metric->SetFixedImageRegion( fixedImage->GetBufferedRegion() );
 
 
   /******************************************************************
  * Set up the registrator.
  ******************************************************************/
 
  // connect up the components
  registration->SetMetric( metric );
  registration->SetOptimizer( optimizer );
  registration->SetTransform( transform );
  registration->SetFixedImage( fixedImage );
  registration->SetMovingImage( movingImage );
  registration->SetInterpolator( interpolator );
 
  // set initial parameters to identity
  RegistrationType::ParametersType initialParameters(
    transform->GetNumberOfParameters() );
 
   initialParameters.Fill( 0.0 );
  initialParameters[0] = 1.0;
  initialParameters[4] = 1.0;
  initialParameters[8] = 1.0;


  itk::Size<2> size;
  size.Fill(100);


   /***********************************************************
   itk::ImageRegion<2> region(start, size);
  * Run the registration - reducing learning rate as we go
   image->SetRegions(region);
  ************************************************************/
   image->Allocate();
  const unsigned int numberOfLoops = 3;
   unsigned int iter[numberOfLoops] = { 300, 300, 350 };
   double      rates[numberOfLoops] = { 1e-3, 5e-4, 1e-4 };


   for ( j = 0; j < numberOfLoops; j++ )
  // Make a square
   for(unsigned int r = 20; r < 80; r++)
     {
     {
 
     for(unsigned int c = 20; c < 80; c++)
     try
      {
      optimizer->SetNumberOfIterations( iter[j] );
      optimizer->SetLearningRate( rates[j] );
      registration->SetInitialTransformParameters( initialParameters );
      registration->Update();
 
      initialParameters = registration->GetLastTransformParameters();
 
      }
    catch( itk::ExceptionObject & e )
       {
       {
       std::cout << "Registration failed" << std::endl;
       itk::Index<2> pixelIndex;
      std::cout << "Reason " << e.GetDescription() << std::endl;
       pixelIndex[0] = r;
      return EXIT_FAILURE;
      pixelIndex[1] = c;
       }
 
    }
 
 
  /***********************************************************
  * Check the results
  ************************************************************/
  RegistrationType::ParametersType solution =
    registration->GetLastTransformParameters();
 
  std::cout << "Solution is: " << solution << std::endl;
 
 
  RegistrationType::ParametersType trueParameters(
    transform->GetNumberOfParameters() );
  trueParameters.Fill( 0.0 );
  trueParameters[ 0] = 1/scale[0];
  trueParameters[ 4] = 1/scale[1];
  trueParameters[ 8] = 1/scale[2];
  trueParameters[ 9] = - displacement[0]/scale[0];
  trueParameters[10] = - displacement[1]/scale[1];
  trueParameters[11] = - displacement[2]/scale[2];
 
  std::cout << "True solution is: " << trueParameters << std::endl;


  for( j = 0; j < 9; j++ )
      image->SetPixel(pixelIndex, 15);
    {
    if( vnl_math_abs( solution[j] - trueParameters[j] ) > 0.025 )
      {
      std::cout << "Failed " << j << std::endl;
      pass = false;
      }
    }
  for( j = 9; j < 12; j++ )
    {
    if( vnl_math_abs( solution[j] - trueParameters[j] ) > 1.0 )
      {
      std::cout << "Failed " << j << std::endl;
      pass = false;
       }
       }
     }
     }
  if( !pass )
    {
    std::cout << "Test failed." << std::endl;
    return EXIT_FAILURE;
    }
  /*************************************************
  * Check for parzen window exception
  **************************************************/
  double oldValue = metric->GetMovingImageStandardDeviation();
  metric->SetMovingImageStandardDeviation( 0.005 );
  try
    {
    pass = false;
    registration->Update();
    }
  catch(itk::ExceptionObject& err)
    {
    std::cout << "Caught expected ExceptionObject" << std::endl;
    std::cout << err << std::endl;
    pass = true;
    }
  if( !pass )
    {
    std::cout << "Should have caught an exception" << std::endl;
    std::cout << "Test failed." << std::endl;
    return EXIT_FAILURE;
    }
  metric->SetMovingImageStandardDeviation( oldValue );
  /*************************************************
  * Check for mapped out of image error
  **************************************************/
  solution[5] = 1000;
  registration->SetInitialTransformParameters( solution );
  try
    {
    pass = false;
    registration->Update();
    }
  catch(itk::ExceptionObject& err)
    {
    std::cout << "Caught expected ExceptionObject" << std::endl;
    std::cout << err << std::endl;
    pass = true;
    }
  if( !pass )
    {
    std::cout << "Should have caught an exception" << std::endl;
    std::cout << "Test failed." << std::endl;
    return EXIT_FAILURE;
    }
  std::cout << "Test passed." << std::endl;
  return EXIT_SUCCESS;
}
namespace
{
/**
* This function defines the test image pattern.
* The pattern is a 3D gaussian in the middle
* and some directional pattern on the outside.
*/
double F( itk::Vector<double,3> & v )
{
  double x = v[0];
  double y = v[1];
  double z = v[2];
  const double s = 50;
  double value = 200.0 * vcl_exp( - ( x*x + y*y + z*z )/(s*s) );
  x -= 8; y += 3; z += 0;
  double r = vcl_sqrt( x*x + y*y + z*z );
  if( r > 35 )
    {
    value = 2 * ( vnl_math_abs( x ) +
      0.8 * vnl_math_abs( y ) +
      0.5 * vnl_math_abs( z ) );
    }
  if( r < 4 )
    {
    value = 400;
    }
  return value;
}
}
}
</source>
</source>


{{ITKCMakeLists|SobelEdgeDetectionImageFilter}}
{{ITKVTKCMakeLists|SobelEdgeDetectionImageFilter}}

Revision as of 17:38, 2 October 2011

SobelEdgeDetectionImageFilter.cxx

<source lang="cpp">

  1. include "itkImage.h"
  2. include "itkSobelEdgeDetectionImageFilter.h"
  3. include "itkImageFileReader.h"
  1. include "itksys/SystemTools.hxx"
  1. include "QuickView.h"

typedef itk::Image<unsigned char, 2> UnsignedCharImageType; typedef itk::Image<float, 2> FloatImageType;

  1. if ITK_VERSION_MAJOR >= 4

typedef UnsignedCharImageType ImageType;

  1. else

typedef FloatImageType ImageType;

  1. endif

static void CreateImage(ImageType::Pointer);

int main(int argc, char *argv[]) {

 ImageType::Pointer image = ImageType::New();
 if (argc < 2)
   {
   CreateImage(image);
   }
 else
   {
   typedef itk::ImageFileReader<ImageType> ReaderType;
   ReaderType::Pointer reader = ReaderType::New();
   reader->SetFileName (argv[1]);
   reader->Update();
   image = reader->GetOutput();
   }
 typedef itk::SobelEdgeDetectionImageFilter <ImageType, FloatImageType>
         SobelEdgeDetectionImageFilterType;
 SobelEdgeDetectionImageFilterType::Pointer sobelFilter
         = SobelEdgeDetectionImageFilterType::New();
 sobelFilter->SetInput(image);
 QuickView viewer;
 if (argc > 1)
   {
   viewer.AddImage<ImageType>(
     image,
     true,
     itksys::SystemTools::GetFilenameName(argv[1]));  
   }
 else
   {
   viewer.AddImage<ImageType>(image);
   }
 viewer.AddImage<FloatImageType>(
   sobelFilter->GetOutput(),
   true,
   "SobelEdgeDetectionImageFilter");
 viewer.Visualize();
 return EXIT_SUCCESS;

}

void CreateImage(ImageType::Pointer image) {

 itk::Index<2> start;
 start.Fill(0);
 itk::Size<2> size;
 size.Fill(100);
 itk::ImageRegion<2> region(start, size);
 image->SetRegions(region);
 image->Allocate();
 // Make a square
 for(unsigned int r = 20; r < 80; r++)
   {
   for(unsigned int c = 20; c < 80; c++)
     {
     itk::Index<2> pixelIndex;
     pixelIndex[0] = r;
     pixelIndex[1] = c;
     image->SetPixel(pixelIndex, 15);
     }
   }

} </source>

CMakeLists.txt

<syntaxhighlight lang="cmake"> cmake_minimum_required(VERSION 3.9.5)

project(SobelEdgeDetectionImageFilter)

find_package(ITK REQUIRED) include(${ITK_USE_FILE}) if (ITKVtkGlue_LOADED)

 find_package(VTK REQUIRED)
 include(${VTK_USE_FILE})

else()

 find_package(ItkVtkGlue REQUIRED)
 include(${ItkVtkGlue_USE_FILE})
 set(Glue ItkVtkGlue)

endif()

add_executable(SobelEdgeDetectionImageFilter MACOSX_BUNDLE SobelEdgeDetectionImageFilter.cxx) target_link_libraries(SobelEdgeDetectionImageFilter

 ${Glue}  ${VTK_LIBRARIES} ${ITK_LIBRARIES})

</syntaxhighlight>

Download and Build SobelEdgeDetectionImageFilter

Click here to download SobelEdgeDetectionImageFilter. and its CMakeLists.txt file. Once the tarball SobelEdgeDetectionImageFilter.tar has been downloaded and extracted,

cd SobelEdgeDetectionImageFilter/build 
  • If ITK is installed:
cmake ..
  • If ITK is not installed but compiled on your system, you will need to specify the path to your ITK build:
cmake -DITK_DIR:PATH=/home/me/itk_build ..

Build the project,

make

and run it:

./SobelEdgeDetectionImageFilter

WINDOWS USERS PLEASE NOTE: Be sure to add the VTK and ITK bin directories to your path. This will resolve the VTK and ITK dll's at run time.

Building All of the Examples

Many of the examples in the ITK Wiki Examples Collection require VTK. You can build all of the the examples by following these instructions. If you are a new VTK user, you may want to try the Superbuild which will build a proper ITK and VTK.

ItkVtkGlue

ITK >= 4

For examples that use QuickView (which depends on VTK), you must have built ITK with Module_ITKVtkGlue=ON.

ITK < 4

Some of the ITK Examples require VTK to display the images. If you download the entire ITK Wiki Examples Collection, the ItkVtkGlue directory will be included and configured. If you wish to just build a few examples, then you will need to download ItkVtkGlue and build it. When you run cmake it will ask you to specify the location of the ItkVtkGlue binary directory.