[ITK] About pipeline recursion

[Nicolas Gallego]

Good morning,

I have been coding a filter that implements an iterative procedure on
images. I have tried the straight forward implementation of updating a
temporal image on the recursion loop but results are wrong. The correct
result is obtained when I explicitly unroll the pipeline although it
becomes harder to code for more complex pipelines, and I am concerned about
the scalability on the number of recursions.

Could you give me some ideas to work this point out. Is there a way to
implement a straight forward way in itk? Here the code that illustrates the
point on recursively adding images:

----- code 1 -------------
        std::cout << "Recursion on image pointer" << std::endl;

    ImageType::Pointer image1 = ImageType::New();

    image1->SetRegions( region );

    image1->Allocate();

    image1->FillBuffer( 1 );

    AdderType::Pointer adder = AdderType::New();

    adder->SetInput1( image1 );

    adder->SetInput2( image1 );

    adder->Update();

    ImageType::Pointer tmpImage = adder->GetOutput();

    std::cout << "start: " << image1->GetPixel( index ) << std::endl;

    const unsigned short N = 10;

    for ( unsigned short k = 0; k < N; k++ ) {

        adder->Update();

        tmpImage = adder->GetOutput();

        std::cout << k << " " << tmpImage->GetPixel( index ) << std::endl;

        adder->SetInput1( tmpImage );

        adder->SetInput2( image1 );

    }

    std::cout << std::endl;

--- end of 1 --------------

Code 1 produces the following output:
-----------------------------------------------
Recursion on image pointer
start: 1
0 2
1 3
2 3
3 3
4 4
5 4
6 4
7 5
8 5
9 5
---------------------------------------------

The explicit pipeline unroll is illustrated in code 2

--- code 2 -----------------------------------------------
std::cout << "Pipeline explicit unroll recursion" << std::endl;

    std::vector< AdderType::Pointer > adders( N );

    adders[0] = AdderType::New();

    adders[0]->SetInput1( image1 );

    adders[0]->SetInput2( image1 );

    std::cout << "start: " << image1->GetPixel( index ) << std::endl;

    // connections

    for ( unsigned short k = 1; k < N; k++ ) {

        adders[k] = AdderType::New();

        adders[k]->SetInput1( adders[k-1]->GetOutput() );

        adders[k]->SetInput2( image1 );

    }

    // execution

    adders[N-1]->Update();

    for ( unsigned short k = 0; k < N; k++ ) {

        std::cout << k  << " " << adders[k]->GetOutput()->GetPixel(
index ) << std::endl;

    }


--- end of code 2 --------------------------------------

Pipeline explicit unroll recursion
start: 1
0 2
1 3
2 4
3 5
4 6
5 7
6 8
7 9
8 10
9 11
------------------------------------

Find attached complete code file, thank you for your help

------
Nicolás Gallego-Ortiz
Université catholique de Louvain, Belgium
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20140620/ba8223ac/attachment-0001.html>
-------------- next part --------------
#include <itkImage.h>
#include <itkAddImageFilter.h>

int main()
{
    typedef unsigned short PixelType;
    const unsigned int Dimension = 3;
    typedef itk::Image< PixelType, Dimension >  ImageType;
    typedef itk::AddImageFilter< ImageType, ImageType, ImageType > AdderType;

    ImageType::RegionType region;
    ImageType::IndexType index;
    index.Fill(0);
    ImageType::SizeType size;
    size.Fill(1);
    region.SetIndex( index );
    region.SetSize( size );

    std::cout << "Recursion on image pointer" << std::endl;

    ImageType::Pointer image1 = ImageType::New();
    image1->SetRegions( region );
    image1->Allocate();
    image1->FillBuffer( 1 );

    AdderType::Pointer adder = AdderType::New();
    adder->SetInput1( image1 );
    adder->SetInput2( image1 );
    adder->Update();
    ImageType::Pointer tmpImage = adder->GetOutput();

    std::cout << "start: " << image1->GetPixel( index ) << std::endl;
    const unsigned short N = 10;
    for ( unsigned short k = 0; k < N; k++ ) {

        adder->Update();
        tmpImage = adder->GetOutput();
        std::cout << k << " " << tmpImage->GetPixel( index ) << std::endl;

        adder->SetInput1( tmpImage );
        adder->SetInput2( image1 );
    }
    std::cout << std::endl;
    std::cout << "Pipeline explicit unroll recursion" << std::endl;

    std::vector< AdderType::Pointer > adders( N );
    adders[0] = AdderType::New();
    adders[0]->SetInput1( image1 );
    adders[0]->SetInput2( image1 );

    std::cout << "start: " << image1->GetPixel( index ) << std::endl;

    // connections
    for ( unsigned short k = 1; k < N; k++ ) {

        adders[k] = AdderType::New();
        adders[k]->SetInput1( adders[k-1]->GetOutput() );
        adders[k]->SetInput2( image1 );
    }

    // execution
    adders[N-1]->Update();
    for ( unsigned short k = 0; k < N; k++ ) {
        std::cout << k  << " " << adders[k]->GetOutput()->GetPixel( index ) << std::endl;
    }

    return 0;
}