[Insight-users] an Image To Image to compute the gradient of the popular l2 approximation of the total variation semi norm
Luis Ibanez
luis.ibanez at kitware.com
Sun Oct 17 10:55:13 EDT 2010
Hi Francois,
> that is to say either I1 or I2 have size 1 on one of
> the coordinates.
I'm not sure that should be the case when you are
computing the gradient for a 2D image.
> Assuming I do a 2nd derivative, I should request a region
> padded by a radius of 2 right ?
Yes, assuming that you are using central differences
to compute the second derivative, you will need at
least 5 values (which is a radius = 2).
Regards,
Luis
-------------------------------------------------------------------------
On Thu, Sep 23, 2010 at 12:13 PM, <devieill at irit.fr> wrote:
> 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();
> }
> }
> }
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.html
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20101017/224c3e3a/attachment.htm>
More information about the Insight-users
mailing list