[Insight-users] an Image To Image to compute the gradient of the popular l2 approximation of the total variation semi norm

[devieill at irit.fr]

Hi all,

I am currently rewriting some of my code, and as a matter of fact it find
it iteresting to write a filter to easily compute the gradient of the
total variation. I have already done this computation with a standard itk
pipeline, at a high cost (swapping) for large image.
So starting from itkGradientMagnitudeImageFilter I have decided to give it
a try. Right now I am bit stuck at computing the properly the second
derivative pass. To explain better from an image I with 2 dimensions, one
has to get the partial derivatives dI/dx, dI/dy and then to form 2 images
:
I1 = (dI/dx)/( (dI/dx)^2 + (dI/dy)^2 + beta^2)
I2 = (dI/dy)/( (dI/dx)^2 + (dI/dy)^2 + beta^2)
And finally to compute the result R:
R = dI1/dx + dI2/dy

this part is a bit tricky to me since I am not sure to have the proper
images to compute it, that is to say either I1 or I2 have size 1 on one of
the coordinates. Assuming I do a 2nd derivative, I should request a region
padded by a radius of 2 right ?

However I am not too sure about how to doing it properly and more
importantly if that is the critical point. Here is the important part of
the code doing this.

Any help would be much appreciated !!

Regards,
de Vieilleville François


template< class TInputImage, class TOutputImage>
void
GradientTotalVariationL2ApproximationImageFilter< TInputImage, TOutputImage >
::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
                       int threadId)
{
  //////////////////////////////////////////////
  //  create copies for storing the returned values
  typename TOutputImage::Pointer ddx[ImageDimension];

  //allocate copies
  for (unsigned int i = 0; i < ImageDimension; ++i) {
    ddx[i] = TOutputImage::New();
    ddx[i]->SetRegions(outputRegionForThread );
    ddx[i]->Allocate();
  }

  OutputPixelType grad_tv;
  ZeroFluxNeumannBoundaryCondition<InputImageType> nbc;

  ConstNeighborhoodIterator<InputImageType> nit;
  ConstNeighborhoodIterator<TInputImage> bit;
  ImageRegionIterator<OutputImageType> it;

  NeighborhoodInnerProduct<InputImageType, RealType> SIP;

  // Get the input and output
  OutputImageType *       outputImage = this->GetOutput();
  const InputImageType *  inputImage  = this->GetInput();

  // Set up operators
  DerivativeOperator<RealType, ImageDimension> op[ImageDimension];

  for (int i = 0; i< ImageDimension; i++)
    {
    op[i].SetDirection(0);
    op[i].SetOrder(1);
    op[i].CreateDirectional();

    // Reverse order of coefficients for the convolution with the image to
    // follow.
    //op[i].FlipAxes();

    // Take into account the pixel spacing if necessary
    if (m_UseImageSpacing == true)
      {
      if ( this->GetInput()->GetSpacing()[i] == 0.0 )
        {
        itkExceptionMacro(<< "Image spacing cannot be zero.");
        }
      else
        {
        op[i].ScaleCoefficients( 1.0 / this->GetInput()->GetSpacing()[i] );
        }
      }
    }

  // Calculate iterator radius
  Size<ImageDimension> radius;
  for (unsigned int i = 0; i < ImageDimension; ++i)
    {
      radius[i]  = op[0].GetRadius()[0];
      //std::cout << "radius["<<i<<"]=" << radius[i] << std::endl;
    }

  // Find the data-set boundary "faces"
  typename
NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<TInputImage>::FaceListType
faceList;
  NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<TInputImage> bC;
  faceList = bC(inputImage, outputRegionForThread, radius);

  typename
NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>::FaceListType::iterator
fit;
  fit = faceList.begin();

  // support progress methods/callbacks
  ProgressReporter progress(this, threadId,
outputRegionForThread.GetNumberOfPixels());

  // Initialize the x_slice array
  nit = ConstNeighborhoodIterator<InputImageType>(radius, inputImage, *fit);

  //////////////////////////////////////////////

  std::slice x_slice[ImageDimension];
  const unsigned long center = nit.Size() / 2;
  for (unsigned int i = 0; i < ImageDimension; ++i)
    {

      x_slice[i] = std::slice( center - nit.GetStride(i) * radius[0],
				 op[i].GetSize()[0], nit.GetStride(i));
    }
  // Process non-boundary face and then each of the boundary faces.
  // These are N-d regions which border the edge of the buffer.

  // set the proper iterators on the copied regions
  ImageRegionIterator<OutputImageType> it_copy[ImageDimension];
  //  typename TInputImage::RegionType region_on_copy_current_face;

  for (fit = faceList.begin(); fit != faceList.end(); ++fit)
    {
      bit = ConstNeighborhoodIterator<InputImageType>(radius,
						      inputImage, *fit);
      it = ImageRegionIterator<OutputImageType>(outputImage, *fit);
      // set the copied regions iterators
      for (unsigned int j =0; j < ImageDimension; ++j) {
	it_copy[j] = ImageRegionIterator<InputImageType>(ddx[j], *fit);
      }

      bit.OverrideBoundaryCondition(&nbc);
      bit.GoToBegin();

      while ( ! bit.IsAtEnd() )
	{
	  RealType a = NumericTraits<RealType>::Zero;
	  RealType deriv[ImageDimension];
	  for (int i = 0; i < ImageDimension; ++i)
	    {
	      const RealType g = SIP(x_slice[i], bit, op[i]);
	      deriv[i] = g;
	      a += g * g;
	      // vcl_pow(val,2.);
	    }
	  it.Value() = vcl_sqrt(a + m_BetaSquare);
	  for (unsigned int j =0; j < ImageDimension; ++j) {
	    it_copy[j].Value() = deriv[j]/it.Get();
	    ++(it_copy[j]);
	  }
	  ++bit;
	  ++it;
	}
    }

  faceList = bC(ddx[0], outputRegionForThread, radius); // assuming they
are all the same for the ddx[i]....
  fit = faceList.begin();
  ConstNeighborhoodIterator<InputImageType> nit_copy[ImageDimension];
  ConstNeighborhoodIterator<InputImageType> bit_copy[ImageDimension];

  for ( unsigned int i =0; i < ImageDimension; ++i) {
    nit_copy[i] = ConstNeighborhoodIterator<InputImageType>(radius,
							    ddx[i], *fit);
  }
  // Initialize the x_slice_copy array
  std::slice x_slice_copy[ImageDimension];
  for (unsigned int i = 0; i < ImageDimension; ++i) {
    const unsigned long center = nit_copy[i].Size() / 2;
    x_slice_copy[i] = std::slice( center - nit_copy[i].GetStride(i) *
radius[0],
				  op[i].GetSize()[0], nit_copy[i].GetStride(i));
  }

  for (fit = faceList.begin(); fit != faceList.end(); ++fit)
    {
      for ( unsigned int i =0; i < ImageDimension; ++i) {
	bit_copy[i] = ConstNeighborhoodIterator<InputImageType>(radius,
								ddx[i], *fit);
	bit_copy[i].OverrideBoundaryCondition(&nbc);
	bit_copy[i].GoToBegin();
      }
      it = ImageRegionIterator<OutputImageType>(outputImage, *fit);
      //bool nit_copy_are_at_end = false;
      while ( !bit_copy[0].IsAtEnd() )
	{
	  grad_tv =0.;
	  for (unsigned int i = 0; i < ImageDimension; ++i)
	    {
	      // compute derivative of previous computed region
	      double val = SIP(x_slice_copy[i], bit_copy[i], op[i]);
	      grad_tv += val;
	    }
	  it.Value() = grad_tv;
	  for (unsigned int j =0; j < ImageDimension; ++j) {
	    ++(bit_copy[j]);
	  }
	  ++it;
	  progress.CompletedPixel();
	}
    }
}