How to rotate a 2D array at non integral angle

Question

Application:

I am doing analysis of a laser beam which is imaged on an IR camera. The goal is the have a real-time measurement, or at least as fast as the camera, which refreshes about once every 125 ms.

The algorithm will sum the intensities between the yellow dashed lines, and outside of them, then find the ratio. The yellow dashed lines are rotated at the same angle as the beam's angle. The beam's angle is almost never a multiple of 45 degrees.

I receive the beam in the form of a 2D array from the camera. It is represented by a 32-bit Integer array of size 1360 x 1024. I want to rotate the array so the yellow dashed lines are perpendicular to the bounds of the array, which would make the summing algorithm simpler.

What I've tried:

I looked into System.Drawing.Drawing2D.Matrix.Rotate(), as has been suggested in answers to similar questions. However, this applies to Graphics objects. Using that method I have rotated images. But I haven't found a way to perform a similar operation on a 2D array.

Edit

The algorithm actually finds w such that the lines encompass n% of the beam's energy. I have tried traversing the entire array without rotation and checking the boundaries using the lines (y = mx + b), however, since I do many iterations of this operation, it is too costly. Again, I have 1.4 million values in the entire image.

The reason I want to rotate is to calculate sums of each row of the rotated array (resulting in a beam profile), then find the narrowest w which produces n% enclosed energy.

The following image shows the rotated beam, the profile (1D array of sums, in red), and the width (blue).

Answer 1

Here is one way to use the drawing transforms to rotate an int[,]. The demo code below needs a new forms application with

• timer1 (Enabled, ~100ms)
• pictureBox1 (1360x1024, hook up paint handler)

This outputs the rotated data to rotArr and completes the operation in ~40ms (on my 4 year-old laptop). Compiling for x86 seems to be faster here. It may be just quick enough if you can manage the rest efficiently and have a reasonably new system to work with. To get better speeds than this you can override the default interpolation mode using :

 g.InterpolationMode = InterpolationMode.NearestNeighbor;

this reduces the quality of the rotated image but it can be faster. A low interpolation quality may add uncertainty to your data - you'll have to test this. Another flag you can test is

g.SmoothingMode = SmoothingMode.AntiAlias;

This adds time to the calculation (not much for me, but your mileage may vary), but may improve the fidelity of the rotated image (reduce artifact errors introduced by the transform).

This has to be C# - you said you'd prefer VB but you have to set "Allow unsafe code" in your project compile options for this work and VB does not support unsafe code.

using System;
using System.Drawing;
using System.Windows.Forms;
using System.Drawing.Drawing2D;
using System.Runtime.InteropServices;
using System.Drawing.Imaging;

namespace WindowsFormsApplication1
{
    public partial class Form1 : Form
    {
        private int[,] imgArr = new int[1360, 1024];  //the 2D source from the camera
        private int[,] rotArr = new int[1360, 1024];  //the rotated output
        private float angle = 0;

        public Form1()
        {
            InitializeComponent();
            //Make an RGB stripe for testing
            for (int i = 0; i < 1360; i++)
            {
                for (int j = 0; j < 1024; j++)
                {
                    if (j < 333)
                    {
                        imgArr[i, j] = -65536; //0xFFFF0000 - red
                    }
                    else if (j < 666)
                    {
                        imgArr[i, j] = -16711936; //0xFF00FF00 - green
                    }
                    else
                    {
                        imgArr[i, j] = -16776961; //0xFF0000FF - blue
                    }                       
                }
            }          

        }

        private void pictureBox1_Paint(object sender, PaintEventArgs e)
        {
            Bitmap drawImg;
            // Copy array to a bitmap - fast!
            unsafe
            {
                fixed (int* intPtr = &imgArr[0, 0])
                {
                    //swap width and height due to [,] row layout
                    drawImg = new Bitmap(1024, 1360, 4096,
                                  PixelFormat.Format32bppArgb, 
                                  new IntPtr(intPtr));                        
                }
            }

            // transform 
            GraphicsPath gp = new GraphicsPath();
            gp.AddPolygon(new Point[]{new Point(0,0), 
                          new Point(drawImg.Width,0), 
                          new Point(0,drawImg.Height)});
            Matrix m = new Matrix();
            m.RotateAt(angle, new PointF((float)drawImg.Width/2, (float)drawImg.Height/2)); 
            gp.Transform(m);
            PointF[] pts = gp.PathPoints;

            //draw rotated image
            Bitmap rotImg = new Bitmap(1024, 1360);
            Graphics g = Graphics.FromImage(rotImg);

            // for speed...default is Bilinear
            //g.InterpolationMode = InterpolationMode.NearestNeighbor;

            // may improve accuracy - default is no AntiAliasing
            // slows calculation
            // g.SmoothingMode = SmoothingMode.AntiAlias;

            g.DrawImage(drawImg, pts);

            //copy array data out 
            BitmapData bData = rotImg.LockBits(
                                        new Rectangle(new Point(), rotImg.Size),
                                        ImageLockMode.ReadOnly,
                                        PixelFormat.Format32bppArgb);
            int byteCount = bData.Stride * (rotImg.Height);
            int[] flatArr = new int[byteCount / 4];
            //Do this in two steps - first to a flat array, then 
            //block copy to the [,] array
            Marshal.Copy(bData.Scan0, flatArr, 0, byteCount / 4);
            Buffer.BlockCopy(flatArr, 0, rotArr, 0, byteCount);

            // unlock the bitmap!!
            rotImg.UnlockBits(bData); 

            // for show... draw the rotated image to the picturebox
            // have to rotate 90deg due to [,] row layout
            rotImg.RotateFlip(RotateFlipType.Rotate90FlipNone);
            e.Graphics.DrawImage(rotImg, new Point(0, 0));

        }

        private void timer1_Tick(object sender, EventArgs e)
        {
            // increment angle and paint
            angle += 1.0f;
            if (angle > 360.0f) { angle = 0; }
            pictureBox1.Invalidate();
        }

    }
 }

Answer 2

I think I have the conversion from pixel to physical coordinates correct (and back). Then what you want is this:

public struct ImageData
{
    /// <summary>
    /// Intensity map
    /// </summary>
    int[,] intensities;
    /// <summary>
    /// Pixel dimensios of image like 1360 × 1024
    /// </summary>
    Size pixel_size;

    /// <summary>
    /// Physical dimensions like 300μ × 260μ
    /// </summary>
    SizeF actual_size;

    /// <summary>
    /// Transforms pixel coordinates to actual dimensions. Assumes center of image is (0,0)
    /// </summary>
    /// <param name="pixel">The pixel coordinates (integer i,j)</param>
    /// <rereturns>The physical coordinates (float x,y)</rereturns>
    public PointF PixelToPoint(Point pixel)
    {
        return new PointF(
            actual_size.Width*((float)pixel.X/(pixel_size.Width-1)-0.5f),
            actual_size.Height*((float)pixel.Y/(pixel_size.Height-1)-0.5f));
    }
    /// <summary>
    /// Transforms actual dimensions to pixels. Assumes center of image is (0,0)
    /// </summary>
    /// <param name="point">The point coordines (float x,y)</param>
    /// <returns>The pixel coordinates (integer i,j)</returns>
    public Point PointToPixel(PointF point)
    {
        return new Point(
            (int)((pixel_size.Width-1)*(point.X/actual_size.Width+0.5f)),
            (int)((pixel_size.Height-1)*(point.Y/actual_size.Height+0.5f)));
    }
    /// <summary>
    /// Get the ratio of intensities included inside the strip (yellow lines). 
    /// Assume beam is located at the center.
    /// </summary>
    /// <param name="strip_width">The strip width in physical dimensions</param>
    /// <param name="strip_angle">The strip angle in degrees</param>
    /// <returns></returns>
    public float GetRatio(float strip_width, float strip_angle)
    {
        // Convert degrees to radians
        strip_angle*=(float)(Math.PI/180);
        // Cache sin() and cos()
        float cos=(float)Math.Cos(strip_angle), sin=(float)Math.Sin(strip_angle);
        //line through (xc,yc) with angle ψ is  (y-yc)*COS(ψ)-(x-xc)*SIN(ψ)=0
        //to offset the line by amount d, by add/subtract d from the equation above
        float inside=0, all=0;
        for(int i=0; i<pixel_size.Width; i++)
        {
            for(int j=0; j<pixel_size.Height; j++)
            {
                Point pixel=new Point(i, j);
                //find distance to strip center line
                PointF point=PixelToPoint(pixel);
                float t=-sin*point.X+cos*pixel.Y;
                if(Math.Abs(t)<=strip_width/2)
                {
                    inside+=intensities[i, j];
                }
                all += intensities[i,j];
            }
        }
        return inside/all;
    }
    public void RotateIntesities(float angle)
    {
        // Convert degrees to radians
        angle*=(float)(Math.PI/180);
        // Cache sin() and cos()
        float cos=(float)Math.Cos(angle), sin=(float)Math.Sin(angle);
        int[,] result=new int[pixel_size.Width, pixel_size.Height];
        for(int i=0; i<pixel_size.Width; i++)
        {
            for(int j=0; j<pixel_size.Height; j++)
            {
                Point pixel=new Point(i, j);
                PointF point=PixelToPoint(pixel);
                PointF point2=new PointF(
                    point.X*cos-point.Y*sin,
                    pixel.X*sin+point.Y*cos);
                Point pixel2=PointToPixel(point2);
                if(pixel2.X>=0&&pixel2.X<pixel_size.Width
                    &&pixel2.Y>=0&&pixel2.Y<pixel_size.Height)
                {
                    result[pixel2.X, pixel2.Y]+=intensities[i, j];
                }
            }
        }

        intensities=result;
    }
}

make sure the intesities are all positive (to avoid possible 1/0 condition). The strip center line is assumed to pass through center of beam.

Answer 3

I have an interesting idea that may be a bit different than what you have above. Here is the algorithm:

• Create two variables for the sums: betweenSum and outsideSum
• Get the equation for both lines (y = mx+b)
  • m = tan((-1)theta) where theta is the angle of the line w/ x axis. b = y intercept
• Go through each entry in the 2D array.
  • if the entry's x and y coordinates put in between both equations, then add the intensity to the betweenSum.
  • else add intensity to the outsideSum.

This approach is quite a bit different, and I'm not sure if this is something that is practical for your project. Just offering an approach that may help avoid doing rotation (and maybe any extra overhead).