BiCubic Interpolation does not match any known filter!

[NicolasRobidoux]

Anthony: Suggestion: I believe that many people expect Catmull-Rom when they ask for a "bicubic". Why don't you simply let me optimize the code for this bicubic so it actually produces Catmull-Rom and not "explicitly" add Catmull-Rom? Nobody in his/her right mind should use the current one over Catmull-Rom. The only situation in which it possibly could be better is for images that are so smooth that they don't exist. So I think it should simply "disappear".
I've done the optimization for VIPS, so should be fast.
-----
As you know, I don't care much for Hermite, but I could optimize it's code too. You'd have to plug the pipes and I'd need to ask you a few simple questions, but this should not take long.
-----
What you say?
P.S. I don't really care because I get what I want in other ways. I'm not being pushy here. Just trying to be helpful. So, only accept the offer if it's a plus value to you.

[anthony]

Anthony: Suggestion: I believe that many people expect Catmull-Rom when they ask for a "bicubic". Why don't you simply let me optimize the code for this bicubic so it actually produces Catmull-Rom and not "explicitly" add Catmull-Rom? Nobody in his/her right mind should use the current one over Catmull-Rom. The only situation in which it possibly could be better is for images that are so smooth that they don't exist. So I think it should simply "disappear".
I've done the optimization for VIPS, so should be fast.
-----
As you know, I don't care much for Hermite, but I could optimize it's code too. You'd have to plug the pipes and I'd need to ask you a few simple questions, but this should not take long.
-----
What you say?
P.S. I don't really care because I get what I want in other ways. I'm not being pushy here. Just trying to be helpful. So, only accept the offer if it's a plus value to you.

Okay. Though from what I saw in IMv6 (where I have been doing this work) The code is even worse! with sub-functions.

Maybe you should look at the 'Spline' code too.

If you can make it so I can simply plugin the constants for Spline, Catrom, and Cardinal. that will be good.
Hermite is different as it only needs 2 image values for the interpolation, not 4.

[NicolasRobidoux]

Anthony: I'll post the "inner loop" bits here starting tomorrow. I won't touch the svn directly.

[NicolasRobidoux]

Anthony: Just to make sure the message got through: I really think that the current bicubic ("Cardinal") should simply disappear. My statement that it may be better that Catmull-Rom when the image is so smooth so as to be unreal is flat out wrong. It won't be better.
The backs of all my envelopes are unanimous in stating that this is a dreadful scheme.
IM is better without it.

[anthony]

You are probably right.

Still I would discontinue the use of 'bicubic' in example documentation as being 'in-exact' definition, even if it is still provided and points to catrom.

ASIDE: In IMv6 I try to not change 'constants' or the API as this cause library incompatibilities. IMv7 has no such restrictions while it remains in Alpha Development. Though I try not to upset your 'padawan's work

[NicolasRobidoux]

Suggestions:
1) Add Mitchell
2) Instead of Spline being cubic B-Spline smoothing, make it quadratic B-Spline smoothing. Tighter, and smoothing enough.

[NicolasRobidoux]

Anthony: I have not triple checked that it gives the right result (I've not run it, just compiled), but here is my rewrite of BicubicInterpolatePixel from pixel.c so that it computes Catmull-Rom efficiently.

Code: Select all

    case BicubicInterpolatePixel:
    {
      p=GetCacheViewVirtualPixels(image_view,x_offset-1,y_offset-1,4,4,
        exception);
      if (p == (const Quantum *) NULL)
        {
          status=MagickFalse;
          break;
        }
      if ((traits & BlendPixelTrait) == 0)
        for (i=0; i < 16; i++)
        {
          alpha[i]=1.0;
          pixels[i]=(MagickRealType) p[i*GetPixelChannels(image)+channel];
        }
      else
        for (i=0; i < 16; i++)
        {
          alpha[i]=QuantumScale*GetPixelAlpha(image,p+i*
            GetPixelChannels(image));
          pixels[i]=alpha[i]*p[i*GetPixelChannels(image)+channel];
        }
      {
        const double deltax = x - x_offset;
        const double deltay = y - y_offset;

        /*
         * Refactoring of the Catmull-Rom computation by Nicolas
         * Robidoux with 55 flops = 28* + 10- + 17+.
         * Originally implemented for the VIPS (Virtual Image
         * Processing System) library.
         */

        const double xtmp1 = 1.0 - deltax;
        const double xtmp2 = -0.5 * deltax;
        const double xtmp3 = xtmp1 * xtmp2;
        const double cxone = xtmp1 * xtmp3;
        const double cxfou = deltax * xtmp3;
        const double xtmp4 = cxfou - cxone;
        const double cxtwo = xtmp1 - cxone + xtmp4;
        const double cxthr = deltax - cxfou - xtmp4;

        const double ytmp1 = 1.0 - deltay;
        const double ytmp2 = -0.5 * deltay;
        const double ytmp3 = ytmp1 * ytmp2;
        const double cyone = ytmp1 * ytmp3;
        const double cyfou = deltay * ytmp3;
        const double ytmp4 = cyfou - cyone;
        const double cytwo = ytmp1 - cyone + ytmp4;
        const double cythr = deltay - cyfou - ytmp4;

        *pixel = cyone * ( cxone * pixels[0] +
                           cxtwo * pixels[1] +
                           cxthr * pixels[2] +
                           cxfou * pixels[3] )
                 +
                 cytwo * ( cxone * pixels[4] +
                           cxtwo * pixels[5] +
                           cxthr * pixels[6] +
                           cxfou * pixels[7] )
                 +
                 cythr * ( cxone * pixels[8] +
                           cxtwo * pixels[9] +
                           cxthr * pixels[10] +
                           cxfou * pixels[11] )
                 +
                 cyfou * ( cxone * pixels[12] +
                           cxtwo * pixels[13] +
                           cxthr * pixels[14] +
                           cxfou * pixels[15] );
      }
      break;
    }

It should run noticeably faster than the original. Also, I've only modified the first instance of it in the file: It appears in InterpolatePixelChannel (the one above) and in InterpolatePixelChannels (which I can similarly mod if you want).

[magick]

Nicolas, don't forget you need to interpolate alpha separately and blend it into the color channel results. Here's the refactored code. You / Anthony can update this method in the Subversion trunk. Don't forget there are 4 locations. Two places in ImageMagick 7.0.0-0 and two in ImageMagick 6.7.7-6.

Code: Select all

    case BicubicInterpolatePixel:
    {
      PointInfo
        delta;

      MagickRealType
        cx[4],
        cy[4],
        temp[4];

      /*
        Refactoring of the Catmull-Rom computation by Nicolas Robidoux with 55
        flops = 28* + 10- + 17+.  Originally implemented for the VIPS
        (Virtual Image Processing System) library.
      */
      p=GetCacheViewVirtualPixels(image_view,x_offset-1,y_offset-1,4,4,
        exception);
      if (p == (const Quantum *) NULL)
        {
          status=MagickFalse;
          break;
        }
      if ((traits & BlendPixelTrait) == 0)
        for (i=0; i < 16; i++)
        {
          alpha[i]=1.0;
          pixels[i]=(MagickRealType) p[i*GetPixelChannels(image)+channel];
        }
      else
        for (i=0; i < 16; i++)
        {
          alpha[i]=QuantumScale*GetPixelAlpha(image,p+i*
            GetPixelChannels(image));
          pixels[i]=alpha[i]*p[i*GetPixelChannels(image)+channel];
        }
      delta.x=x-x_offset;
      delta.y=y-y_offset;
      temp[0]=1.0-delta.x;
      temp[1]=(-0.5)*delta.x;
      temp[2]=temp[0]*temp[1];
      cx[0]=temp[0]*temp[2];
      cx[3]=delta.x*temp[2];
      temp[3]=cx[3]-cx[0];
      cx[1]=temp[0]-cx[0]+temp[3];
      cx[2]=delta.x-cx[3]-temp[3];
      temp[0]=1.0-delta.y;
      temp[1]=(-0.5)*delta.y;
      temp[2]=temp[0]*temp[1];
      cy[0]=temp[0]*temp[2];
      cy[3]=delta.y*temp[2];
      temp[3]=cy[3]-cy[0];
      cy[1]=temp[0]-cy[0]+temp[3];
      cy[2]=delta.y-cy[3]-temp[3];
      gamma=cy[0]*(cx[0]*alpha[0]+cx[1]*alpha[1]+cx[2]*alpha[2]+
        cx[3]*alpha[3])+cy[1]*(cx[0]*alpha[4]+cx[1]*alpha[5]+
        cx[2]*alpha[6]+cx[3]*alpha[7])+cy[2]*(cx[0]*alpha[8]+
        cx[1]*alpha[9]+cx[2]*alpha[10]+cx[3]*alpha[11])+cy[3]*(
        cx[0]*alpha[12]+cx[1]*alpha[13]+cx[2]*alpha[14]+cx[3]*alpha[15]);
      gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
      *pixel=gamma*(cy[0]*(cx[0]*pixels[0]+cx[1]*pixels[1]+cx[2]*pixels[2]+
        cx[3]*pixels[3])+cy[1]*(cx[0]*pixels[4]+cx[1]*pixels[5]+
        cx[2]*pixels[6]+cx[3]*pixels[7])+cy[2]*(cx[0]*pixels[8]+
        cx[1]*pixels[9]+cx[2]*pixels[10]+cx[3]*pixels[11])+cy[3]*(
        cx[0]*pixels[12]+cx[1]*pixels[13]+cx[2]*pixels[14]+cx[3]*pixels[15]));
      break;
    }

[fmw42]

Anthony: I have not triple checked that it gives the right result (I've not run it, just compiled), but here is my rewrite of BicubicInterpolatePixel from pixel.c so that it computes Catmull-Rom efficiently.

Code: Select all

    case BicubicInterpolatePixel:
    {
      p=GetCacheViewVirtualPixels(image_view,x_offset-1,y_offset-1,4,4,
        exception);
      if (p == (const Quantum *) NULL)
        {
          status=MagickFalse;
          break;
        }
      if ((traits & BlendPixelTrait) == 0)
        for (i=0; i < 16; i++)
        {
          alpha[i]=1.0;
          pixels[i]=(MagickRealType) p[i*GetPixelChannels(image)+channel];
        }
      else
        for (i=0; i < 16; i++)
        {
          alpha[i]=QuantumScale*GetPixelAlpha(image,p+i*
            GetPixelChannels(image));
          pixels[i]=alpha[i]*p[i*GetPixelChannels(image)+channel];
        }
      {
        const double deltax = x - x_offset;
        const double deltay = y - y_offset;

        /*
         * Refactoring of the Catmull-Rom computation by Nicolas
         * Robidoux with 55 flops = 28* + 10- + 17+.
         * Originally implemented for the VIPS (Virtual Image
         * Processing System) library.
         */

        const double xtmp1 = 1.0 - deltax;
        const double xtmp2 = -0.5 * deltax;
        const double xtmp3 = xtmp1 * xtmp2;
        const double cxone = xtmp1 * xtmp3;
        const double cxfou = deltax * xtmp3;
        const double xtmp4 = cxfou - cxone;
        const double cxtwo = xtmp1 - cxone + xtmp4;
        const double cxthr = deltax - cxfou - xtmp4;

        const double ytmp1 = 1.0 - deltay;
        const double ytmp2 = -0.5 * deltay;
        const double ytmp3 = ytmp1 * ytmp2;
        const double cyone = ytmp1 * ytmp3;
        const double cyfou = deltay * ytmp3;
        const double ytmp4 = cyfou - cyone;
        const double cytwo = ytmp1 - cyone + ytmp4;
        const double cythr = deltay - cyfou - ytmp4;

        *pixel = cyone * ( cxone * pixels[0] +
                           cxtwo * pixels[1] +
                           cxthr * pixels[2] +
                           cxfou * pixels[3] )
                 +
                 cytwo * ( cxone * pixels[4] +
                           cxtwo * pixels[5] +
                           cxthr * pixels[6] +
                           cxfou * pixels[7] )
                 +
                 cythr * ( cxone * pixels[8] +
                           cxtwo * pixels[9] +
                           cxthr * pixels[10] +
                           cxfou * pixels[11] )
                 +
                 cyfou * ( cxone * pixels[12] +
                           cxtwo * pixels[13] +
                           cxthr * pixels[14] +
                           cxfou * pixels[15] );
      }
      break;
    }

It should run noticeably faster than the original. Also, I've only modified the first instance of it in the file: It appears in InterpolatePixelChannel (the one above) and in InterpolatePixelChannels (which I can similarly mod if you want).

Nicolas,

Sorry to stick my nose in here. I don't want to make a lot of work for you here, but just wondering if this is limited to catmull-rom or can be done for the whole cubic family with B,C control. Would each B,C combination need separate optimization code?

Fred

[NicolasRobidoux]

just wondering if this is limited to catmull-rom or can be done for the whole cubic family with B,C control. Would each B,C combination need separate optimization code?

Generality adds flops. I did this work a long time ago (and it's not published: TODO list (very far from the top): intended title of article: Image Resampling Hacker's Delight) and I'm not sure I can track my notes, but basically:

BC-spline > Keys > individual ones (">" means "more flops")

and some of the individual ones optimize more than others.

If I don't find my notes, I'll have to go at it with Axiom (the CAS) and re-figure things out. Catmull-Rom was quick because I had programmed this into the VIPS library.

[NicolasRobidoux]

Nicolas, don't forget you need to interpolate alpha separately and blend it into the color channel results.

I was just about to post here RE: that it's weird that bilinear appears to gamma correct but that the current bicubic does not.
I'll look at things more carefully. (Multiple channels going on right now: my grad student with the EXQUIRES test suite, a test image bank contributor who wants the histogram to be tweaked more, this trying to improve the interpolate code... and my wife who called me about what to cook for dinner!).

[NicolasRobidoux]

I was just about to post here RE: that it's weird that bilinear appears to gamma correct but that the current bicubic does not.

I was confused by the use of the variable gamma: In this code, gamma has nothing to do with gamma correction, it's an interpolated transparency.

[NicolasRobidoux]

(deleted) I mistakenly edited this post instead of posting a reply. (Moral: Don't post before coffee.)

[NicolasRobidoux]

The reason why continuity matters is that otherwise round off error can have a large impact on pixel values. The current version is "unstable" (in the numerical analysis sense).

(magick: I did not write the above for you.)

[NicolasRobidoux]

Code self-documentation suggestion: Use "beta" instead of "gamma" so that the confusion is removed.