Optimal width and height after image rotation

A while ago, a blog reader Fabien sent me some code (you can read it here. Thanks Fabien!). The PHP code is a modification of my image rotation code with some upgrades.

I was looking through his code (French variable names!) and was puzzled by the initial section. I believe he based his code on my code where the resulting image wasn’t clipped after rotation, meaning the whole image was still in the picture/bitmap (though rotated).

In that piece of code, I just used the diagonal length of the image (from top-left corner to bottom-right corner) as the final length and breadth of the resulting image. This gave the simplest resulting image dimension without doing complicated maths calculations (a square in this case).

However, what if you want to know the optimal width and height of the resulting image after rotation? Meaning the best-fit width and height that just manages to contain the resulting rotated image. For that, I need to tell you some basic trigonometry and geometry.

Suppose you have a rectangle with L as the length and H as the height. It is rotated t angles. I’m not going to explain the maths behind it. It involves complementary angles, supplementary angles, rotation symmetry and trigonometry with sines and cosines. Convince yourself that the diagram is true.

So after rotating t angles, the optimal width is L * cos(t) + H * cos(90 – t)

The optimal height is L * sin(t) + H * sin(90 – t)

Short digression: You might notice that any lengths that lie parallel to the x-axis usually involve cosines, and lengths that lie parallel to the y-axis usually involve sines. It’s just the way trigonometry works.

Now, although the image rotation is carried out with respect to the image’s centre, rotating by the top-left corner will result in the same optimal width and height. Again, this is basic maths so you’ll just have to convince yourself it’s true (and that I don’t really want to explain it…).

But that’s if t is an acute angle. What about other angles?

For those angles, we just need to calculate the acute angle based on the initial rotation angle. After that, just substitute that calculated acute angle into our formula above. I have absolute confidence in your ability to check which quadrant in the Cartesian coordinate system does your rotation angle lie in.

UPDATE: In case you are unable to view images, if your rotation angle is in the 2nd quadrant, the calculated angle is (180 – t). If in the 3rd quadrant, it’s (t – 180). And if in the 4th quadrant, it’s (360 – t).

In practice, you might still want to pad a couple of pixels around. But that should give you the smallest image dimension which can still contain your rotated image.

Image rotation with bilinear interpolation and alpha progressive borders

So a blog reader, Fabien AurÃ©jac, emailed me an improvement over the code I posted on image rotation. Here’s the one with bilinear interpolation, and here’s the one with bilinear interpolation and no clipping.

You can find Fabien here and here (warning: French ahead), and also at Stack Overflow.

Fabien translated the core of my code into PHP. The improvement was on assigning alpha values to the edge pixels after rotation. Edge pixels are pixels beside the “blank” pixels (I used black in my code, for instance). The alpha values mean the edge pixels are “softer” and thus the resulting image looks smoother.

I suppose if you really want to, you could also “dumb down” the values of the red, green and blue colour components for more softening (in addition to the alpha component). I say “dumb down” because the blank pixels I used are black (meaning zero for the RGB values). You’re free to go ahead and do more interpolation.

Fabien has given permission for me to post the code here. I’ll leave it as an exercise for you to translate to your programming language.

```\$distOmbre=3;
\$flouOmbre=4;
\$angleRot=60;
\$img=imagecreatefromjpeg("media/diapo-Chinon.jpg");
\$size=getimagesize("media/diapo-Chinon.jpg");
\$LsupH=(\$size[0]>\$size[1])?1:0;
\$angleBool=(int)(\$angleRot/90)%2==0?0:1;
if ((\$angleBool+\$LsupH)%2==0) {
\$largeur=round(abs(\$size[0]*sin(\$angleRot%90*pi()/180))+abs(\$size[1]*sin((90-\$angleRot%90)*pi()/180)));
\$hauteur=round(abs(\$size[0]*cos(\$angleRot%90*pi()/180))+abs(\$size[1]*cos((90-\$angleRot%90)*pi()/180)));
} else {
\$largeur=round(abs(\$size[0]*cos(\$angleRot%90*pi()/180))+abs(\$size[1]*cos((90-\$angleRot%90)*pi()/180)));
\$hauteur=round(abs(\$size[0]*sin(\$angleRot%90*pi()/180))+abs(\$size[1]*sin((90-\$angleRot%90)*pi()/180)));
}
\$largeur+=\$distOmbre+\$flouOmbre*2;
\$hauteur+=\$distOmbre+\$flouOmbre*2;
\$angleRot*=pi()/180;
\$imgRot=imagecreatetruecolor(\$largeur, \$hauteur);
imagealphablending(\$imgRot, true);
imageantialias(\$imgRot, true);
for (\$i=0; \$i<\$hauteur; \$i++) {
for (\$j=0; \$j<\$largeur; \$j++) {
// convert raster to Cartesian
\$x = \$j - \$largeur*0.5;
\$y = \$hauteur*0.5 - \$i;

// convert Cartesian to polar
\$fDistance = sqrt(\$x * \$x + \$y * \$y);
\$fPolarAngle = atan2(\$y, \$x);

// the crucial rotation part
// "reverse" rotate, so minus instead of plus
\$fPolarAngle -= \$angleRot;
// convert polar to Cartesian
\$fTrueX = \$fDistance * cos(\$fPolarAngle);
\$fTrueY = \$fDistance * sin(\$fPolarAngle);

// convert Cartesian to raster
\$fTrueX = \$fTrueX + \$size[0]*0.5;
\$fTrueY = \$size[1]*0.5 - \$fTrueY;

\$iFloorX = (int)(floor(\$fTrueX));
\$iFloorY = (int)(floor(\$fTrueY));
\$iCeilingX = (int)(ceil(\$fTrueX));
\$iCeilingY = (int)(ceil(\$fTrueY));
//echo \$fTrueX." ".\$fTrueY." ".\$iFloorX." ".\$iCeilingX." ".\$iFloorY." ".\$iCeilingY."<br>";
if (\$iFloorX >= 0 && \$iCeilingX >= 0 && \$iFloorX < \$size[0] && \$iCeilingX < \$size[0] && \$iFloorY >= 0 && \$iCeilingY >= 0 && \$iFloorY < \$size[1] && \$iCeilingY < \$size[1]) {
\$fDeltaX = \$fTrueX - \$iFloorX;
\$fDeltaY = \$fTrueY - \$iFloorY;
\$clrTopLeft = imagecolorat(\$img, \$iFloorX, \$iFloorY);
\$colorsTopLeft = imagecolorsforindex(\$img, \$clrTopLeft);
\$clrTopRight = imagecolorat(\$img, \$iCeilingX, \$iFloorY);
\$colorsTopRight = imagecolorsforindex(\$img, \$clrTopRight);
\$clrBottomLeft = imagecolorat(\$img, \$iFloorX, \$iCeilingY);
\$colorsBottomLeft = imagecolorsforindex(\$img, \$clrBottomLeft);
\$clrBottomRight = imagecolorat(\$img, \$iCeilingX, \$iCeilingY);
\$colorsBottomRight = imagecolorsforindex(\$img, \$clrBottomRight);
// linearly interpolate horizontally between top neighbours
\$fTopRed = (1 - \$fDeltaX) * \$colorsTopLeft['red'] + \$fDeltaX * \$colorsTopRight['red'];
\$fTopGreen = (1 - \$fDeltaX) * \$colorsTopLeft['green'] + \$fDeltaX * \$colorsTopRight['green'];
\$fTopBlue = (1 - \$fDeltaX) * \$colorsTopLeft['blue'] + \$fDeltaX * \$colorsTopRight['blue'];
// linearly interpolate horizontally between bottom neighbours
\$fBottomRed = (1 - \$fDeltaX) * \$colorsBottomLeft['red'] + \$fDeltaX * \$colorsBottomRight['red'];
\$fBottomGreen = (1 - \$fDeltaX) * \$colorsBottomLeft['green'] + \$fDeltaX * \$colorsBottomRight['green'];
\$fBottomBlue = (1 - \$fDeltaX) * \$colorsBottomLeft['blue'] + \$fDeltaX * \$colorsBottomRight['blue'];
// linearly interpolate vertically between top and bottom interpolated results
\$iRed = (int)(round((1 - \$fDeltaY) * \$fTopRed + \$fDeltaY * \$fBottomRed));
\$iGreen = (int)(round((1 - \$fDeltaY) * \$fTopGreen + \$fDeltaY * \$fBottomGreen));
\$iBlue = (int)(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;
if (\$iFloorX > 0 && \$iCeilingX > 0 && \$iFloorX < \$size[0]-1 && \$iCeilingX < \$size[0]-1 && \$iFloorY > 0 && \$iCeilingY > 0 && \$iFloorY < \$size[1]-1 && \$iCeilingY < \$size[1]-1) {
\$colorallocation=imagecolorallocate(\$imgRot, \$iRed, \$iGreen, \$iBlue);
imagesetpixel(\$imgRot, \$j, \$i, \$colorallocation);
} else if (\$iFloorX == 0 && \$iFloorY >= 0 && \$iCeilingY >= 0 && \$iFloorY < \$size[1] && \$iCeilingY < \$size[1]) {//left
\$alpha=round((1-abs(\$fDeltaX))*127);
\$colorallocation=imagecolorallocatealpha(\$imgRot, \$iRed, \$iGreen, \$iBlue, \$alpha);
imagesetpixel(\$imgRot, \$j, \$i, \$colorallocation);
} else if (\$iFloorX >= 0 && \$iCeilingX >= 0 && \$iFloorX < \$size[0] && \$iCeilingX < \$size[0] && \$iFloorY == 0) {//top
\$alpha=round((1-abs(\$fDeltaY))*127);
\$colorallocation=imagecolorallocatealpha(\$imgRot, \$iRed, \$iGreen, \$iBlue, \$alpha);
imagesetpixel(\$imgRot, \$j, \$i, \$colorallocation);
} else if (\$iCeilingX == \$size[0]-1 && \$iFloorY >= 0 && \$iCeilingY >= 0 && \$iFloorY < \$size[1] && \$iCeilingY < \$size[1]) {//right
\$alpha=round(abs(\$fDeltaX)*127);
\$colorallocation=imagecolorallocatealpha(\$imgRot, \$iRed, \$iGreen, \$iBlue, \$alpha);
imagesetpixel(\$imgRot, \$j, \$i, \$colorallocation);
} else if (\$iFloorX >= 0 && \$iCeilingX >= 0 && \$iFloorX < \$size[0] && \$iCeilingX < \$size[0] && \$iCeilingY == \$size[1]-1) {//bottom
\$alpha=round(abs(\$fDeltaY)*127);
\$colorallocation=imagecolorallocatealpha(\$imgRot, \$iRed, \$iGreen, \$iBlue, \$alpha);
imagesetpixel(\$imgRot, \$j, \$i, \$colorallocation);
}
}
}
}
```

Fabien is French (I think), which is why you get variable names such as distOmbre (shadow distance?), flouOmbre (fuzzy shadow?), largeur (width), hauteur (height). And this one took me a bit more time to translate… LsupH is probably “width greater than height?”. The “L” probably refers to “largeur”, and “H” refers to “hauteur”.

Reading international programming code is fun. *smile*

There’s also an interesting piece of code:

```\$size=getimagesize("media/diapo-Chinon.jpg");
\$LsupH=(\$size[0]>\$size[1])?1:0;
\$angleBool=(int)(\$angleRot/90)%2==0?0:1;
if ((\$angleBool+\$LsupH)%2==0) {
\$largeur=round(abs(\$size[0]*sin(\$angleRot%90*pi()/180))+abs(\$size[1]*sin((90-\$angleRot%90)*pi()/180)));
\$hauteur=round(abs(\$size[0]*cos(\$angleRot%90*pi()/180))+abs(\$size[1]*cos((90-\$angleRot%90)*pi()/180)));
} else {
\$largeur=round(abs(\$size[0]*cos(\$angleRot%90*pi()/180))+abs(\$size[1]*cos((90-\$angleRot%90)*pi()/180)));
\$hauteur=round(abs(\$size[0]*sin(\$angleRot%90*pi()/180))+abs(\$size[1]*sin((90-\$angleRot%90)*pi()/180)));
}
```

So here’s your mission, should you choose to accept it (I recently watched Mission Impossible…). What is Fabien trying to accomplish in that section of code? Hint: it has something to do with getting a “nice” resulting image width and height.

I’ll tell you a more “elegant” alternative to that code section. But it’ll involve some mathematics. And drawings. Prepare for poorly drawn diagrams…

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.

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.

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.

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.

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:

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.

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:

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.

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).

Converting between raster, Cartesian and polar coordinates

As with the article on linear interpolation, this article is in preparation for the upcoming discussion on image rotation. I just realised that if I don’t write this, I’m going to take a long time explaining the image rotation code.

We’ve already covered Cartesian coordinates. So what are raster and polar coordinates?

Raster coordinates

I did some research, and to my surprise, there’s no such thing as a raster coordinate system! So where did this term enter my memory? Hmm…

That’s fine, we’ll define it here then. Raster coordinates are actually very simple.

The horizontal axis is similar to the usual x-axis, just only with non-negative values. The vertical axis is an inverted y-axis (so values increase downwards instead of upwards), and also has only non-negative values. I’ve only used raster coordinates with images, so negative indices don’t make sense.

For illustration, suppose w is the width of the image and h is the height of the image (in pixels). Then (0,0) is the top left corner of the image. (w,0) is the top right corner of the image. (0,h) is the bottom left corner of the image. And (w,h) is the bottom right corner of the image.

You’ll encounter raster coordinates when you deal with texture mapping. But that’s another story for another day…

Polar coordinates

Polar coordinates have two components, a length and an angle. 2D Cartesian coordinates also have two components, an x and a y. So what’s the difference?

Note that angles are measured from the positive x-axis and goes anti-clockwise. Remember the quadrants of the 2D Cartesian plane? I’ve included an example in the illustration with the line in the 3rd quadrant.

So how are polar coordinates related to Cartesian coordinates?

You’ll have to revise your trigonometry lessons. I’ll leave it to you to find out the derivation.

Converting from raster to Cartesian to polar coordinates. And back.

Why would anyone convert from raster to Cartesian to polar, only to convert back from polar to Cartesian to raster? Ahh… let’s look at a diagram.

We can only do proper rotation at the polar coordinate stage. But we start with an image, with raster coordinates. So we convert from (image) raster coordinates to Cartesian, then to polar, do the rotation part, convert back to Cartesian, and then back to raster coordinates.

To determine the formula for raster-Cartesian conversion, let’s look at the four corners of the image. We want
raster (0,0) -> Cartesian (-w/2,h/2)
raster (w,0) -> Cartesian (w/2,h/2)
raster (0,h) -> Cartesian (-w/2,-h/2)
raster (w,h) -> Cartesian (w/2,-h/2)

Based on that, the formula is
x = rasterX – w/2
y = h/2 – rasterY

For Cartesian-polar conversion, we have
r = sqrt(x*x + y*y)
theta = PI/2 if x=0 and y>0
theta = 3*PI/2 if x=0 and y<0 theta = arctan(y/x) otherwise I don't need to restrict theta to be within [0,2*PI) interval, even though it's mentioned here.

[short digression]
The square bracket [ means 0 is included in the interval. The round bracket ) means 2*PI is not included. The sine and cosine functions take in any real values. 2*PI + 1 radians automatically wraps to 1 radian by sine and cosine.
[end digression]

I don’t really need to explain why there are 3 conditional assignments for theta, right?

Alright, fine. x is a denominator as a parameter in the arctan function. That means you need to check for x equal to zero. Next, if x=0, the angle can either be 90 degrees (positive y-axis) or 270 degrees (negative y-axis). Hence spawning the other 2 conditions.

For polar-Cartesian conversion, we have
x = r * cos(theta)
y = r * sin(theta)

For Cartesian-raster conversion, we have
rasterX = x – w/2
rasterY = h/2 – x

As for some of the edge cases, we’ll look at them when we get to the code. Oh yes, there’s code. Lot’s of it. Stay tuned.

UPDATE: “I’m dying to look at the code for rotating images with bilinear interpolation! Bring me there!”

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.

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:

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

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.

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.