Image rotation with bilinear interpolation and no clipping

I wrote an article on image rotation with bilinear interpolation to calculate the pixel information. Commenter Jahn wanted to know if the resulting image could be obtained without clipping. I replied with the required code changes (it was fairly minimal), then I thought it might be useful to you too.

The code was cleaned up a bit, and here it is:

const int cnSizeBuffer = 20;
// 30 deg = PI/6 rad
// rotating clockwise, so it's negative relative to Cartesian quadrants
const double cnAngle = -0.52359877559829887307710723054658;
// use whatever image you fancy
Bitmap bm = new Bitmap("slantedtreetopsky.jpg");
// general iterators
int i, j;
// calculated indices in Cartesian coordinates
int x, y;
double fDistance, fPolarAngle;
// for use in neighbouring indices in Cartesian coordinates
int iFloorX, iCeilingX, iFloorY, iCeilingY;
// calculated indices in Cartesian coordinates with trailing decimals
double fTrueX, fTrueY;
// for interpolation
double fDeltaX, fDeltaY;
// pixel colours
Color clrTopLeft, clrTopRight, clrBottomLeft, clrBottomRight;
// interpolated "top" pixels
double fTopRed, fTopGreen, fTopBlue;
// interpolated "bottom" pixels
double fBottomRed, fBottomGreen, fBottomBlue;
// final interpolated colour components
int iRed, iGreen, iBlue;
int iCentreX, iCentreY;
int iDestCentre;
int iWidth, iHeight;
int iDiagonal;
iWidth = bm.Width;
iHeight = bm.Height;
iDiagonal = (int)(Math.Ceiling(Math.Sqrt((double)(bm.Width * bm.Width + bm.Height * bm.Height)))) + cnSizeBuffer;

iCentreX = iWidth / 2;
iCentreY = iHeight / 2;
iDestCentre = iDiagonal / 2;

Bitmap bmBilinearInterpolation = new Bitmap(iDiagonal, iDiagonal);

for (i = 0; i < iDiagonal; ++i)
{
	for (j = 0; j < iDiagonal; ++j)
	{
		bmBilinearInterpolation.SetPixel(j, i, Color.Black);
	}
}

// assigning pixels of destination image from source image
// with bilinear interpolation
for (i = 0; i < iDiagonal; ++i)
{
	for (j = 0; j < iDiagonal; ++j)
	{
		// convert raster to Cartesian
		x = j - iDestCentre;
		y = iDestCentre - i;

		// convert Cartesian to polar
		fDistance = Math.Sqrt(x * x + y * y);
		fPolarAngle = 0.0;
		if (x == 0)
		{
			if (y == 0)
			{
				// centre of image, no rotation needed
				bmBilinearInterpolation.SetPixel(j, i, bm.GetPixel(iCentreX, iCentreY));
				continue;
			}
			else if (y < 0)
			{
				fPolarAngle = 1.5 * Math.PI;
			}
			else
			{
				fPolarAngle = 0.5 * Math.PI;
			}
		}
		else
		{
			fPolarAngle = Math.Atan2((double)y, (double)x);
		}

		// the crucial rotation part
		// "reverse" rotate, so minus instead of plus
		fPolarAngle -= cnAngle;

		// convert polar to Cartesian
		fTrueX = fDistance * Math.Cos(fPolarAngle);
		fTrueY = fDistance * Math.Sin(fPolarAngle);

		// convert Cartesian to raster
		fTrueX = fTrueX + (double)iCentreX;
		fTrueY = (double)iCentreY - fTrueY;

		iFloorX = (int)(Math.Floor(fTrueX));
		iFloorY = (int)(Math.Floor(fTrueY));
		iCeilingX = (int)(Math.Ceiling(fTrueX));
		iCeilingY = (int)(Math.Ceiling(fTrueY));

		// check bounds
		if (iFloorX < 0 || iCeilingX < 0 || iFloorX >= iWidth || iCeilingX >= iWidth || iFloorY < 0 || iCeilingY < 0 || iFloorY >= iHeight || iCeilingY >= iHeight) continue;

		fDeltaX = fTrueX - (double)iFloorX;
		fDeltaY = fTrueY - (double)iFloorY;

		clrTopLeft = bm.GetPixel(iFloorX, iFloorY);
		clrTopRight = bm.GetPixel(iCeilingX, iFloorY);
		clrBottomLeft = bm.GetPixel(iFloorX, iCeilingY);
		clrBottomRight = bm.GetPixel(iCeilingX, iCeilingY);

		// linearly interpolate horizontally between top neighbours
		fTopRed = (1 - fDeltaX) * clrTopLeft.R + fDeltaX * clrTopRight.R;
		fTopGreen = (1 - fDeltaX) * clrTopLeft.G + fDeltaX * clrTopRight.G;
		fTopBlue = (1 - fDeltaX) * clrTopLeft.B + fDeltaX * clrTopRight.B;

		// linearly interpolate horizontally between bottom neighbours
		fBottomRed = (1 - fDeltaX) * clrBottomLeft.R + fDeltaX * clrBottomRight.R;
		fBottomGreen = (1 - fDeltaX) * clrBottomLeft.G + fDeltaX * clrBottomRight.G;
		fBottomBlue = (1 - fDeltaX) * clrBottomLeft.B + fDeltaX * clrBottomRight.B;

		// linearly interpolate vertically between top and bottom interpolated results
		iRed = (int)(Math.Round((1 - fDeltaY) * fTopRed + fDeltaY * fBottomRed));
		iGreen = (int)(Math.Round((1 - fDeltaY) * fTopGreen + fDeltaY * fBottomGreen));
		iBlue = (int)(Math.Round((1 - fDeltaY) * fTopBlue + fDeltaY * fBottomBlue));

		// make sure colour values are valid
		if (iRed < 0) iRed = 0;
		if (iRed > 255) iRed = 255;
		if (iGreen < 0) iGreen = 0;
		if (iGreen > 255) iGreen = 255;
		if (iBlue < 0) iBlue = 0;
		if (iBlue > 255) iBlue = 255;

		bmBilinearInterpolation.SetPixel(j, i, Color.FromArgb(iRed, iGreen, iBlue));
	}
}
bmBilinearInterpolation.Save("rotatedslantedtreetopsky.jpg", System.Drawing.Imaging.ImageFormat.Jpeg);

In the original article, the resulting image had the same dimensions as the source image. This meant that after rotation, some pixel information would be lost. To counter that, we make the dimensions of the resulting image “large enough”.

Considering that this is a rotation operation with an arbitrary angle, when all possible rotations of the source image are placed on top of one another, you get a smudgy circle. Not sure if you can visualise that. So the resulting image must be able to contain a circle, and because our images are stored in rectangular format, we need a square. That’s why in the code, I didn’t use a separate X and Y component of the dimensions, because they’re the same.

The width of the square has to be, at the minimum, the diagonal of the source image (which you can use Pythagoras’ Theorem to calculate). I added a buffer (of 20 pixels) to ensure that our square can contain the resulting image.

Here’s the source image:
Slanted treetop and sky

And here’s the resulting image:
Rotated slanted treetop with sky

Alright, have fun with the code.

Bilinear interpolation article referenced in another language

I’m thrilled. I’m also confused. I wrote an article on bilinear interpolation, specifically for image rotations, and it was provided as a link in a comment to this post. Thanks to vottini for referencing my article (it was used in a good light, right?)

I’ve got a problem. I don’t know what’s the language of the site! I think it’s French or Spanish, but I can’t be sure. If you know what’s the language, please let me know. Better yet, translate and summarise that article and tell me, because I only understood it as a example of how to expand an image and fill in the extra pixels. Thanks.

UPDATE: Thanks to an anonymous commentor, I finally found that the article (and site) is in Portuguese. I tried online translation tools and Portuguese-English yielded enough for me to understand what was going on. The commentor recommended the tool by Google, but I got barred with “automated request” error.

So I used the Babel Fish.

Excerpt of translated comment:

the solution, after all, was well simple. Instead of rotacionar pixels of the source, it rotacionava pixels of the destination and discovered in which place of the source would go to be. Joining with the technique of the bilinear interpolation or some most advanced one, the result is really impressive.

In the article, I suggested starting from the destination image, and find out what pixels to use from the source. This sort of assumes that both source and destination speak the same “language” (both are talking about RGB pixels).

The irony? In this incident, I don’t know anything about the source language (Portuguese), so I don’t know how to start from the destination (translated article in English).

UPDATE: Actually Christopher also commented that it’s in Portuguese (missed his comment in the diligent black hole SPAM processor…). Thanks Christopher!

And the owner of the site came to confirm that it’s in Portuguese! Wow. He (“o velho” means “old man” in Portuguese) apparently posted an English post about this, specially for me. Wow. Thanks.

Image rotation with bilinear interpolation

In this article, I’ll show you how to rotate an image about its centre. 3 assignment methods will be shown,

  • assign source pixels to destination pixels
  • assign destination pixels from source pixels
  • assign destination pixels from source pixels with bilinear interpolation

I’ll show you the code, then the results for comparison. So what’s bilinear interpolation?

Bilinear interpolation

Read up on linear interpolation first if you haven’t done so. “Bilinear” means there are 2 directions to interpolate. Let me illustrate.

Bilinear interpolation

In our case, we’re interpolating between 4 pixels. Visualise each pixel as a single point. Linearly interpolate between the top 2 pixels. Linearly interpolate between the bottom 2 pixels. Then linearly interpolate between the calculated results of the previous two.

You can expand on this concept to get trilinear interpolation.

Trilinear interpolation

LERPs is a short form of linear interpolations. When would trilinear interpolation be useful? Voxels, which is out of scope in this article.

Defining the centre of an image

I’m going to be fuzzy about this. I’m going to just take one pixel in the image and define it as the centre. This pixel is defined as having a horizontal index equal to half of its width (rounded down), and a vertical index equal to half its height (rounded down).

This means the image isn’t rotated about its “true” centre, but with a relatively large size, it won’t matter anyway. It’s not like you’re rotating an image of 5 pixel width and 3 pixel height, right?

The preparation part

The actual code is quite long, so I’m separating it into 4 parts.

  • Initialisation and variable declaration
  • Assigning source pixels to destination pixels
  • Assigning destination pixels from source pixels
  • Assigning destination pixels from source pixels with bilinear interpolation

It’s hard-coded with -30 degrees as the angle of rotation, but you can easily write it into a function.

// 30 deg = PI/6 rad
// rotating clockwise, so it's negative relative to Cartesian quadrants
const double cnAngle = -0.52359877559829887307710723054658;
// use whatever image you fancy
Bitmap bm = new Bitmap("rotationsource.jpg");
// general iterators
int i, j;
// calculated indices in Cartesian coordinates
int x, y;
double fDistance, fPolarAngle;
// for use in neighbouring indices in Cartesian coordinates
int iFloorX, iCeilingX, iFloorY, iCeilingY;
// calculated indices in Cartesian coordinates with trailing decimals
double fTrueX, fTrueY;
// for interpolation
double fDeltaX, fDeltaY;
// pixel colours
Color clrTopLeft, clrTopRight, clrBottomLeft, clrBottomRight;
// interpolated "top" pixels
double fTopRed, fTopGreen, fTopBlue;
// interpolated "bottom" pixels
double fBottomRed, fBottomGreen, fBottomBlue;
// final interpolated colour components
int iRed, iGreen, iBlue;
int iCentreX, iCentreY;
int iWidth, iHeight;
iWidth = bm.Width;
iHeight = bm.Height;

iCentreX = iWidth / 2;
iCentreY = iHeight / 2;

Bitmap bmSourceToDestination = new Bitmap(iWidth, iHeight);
Bitmap bmDestinationFromSource = new Bitmap(iWidth, iHeight);
Bitmap bmBilinearInterpolation = new Bitmap(iWidth, iHeight);

for (i = 0; i < iHeight; ++i)
{
    for (j = 0; j < iWidth; ++j)
    {
        // initialise when "throwing" values
        bmSourceToDestination.SetPixel(j, i, Color.Black);
        // since we're looping, we might as well do for the others
        bmDestinationFromSource.SetPixel(j, i, Color.Black);
        bmBilinearInterpolation.SetPixel(j, i, Color.Black);
    }
}

Some of it might not mean anything to you yet. Just wait for the rest of the code. You might want to read up on converting between raster, Cartesian and polar coordinates first before moving on.

Throwing values from source to destination

// assigning pixels from source image to destination image
for (i = 0; i < iHeight; ++i)
{
    for (j = 0; j < iWidth; ++j)
    {
        // convert raster to Cartesian
        x = j - iCentreX;
        y = iCentreY - i;

        // convert Cartesian to polar
        fDistance = Math.Sqrt(x * x + y * y);
        fPolarAngle = 0.0;
        if (x == 0)
        {
            if (y == 0)
            {
                // centre of image, no rotation needed
                bmSourceToDestination.SetPixel(j, i, bm.GetPixel(j, i));
                continue;
            }
            else if (y < 0)
            {
                fPolarAngle = 1.5 * Math.PI;
            }
            else
            {
                fPolarAngle = 0.5 * Math.PI;
            }
        }
        else
        {
            fPolarAngle = Math.Atan2((double)y, (double)x);
        }

        // the crucial rotation part
        fPolarAngle += cnAngle;

        // convert polar to Cartesian
        x = (int)(Math.Round(fDistance * Math.Cos(fPolarAngle)));
        y = (int)(Math.Round(fDistance * Math.Sin(fPolarAngle)));

        // convert Cartesian to raster
        x = x + iCentreX;
        y = iCentreY - y;

        // check bounds
        if (x < 0 || x >= iWidth || y < 0 || y >= iHeight) continue;

        bmSourceToDestination.SetPixel(x, y, bm.GetPixel(j, i));
    }
}
bmSourceToDestination.Save("rotationsrctodest.jpg", System.Drawing.Imaging.ImageFormat.Jpeg);

It should be fairly easy to read. Note the part about checking for the central pixel of the image. No rotation calculation necessary, so we assign and move to the next pixel. Note also the part about checking boundaries.

Finding values from the source

// assigning pixels of destination image from source image
for (i = 0; i < iHeight; ++i)
{
    for (j = 0; j < iWidth; ++j)
    {
        // convert raster to Cartesian
        x = j - iCentreX;
        y = iCentreY - i;

        // convert Cartesian to polar
        fDistance = Math.Sqrt(x * x + y * y);
        fPolarAngle = 0.0;
        if (x == 0)
        {
            if (y == 0)
            {
                // centre of image, no rotation needed
                bmDestinationFromSource.SetPixel(j, i, bm.GetPixel(j, i));
                continue;
            }
            else if (y < 0)
            {
                fPolarAngle = 1.5 * Math.PI;
            }
            else
            {
                fPolarAngle = 0.5 * Math.PI;
            }
        }
        else
        {
            fPolarAngle = Math.Atan2((double)y, (double)x);
        }

        // the crucial rotation part
        // "reverse" rotate, so minus instead of plus
        fPolarAngle -= cnAngle;

        // convert polar to Cartesian
        x = (int)(Math.Round(fDistance * Math.Cos(fPolarAngle)));
        y = (int)(Math.Round(fDistance * Math.Sin(fPolarAngle)));

        // convert Cartesian to raster
        x = x + iCentreX;
        y = iCentreY - y;

        // check bounds
        if (x < 0 || x >= iWidth || y < 0 || y >= iHeight) continue;

        bmDestinationFromSource.SetPixel(j, i, bm.GetPixel(x, y));
    }
}
bmDestinationFromSource.Save("rotationdestfromsrc.jpg", System.Drawing.Imaging.ImageFormat.Jpeg);

The key difference here is the use of the rotation angle. Instead of adding it, we subtract it. The reason is, we rotate source pixels 30 degrees clockwise and assign it to destination pixels. But from destination pixels, we get source pixels which are rotated 30 degrees anticlockwise. Either way, we get a destination image that's the source image rotated 30 degrees clockwise.

Rotating the other way

Also compare the assignment, noting the indices:

bmSourceToDestination.SetPixel(x, y, bm.GetPixel(j, i));
bmDestinationFromSource.SetPixel(j, i, bm.GetPixel(x, y));

x and y variables are calculated and thus "messy". I prefer my messy indices on the right. There's a practical reason for it too, which will be evident when I show you the rotation results.

Image rotation code with bilinear interpolation

// assigning pixels of destination image from source image
// with bilinear interpolation
for (i = 0; i < iHeight; ++i)
{
    for (j = 0; j < iWidth; ++j)
    {
        // convert raster to Cartesian
        x = j - iCentreX;
        y = iCentreY - i;

        // convert Cartesian to polar
        fDistance = Math.Sqrt(x * x + y * y);
        fPolarAngle = 0.0;
        if (x == 0)
        {
            if (y == 0)
            {
                // centre of image, no rotation needed
                bmBilinearInterpolation.SetPixel(j, i, bm.GetPixel(j, i));
                continue;
            }
            else if (y < 0)
            {
                fPolarAngle = 1.5 * Math.PI;
            }
            else
            {
                fPolarAngle = 0.5 * Math.PI;
            }
        }
        else
        {
            fPolarAngle = Math.Atan2((double)y, (double)x);
        }

        // the crucial rotation part
        // "reverse" rotate, so minus instead of plus
        fPolarAngle -= cnAngle;

        // convert polar to Cartesian
        fTrueX = fDistance * Math.Cos(fPolarAngle);
        fTrueY = fDistance * Math.Sin(fPolarAngle);

        // convert Cartesian to raster
        fTrueX = fTrueX + (double)iCentreX;
        fTrueY = (double)iCentreY - fTrueY;

        iFloorX = (int)(Math.Floor(fTrueX));
        iFloorY = (int)(Math.Floor(fTrueY));
        iCeilingX = (int)(Math.Ceiling(fTrueX));
        iCeilingY = (int)(Math.Ceiling(fTrueY));

        // check bounds
        if (iFloorX < 0 || iCeilingX < 0 || iFloorX >= iWidth || iCeilingX >= iWidth || iFloorY < 0 || iCeilingY < 0 || iFloorY >= iHeight || iCeilingY >= iHeight) continue;

        fDeltaX = fTrueX - (double)iFloorX;
        fDeltaY = fTrueY - (double)iFloorY;

        clrTopLeft = bm.GetPixel(iFloorX, iFloorY);
        clrTopRight = bm.GetPixel(iCeilingX, iFloorY);
        clrBottomLeft = bm.GetPixel(iFloorX, iCeilingY);
        clrBottomRight = bm.GetPixel(iCeilingX, iCeilingY);

        // linearly interpolate horizontally between top neighbours
        fTopRed = (1 - fDeltaX) * clrTopLeft.R + fDeltaX * clrTopRight.R;
        fTopGreen = (1 - fDeltaX) * clrTopLeft.G + fDeltaX * clrTopRight.G;
        fTopBlue = (1 - fDeltaX) * clrTopLeft.B + fDeltaX * clrTopRight.B;

        // linearly interpolate horizontally between bottom neighbours
        fBottomRed = (1 - fDeltaX) * clrBottomLeft.R + fDeltaX * clrBottomRight.R;
        fBottomGreen = (1 - fDeltaX) * clrBottomLeft.G + fDeltaX * clrBottomRight.G;
        fBottomBlue = (1 - fDeltaX) * clrBottomLeft.B + fDeltaX * clrBottomRight.B;

        // linearly interpolate vertically between top and bottom interpolated results
        iRed = (int)(Math.Round((1 - fDeltaY) * fTopRed + fDeltaY * fBottomRed));
        iGreen = (int)(Math.Round((1 - fDeltaY) * fTopGreen + fDeltaY * fBottomGreen));
        iBlue = (int)(Math.Round((1 - fDeltaY) * fTopBlue + fDeltaY * fBottomBlue));

        // make sure colour values are valid
        if (iRed < 0) iRed = 0;
        if (iRed > 255) iRed = 255;
        if (iGreen < 0) iGreen = 0;
        if (iGreen > 255) iGreen = 255;
        if (iBlue < 0) iBlue = 0;
        if (iBlue > 255) iBlue = 255;

        bmBilinearInterpolation.SetPixel(j, i, Color.FromArgb(iRed, iGreen, iBlue));
    }
}
bmBilinearInterpolation.Save("rotationbilinearinterpolation.jpg", System.Drawing.Imaging.ImageFormat.Jpeg);

This part is similar to the destination-from-source part, with a lot more calculations. We have to find the 4 pixels that surrounds our "true" position-calculated pixel. Then we perform linear interpolation on the 4 neighbouring pixels.

We need to interpolate for the red, green and blue components individually. Refer to my article on colour theory for a refresher.

Pictures, pictures!

After doing all that, we're finally done. Let me show you my source image first.

Rotation source

I added the marble cylinder for emphasising image quality. I needed something that's straight (vertically or horizontally) in the source image.

Here's what we get after rotating with the source-to-destination method:

Rotation - source to destination

Note the speckled black pixels dotted all over. This is because some of the destination pixels (which are within the image bounds) were unassigned.

Note also that the specks even form patterns. This is due to the sine and cosine functions, and the regularity of pixel width and height. Sine and cosine are periodic functions. Since pixel indices are regular, therefore sine and cosine results are regular too. Hence, calculations regularly fail to assign pixel values.

There might be source pixels that have the same calculated destination pixel (due to sine and cosine and rounding). This also implies that there might be anti-gravity destination pixels that no source pixel can ever matched to! I haven't verified this, but it seems a possibility.

Anti-gravity destination pixel

Still think you should iterate over the source (image/array) instead of over the destination?

Next, we have the image result of the destination-from-source method:

Rotation - destination from source

Compare the quality with the source-to-destination part. No missing pixels. It's still sort of grainy though. This is because some of the destination pixels get their values from the same source pixel, so there might be 2 side-by-side destination pixels with the same colour. This gives mini blocks of identical colour in the result, which on the whole, gives an unpolished look.

Now, we have the bilinear interpolation incorporated version.

Rotation - destination from source with bilinear interpolation

It looks smoother, right? Note the straight edge of the marble cylinder. Compare with the image result without bilinear interpolation.

I might do a version with cubic interpolation for even smoother results, but I feel the bilinear version is good enough for now. Have fun!

[UPDATE: Now, there's a version without the resulting image being clipped.]

If you enjoyed this article and found it useful, please share it with your friends. You should also subscribe to get the latest articles (it's free).

Linear and cubic interpolation

Interpolation is a method of calculating a value from a set of given values. We’ll be looking at interpolation with a bias towards image processing, but the theory can be generalised for other purposes. You’ve probably already solved some interpolation problems without knowing it. Let me give you an example.

A distance problem

Suppose there are 3 towns A, B, C and they happen to lie on a straight line, in that order. B is 5 kilometres away from A, and C is 15 kilometres away from A. If you travel one quarter of the way from town B to town C, how far are you from town A?

Distance problem

To solve it, you can figure out the distance between B and C, which is 15 – 5 = 10 km. One quarter of the way means 1/4 * 10 = 2.5 km. Then add the distance between A and B to this and you have 5 + 2.5 = 7.5 km.

Linear interpolation

If you visualise the problem as interpolating between 2 points, then B becomes the point p0 with a value of 5 (km) and C becomes the point p1 with a value of 15 (km). The usual variable used is t, so the generic formula is:
f(t) = (1 – t) * p0 + t * p1, where t lies between 0 and 1 inclusive.

Using this, we have
f(1/4) = (1 – 1/4) * 5 + 1/4 * 15
= 3/4 * 5 + 1/4 * 15
= 7.5

This is linear interpolation. Linearity refers to the power of the variable t, which is 1. Note that there’s no stopping you from using negative values of t or values greater than 1.

Suppose you travelled from B to A one quarter of the distance between B and C. How far are you from town A?
f(-1/4) = (1 – (-1/4)) * 5 + (-1/4) * 15
= 5/4 * 5 – 1/4 * 15
= 2.5

Suppose you travelled from B to C and went past C by a quarter of the way. How far are you from town A?
f(5/4) = (1 – 5/4) * 5 + 5/4 * 15
= -1/4 * 5 + 5/4 * 15
= 17.5

What happens if you get a negative result?
f(-1) = (1 – (-1)) * 5 + (-1) * 15
= 2 * 5 – 15
= -5

It means you’re 5 kilometres away from town A. You’re just in the opposite direction from towns B and C. The calculation result is correct. It’s how you interpret the value.

Applications in image processing

A common operation in image processing is manipulating height maps. Height maps are usually greyscale bitmap files where a white pixel (RGB values are 255 for all 3) is the highest point, and a black pixel (RGB values are 0 for all 3) is the lowest point.

Terrain editor in Bryce

You know enlarging photographs can give you some weird results. What happens is you’re trying to fill in the blanks in a larger image using values from the original image. Where do you think the image editing software comes up with values? Interpolation.

If you think of the red, green and blue values of image pixels as 3 different “height maps”, then you’re just performing interpolation on 3 values. Suppose we’re talking about linear interpolation between two pixels. You’ll interpolate between the red component of the 2 pixels and get a value. Similarly you do it for the green and blue components. The calculated results of the red, green and blue become the interpolated colour.

Cubic Bezier interpolation

There are all kinds of cubic curves available. The Catmull–Rom spline, the non-uniform rational B-spline (NURBS) and I didn’t really want to write anything on the subject after I remember my Hermite splines… I love Bezier curves though, so I thought maybe I can write something with that.

Instead of 2 points used in linear interpolation, cubic interpolation uses 4 points. To illustrate, suppose you’re on an undulating plain with small hills undulating in their usual carefree manner. You’re in between two such (undulating) hills and you want to find out how high you are.

Instead of linear interpolating your way through these two (undulating) hills, the better way will be to interpolate with another 2 (undulating) hills! Ok, I’m stopping with the undulating thing…

Undulating hill interpolation

The Bezier curve equation looks like this:
B(t) = (1-t)^3 * p0 + 3*(1-t)^2 * t * p1 + 3*(1-t)* t^2 * p2 + t^3 * p3
where p0, p1, p2, p3 are the (height) values, and t lies between 0 and 1 inclusive.

You will be between p1 and p2. Let’s also assume that the hills are equidistant from each other. Like the pixels on an image, the hills shall be of equal distance from its neighbour.

Because of this equidistant property, p1 is 0.33 (roughly 1/3) units away from p0, p2 is 0.67 (roughly 2/3) units away from p0 and p3 is 1 unit away from p0.

How do you know what’s the value of t to use? You might be able to calculate the t if you do linear interpolation between p1 and p2. But that t value is different from the t value in the Bezier curve.

Ahhh… once you get the t-linear value, you interpolate with 0.33 and 0.67 to get the t-Bezier value. Confused? Suppose you’re one quarter way from p1 to p2. Your t-linear value is 1/4. Interpolate that with 0.33 and 0.67 to get
f(1/4) = (1 – 1/4) * 0.33 + 1/4 * 0.67
= 0.415

And 0.415 is your t-Bezier value. Voila!

You skipped the quadratic power!

I know. It’s logical to think that there’s a power 2 somewhere. But there isn’t. There is one fundamental flaw with quadratic interpolation. Which segment do you use?

Quadratic interpolation has issues...

In closing

Interpolation is just a method of creating data values using a set of existing data. What those created values mean is up to you to interpret.

In image processing, interpolation can be used to fill in blanks when enlarging an image. It doesn’t guarantee that the enlarged image looks good. Image processing is very much an aesthetic-based operation. I’ll talk a bit more on this when I get to writing code to rotate images.

Messy indices on the right please

In case you don’t know, indices are the plural form of index. Wait, isn’t the plural form of index, indexes? Well, indexes is the plural form, if you’re talking about database table indexes (I’m confused at first too). I’m referring to the offset of an array, the index of an array.

So what do I mean by messy indices on the right? Variable assignment. I mentioned something of this in the matrix rotation article, and I want to talk more on it here. And I will tell you why I prefer them on the right side of the assignment. First, I need to mention raytracing.

Tracing a light ray from the eye back to the scene

In 3D graphics, a scene is rendered and you need to know what colour each resulting pixel is. You can throw light rays everywhere, let them hit an object, bounce off that object, and continue bouncing and if that ray happen to shoot through to the camera, render it. That doesn’t sound very efficient, does it?

There are millions of light rays to shoot from a light source. They can bounce off objects in any number of ways. And only a small percentage of them will ever enter the camera.

Raytracing simplified

What you can do is, start from the camera. Trace a ray from the camera through the resulting pixel, back out into the scene and see what it hits. If it doesn’t hit an object, it will travel all the way into the sky or space. Render accordingly. If it does hit an object, continue to track the bounced (or reflected) ray until it hits nothing. Or hit a light source. Either way, you know how to render the colour. It’s a bit more involved than that, and I’m simplifying the explanation.

Basically you start from the end.

Rotating the square matrix

Let’s revisit that rotating-the-matrix problem.

For Y = 0 to 3
    For X = 0 to 3
        Destination(3-Y,X) = Source(X,Y)
    Next
Next

Essentially, the author concentrated on the source, iterating through it with appropriate indices and assigning those values to the destination. See how “messy” the indices on left side are? I call it, the “throw values over” assignment:

Throwing values over

Then there’s the raytracing way, the “finding values” way:

Finding values assignment

There doesn’t seem to be any difference between the two. Until I tell you about rotating images…

Rotating an image 30 degrees clockwise

Suppose you want to rotate an image in non-right-angled degrees, say 30 degrees clockwise.

Rotated image outside rectangle

We can discuss the code to rotate the image some other time. Hmm… new article fodder…

Anyway, you can start from the source image, iterating through each pixel and “throw values” to the destination image. Or you can start iterating the destination image and “find values” from the source image. Each method requires you to check the boundaries of the image (index must be within width and height).

If both seem equivalent, what’s the problem? Image quality.

You see, with the algorithm involved, there’s the possibility that you’d throw different pixels from the source image onto the same pixel on the destination image. You’d be doing sine’s and cosine’s and rounding the calculation into an integer so you can map it to an index. This isn’t good, since you can override a previous calculation.

With the “find value” way, you have a similar problem. There’s the possibility that for different pixels on the destination image, you’d map back to the same pixel on the source image. However, you don’t have to use the exact colour of that source pixel. You can use interpolation to get a sort of “average” colour from the surrounding pixels.

You know, I should write something on the image rotation part. Alright, coming up soon, an image processing article focused on image rotation!

And that’s the reason why…

I prefer the messy indices on the right side of the assignment. Because I start from the end (or destination), the iteration indices are “clean” on the left side.

And now, homework! See if you can come up with code (or algorithm) to rotate an image any number of degrees, assuming the rotation is about the image’s centre. You can make it easier by rotating about the top-left corner. And you have to iterate through the pixels, not use some in-built function to manipulate the image. You can also make the image black and white, so you focus on grey values instead of an RGB triplet, to make the code simpler to write.

If you’re adventurous, and want to try the interpolation thingy to make the resulting image look better, you can research on linear interpolation and cubic splines.