[Insight-users] Example on how to parallelize a simple neighborhood filter

[lynx.abraxas at freenet.de]

Hello again!


Since  my  last  post below I got a bit further. Since I skip the neighborhood
evaluation where voxels of the mask are I don't think the  parallelization  by
subregion (as mentioned in the users guide) is a good choise for my problem.
So  I've put my voxel operation in the iteration into a separate function. I'm
now wondering:

- Should I use pthreads or just a simple fork/child process to distribute each
voxel analysis on a separate core (of a 8 core machine)?

-  Would  I have to make a full copy of each iterator since they get increased
during the execution of other threads?

- Is there a much better way to go than mine?

Thanks for any help or hints on this.
Lynx
-------------- next part --------------
//calculate the mean of a masked binary image within a neighborhood
//optimized to only look at the neighborhood if the centre pixel is not from the mask
//mean_masked parallelized


#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkNeighborhoodAlgorithm.h"
#include <math.h>

#include "itkConstantBoundaryCondition.h"
#include "itkConstShapedNeighborhoodIterator.h"
#include "itkImageRegionIterator.h"

const unsigned int Dimension = 3;

typedef unsigned char  IPixelType;
//typedef float          OPixelType; //code needs to be modified for float!!!
typedef unsigned char  OPixelType; //easier to use

const IPixelType foreground_value, background_value, maskground_value; //*ground_value does not change during iteration -> global vars.

int operate_on_voxel(it, out, np){//*ground_value does not change during iteration -> global vars.
    unsigned int sum_f, sum_b, tot;
    IPixelType value;

    ShapedNeighborhoodIteratorType::ConstIterator ci;

    if (it.GetCenterPixel() == maskground_value){
        out.Set(0);
        return(0);
        }

    sum_f= 0;
    sum_b= 0;
    for (ci = it.Begin(); ci != it.End(); ci++){
        value= ci.Get();
        if (value == foreground_value)
            sum_f++;
        else if (value == background_value) //else if is quicker than a pure if!!!
            sum_b++;
        }
            
    tot= sum_f + sum_b;
    if (tot)
        //out.Set(static_cast<OPixelType>( static_cast<float>(sum_f) / tot * std::numeric_limits<OPixelType>::max())); //include scaling
        out.Set(static_cast<OPixelType>( sum_f * std::numeric_limits<OPixelType>::max() / tot)); //quicker?
    else out.Set(0);

    return(1);
    }


int main( int argc, char ** argv )
    {
    if ( argc < 7 ){
        std::cerr << "Missing parameters. " << std::endl;
        std::cerr << "Usage: " << std::endl;
        std::cerr << argv[0]
                  << " inputImageFile outputImageFile element_radius foreground background maskground"
                  << std::endl;
        return -1;
        }


    foreground_value = atoi(argv[4]);
    background_value = atoi(argv[5]);
    maskground_value = atoi(argv[6]); //maskground is just needed for const boundary cond.

    typedef itk::Image<IPixelType, Dimension>  ImageType;
    typedef itk::Image<OPixelType, Dimension>  OImageType;

    itk::ConstantBoundaryCondition<ImageType> BC;
    //BoundaryConditionType BoundaryCondition=  BoundaryConditionType::New();
    BC.SetConstant(maskground_value);

    //typedef itk::ConstShapedNeighborhoodIterator<ImageType, BoundaryCondition> ShapedNeighborhoodIteratorType;
    typedef itk::ConstShapedNeighborhoodIterator<ImageType> ShapedNeighborhoodIteratorType;

    typedef itk::ImageRegionIterator<ImageType> IteratorType;
    typedef itk::ImageRegionIterator<OImageType> OIteratorType;
    OIteratorType out;
  
    typedef itk::ImageFileReader< ImageType > ReaderType;
    ReaderType::Pointer reader = ReaderType::New();
    reader->SetFileName( argv[1] );
  
    try {
        reader->Update();
        }
    catch ( itk::ExceptionObject &err){
        std::cout << "ExceptionObject caught !" << std::endl; 
        std::cout << err << std::endl; 
        return -1;
        }

    OImageType::Pointer output = OImageType::New();
    output->SetRegions(reader->GetOutput()->GetRequestedRegion());
    output->Allocate();

    unsigned int element_radius = ::atoi( argv[3] );
    ShapedNeighborhoodIteratorType::RadiusType radius;
    radius.Fill(element_radius);

    typedef itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<ImageType> FaceCalculatorType;
  
    FaceCalculatorType faceCalculator;
    FaceCalculatorType::FaceListType faceList;
    FaceCalculatorType::FaceListType::iterator fit;
  
    faceList = faceCalculator( reader->GetOutput(), output->GetRequestedRegion(), radius );
  
    const float rad = static_cast<float>(element_radius);

    unsigned int nop, count, ll; //for progress info
    count= 0;
    nop = reader->GetOutput()->GetRequestedRegion().GetNumberOfPixels();
    ll = nop / 10000;

    std::cout << "Total number of voxels: " << nop << std::endl << std::flush;


    //split up the input region into border regions and a (big) none border region

    for ( fit=faceList.begin(); fit != faceList.end(); ++fit){
        
        ShapedNeighborhoodIteratorType it( radius, reader->GetOutput(), *fit );//iterate over each subregion
        it.OverrideBoundaryCondition(&BC); // assign the boundary condition
        out = OIteratorType( output, *fit );

        // Creates a circular structuring element by activating all the pixels less
        // than radius distance from the center of the neighborhood.
        ShapedNeighborhoodIteratorType::OffsetType off;
  
        for (float z = -rad; z <= rad; z++){
            for (float y = -rad; y <= rad; y++){
                for (float x = -rad; x <= rad; x++){
                    float dis = vcl_sqrt( x*x + y*y + z*z);
                    if (dis <= rad){
                        off[0] = static_cast<int>(x);
                        off[1] = static_cast<int>(y);
                        off[2] = static_cast<int>(z);
                        it.ActivateOffset(off); //if segfaults comment this! Why???
                        //because 3 dim (z, off[2]) was missing!!!!!!!
                        }
                    }
                }
            }


        ////now do the actual voxel iteration

        it.GoToBegin();
        out.GoToBegin();
        np= 0; //reset number of processes

        while(!it.IsAtEnd()){
            count++;
            if (np < max_p){
                np++;
                operate_on_voxel(it, out, np);//does it need a full copy of it and out since the pointers change during execution of thread?!?
                ++it;
                ++out;
                }
            else
                wait(); //waiting sensible?
            //print out some progess info
            if ( count % ll == 0 ){
                fprintf(stderr, "\rFilter progress: %6.2f%%", 100.0 * count / nop);
                std::cerr.flush();
                }
            }

        }
    }
        
    typedef itk::ImageFileWriter<OImageType> WriterType;
  
    WriterType::Pointer writer = WriterType::New();
    writer->SetFileName( argv[2] );
    writer->SetInput( output );
    try {
        writer->Update();
        }
    catch ( itk::ExceptionObject &err) {
        std::cout << "ExceptionObject caught !" << std::endl;
        std::cout << err << std::endl;
        return -1;
        }

    return 0;
    }