[Insight-users] Modifying the output (buffer) of an image filter

[Stefan Engstrom]

I have an image processing task that is fairly well described by
existing filters in ITK. It takes an input image and a kernel, finds a
locally smoothed version of the input image and tries to subtract that
from the original input before that is passed on to a detection
correlation with the kernel.

The local average is over the same size as the input kernel which can
be large so the calculation is done as a correlation in Fourier space.
I accomplish that with the fftw filters, by multiplying the complex
transforms using image iterators. The essential part of the code that
does that looks some thing like this:

  workRegStart[0]=0;
  workRegStart[1]=0;
  workReg.SetSize(
        fftlmInput->GetOutput()->GetLargestPossibleRegion().GetSize());
  workReg.SetIndex(workRegStart);

  ComplexConstIteratorType iterfftlmKernel(fftlmKernel->GetOutput(),workReg);
  ComplexIteratorType iterfftlmInput(fftlmInput->GetOutput(),workReg);

  for(iterfftlmKernel.GoToBegin(),iterfftlmInput.GoToBegin();
        !iterfftlmKernel.IsAtEnd(); ++iterfftlmKernel,++iterfftlmInput) {
    iterfftlmInput.Set(iterfftlmKernel.Get()*iterfftlmInput.Get());
  }

The point of showing this here is that the fflmInput buffer is
manipulated directly by the iterator. A little later I have the input
image and a smoothed (shifted) version of the local mean calculation
that I am trying to calculate the difference of and again put back
into the filter buffer, but this time it does not work. The relevant
code looks like this:

  inputResampler->Update(); // we just converted the input image to
something useful

// Write a snapshot of what we are working on
  tmprescaler->SetInput( inputResampler->GetOutput() );
  tmprescaler->SetOutputMinimum(  0  );
  tmprescaler->SetOutputMaximum( 255 );
  tmpwriter->SetFileName("tmp-inputimage.png");
  tmpwriter->SetInput(tmprescaler->GetOutput());
  tmpwriter->Update();

  WorkImageType::RegionType wr1,wr2;
  WorkImageType::RegionType::IndexType wrStart1,wrStart2;
  WorkImageType::RegionType::SizeType wrSize;


  wrStart1[0]=kernelsize[0];
  wrStart1[1]=kernelsize[1];
  wrStart2[0]=kernelsize[0]/2; // The evaluation point for the local
mean is shifted
  wrStart2[1]=kernelsize[0]/2; // by half the kernelsize
  wrSize[0]=inputsize[0];
  wrSize[1]=inputsize[1];
  wr1.SetSize(wrSize);
  wr2.SetSize(wrSize);
  wr1.SetIndex(wrStart1);
  wr2.SetIndex(wrStart2);

  WorkIteratorWithIndexType iter1(inputResampler->GetOutput(),wr1);
  WorkIteratorWithIndexType iter2(lmFlipper->GetOutput(),wr2);

  WorkPixelType
kernelarea=static_cast<WorkPixelType>(kernelsize[0]*kernelsize[1]);
  WorkPixelType a,lmean;
  for(iter1.GoToBegin(),iter2.GoToBegin();
        !iter1.IsAtEnd(); ++iter1,++iter2) {
    lmean=iter2.Get()/kernelarea;
    a=iter1.Get();
    iter1.Set(a-lmean);
      std::cout << "subtract: idx1=" << iter1.GetIndex() <<
        ", idx2=" << iter2.GetIndex() <<", a=" << a <<
        ", b=" << lmean << ", diff=" << iter1.Get() << std::endl;
  }

The indicies look good. I get values for 'a' and iter1.Get() is
properly set in the loop. When I try to use this image I get the same
thing I had coming in.

Is the Fourier routine special in some way in that it lets me modify
the buffer directly?
What is the recommended approach for this? I am already worried about
excessive memory use by having a longish string of filters connected
up. It seems like very poor use of memory to create so many instances
of the images for the filter flow.
An image adapter does not seem to be the answer since I need to
evaluate the result of two images (which have different origins).

Any leads, practial or philosophical to this would be greatly appreciated.

Stefan