Circular Fisheye Image dewarp to flat image

Question

UPDATE as on 12 Nov 2015

I used PanoTools plugin with Photoshop and Hugin and played with all those parameters. End up i found the parameters for projection, HFOV and image output size that fulfill my lowest requirement.

Parameteres:

Processed Output:

My question is then how can i convert all these parameters and values into C# algorithm coding so that when I provide the original image, i will get the corrected output image?

Thanks a lot.

---

I have a square image captured from a circular fisheye camera. The size is 2650 * 2650 pixels.

Now, i will need to programmatically dewarp the image to a flat panorama image using C# language. I had look around from internet with different algorithm example from Link for code below , Link1 and Link2 but just can't make it success. My maths sincerely sucks and can't help me with that. Hopefully someone able to guide me through this. Thanks a lot.

Example of image output from the camera:

--Image grabbed from Wikipedia Fisheye Lens & size modified to fit my sample pixel.

The code i tried to dewarp it but no luck:

        Bitmap sourceImage = (Bitmap)Bitmap.FromFile("circularfisheye.jpg");
        double factor = 0.5;
        Boolean autoCrop = false;
        Color backgroundColor = Color.White;

        Bitmap StartImage = null;
        BitmapData srcBitmapData = null;
        Byte[] srcPixels = null;
        Byte[] dstPixels = null;
        Bitmap NewImage = null;
        BitmapData dstBitmapData = null;

        try
        {

            // Checks whether bpp ​​( Bits Per Pixel ) is 8 , 24, or 32
            int Depth = System.Drawing.Bitmap.GetPixelFormatSize(sourceImage.PixelFormat);
            if (Depth != 8 && Depth != 24 && Depth != 32)
            {
                throw new ArgumentException("Only 8, 24 and 32 bpp images are supported.");
            }

            // Retrieves the count of the color components
            int cCount = Depth / 8;

            Size baseSize = new Size(sourceImage.Width, sourceImage.Height);

            // check if a low image resize and need to improve the quality
            // and not generate image aliasing
            Int32 maxSize = Math.Max(sourceImage.Width, sourceImage.Height);
            if (maxSize < 3000)
            {
                float percent = 3000F / (float)maxSize;
                baseSize = new Size((Int32)((float)sourceImage.Width * percent), (Int32)((float)sourceImage.Height * percent));
            }

            StartImage = new Bitmap(baseSize.Width, baseSize.Height, sourceImage.PixelFormat);
            StartImage.SetResolution(sourceImage.HorizontalResolution, sourceImage.VerticalResolution);

            // Create the drawing object and white background
            Graphics g = Graphics.FromImage(StartImage);
            g.SmoothingMode = SmoothingMode.AntiAlias;
            g.InterpolationMode = InterpolationMode.HighQualityBicubic;
            g.PixelOffsetMode = PixelOffsetMode.HighQuality;
            g.DrawImage(sourceImage, new Rectangle(-1, -1, baseSize.Width + 1, baseSize.Height + 1), 0, 0, sourceImage.Width, sourceImage.Height, GraphicsUnit.Pixel);
            g.Dispose();
            // Locks the source image and copies it to the byte array and releases the source image
            srcBitmapData = StartImage.LockBits(new Rectangle(0, 0, StartImage.Width, StartImage.Height), ImageLockMode.ReadOnly, StartImage.PixelFormat);
            srcPixels = new byte[StartImage.Width * StartImage.Height * (Depth / 8)];
            Marshal.Copy(srcBitmapData.Scan0, srcPixels, 0, srcPixels.Length);
            StartImage.UnlockBits(srcBitmapData);
            srcBitmapData = null;

            // Create the target image byte array
            dstPixels = new Byte[srcPixels.Length];

            // Fill the entire frame with the selected background color
            Int32 index = ((1 * StartImage.Width) + 1) * cCount; //index = ((Y * Width) + X) * cCount
            do
            {
                if (Depth == 32) //For 32 bpp defines Red , Green, Blue and Alpha
                {
                    dstPixels[index++] = backgroundColor.B;
                    dstPixels[index++] = backgroundColor.G;
                    dstPixels[index++] = backgroundColor.R;
                    dstPixels[index++] = backgroundColor.A; // a
                }
                if (Depth == 24) //For 24 bpp defines Red , Green and Blue
                {
                    dstPixels[index++] = backgroundColor.B;
                    dstPixels[index++] = backgroundColor.G;
                    dstPixels[index++] = backgroundColor.R;
                }
                if (Depth == 8)
                // For 8 bpp defines the value of color ( Red , Green and Blue to be the same thing)
                {
                    dstPixels[index++] = backgroundColor.B;
                }

            } while (index < srcPixels.Length);
            // Calculate the maximum possible extent for the image and multiply by the desired factor
            double amp = 0;
            double ang = Math.PI * 0.5;
            for (Int32 a = 0; a < StartImage.Height; a++)
            {
                int y = (int)((StartImage.Height / 2) - amp * Math.Sin(ang));
                if ((y < 0) || (y > StartImage.Height))
                    break;
                amp = a;
            }
            amp = (amp - 2) * (factor < -1 ? -1 : (factor > 1 ? 1 : factor));
            // Define variables that calculates the cutoff points (if any)
            Int32 x1, y1, x2, y2;
            x1 = StartImage.Width;
            y1 = StartImage.Height;
            x2 = 0;
            y2 = 0;


            // Copy pixel by pixel for the new positions
            index = ((1 * StartImage.Width) + 1) * cCount;
            do
            {

                Int32 y = (Int32)((index / cCount) / StartImage.Width);
                Int32 x = (index / cCount) - (y * StartImage.Width);

                Point pt = NewPoint(new Point(x, y), StartImage.Width, StartImage.Height, amp, factor < 0);

                //Values ​​for crop
                if (factor >= 0)
                {
                    if (x == StartImage.Width / 2)
                    {
                        if (pt.Y < y1)
                            y1 = pt.Y;

                        if (pt.Y > y2)
                            y2 = pt.Y;
                    }

                    if (y == StartImage.Height / 2)
                    {
                        if (pt.X < x1)
                            x1 = pt.X;

                        if (pt.X > x2)
                            x2 = pt.X;
                    }
                }
                else
                {
                    if ((x == 1) && (y == 1))
                    {
                        y1 = pt.Y;
                        x1 = pt.X;
                    }

                    if ((x == StartImage.Width - 1) && (y == StartImage.Height - 1))
                    {
                        y2 = pt.Y;
                        x2 = pt.X;
                    }
                }

                //Bytes Index which will apply the pixel
                Int32 dstIndex = ((pt.Y * StartImage.Width) + pt.X) * cCount;

                if (Depth == 32)
                {
                    dstPixels[dstIndex] = srcPixels[index++];
                    dstPixels[dstIndex + 1] = srcPixels[index++];
                    dstPixels[dstIndex + 2] = srcPixels[index++];
                    dstPixels[dstIndex + 3] = srcPixels[index++]; // a
                }
                if (Depth == 24)
                {
                    dstPixels[dstIndex] = srcPixels[index++];
                    dstPixels[dstIndex + 1] = srcPixels[index++];
                    dstPixels[dstIndex + 2] = srcPixels[index++];
                }
                if (Depth == 8)
                {
                    dstPixels[dstIndex] = srcPixels[index++];
                }

            } while (index < srcPixels.Length);

            //Creates a new image based on the byte array previously created
            NewImage = new Bitmap(StartImage.Width, StartImage.Height, StartImage.PixelFormat);
            NewImage.SetResolution(StartImage.HorizontalResolution, StartImage.VerticalResolution);
            dstBitmapData = NewImage.LockBits(new Rectangle(0, 0, StartImage.Width, StartImage.Height), ImageLockMode.WriteOnly, StartImage.PixelFormat);
            Marshal.Copy(dstPixels, 0, dstBitmapData.Scan0, dstPixels.Length);
            NewImage.UnlockBits(dstBitmapData);


            //Generates the final image to crop or resize the real coo
            Bitmap FinalImage = new Bitmap(sourceImage.Width + 1, sourceImage.Height, StartImage.PixelFormat);
            NewImage.SetResolution(StartImage.HorizontalResolution, StartImage.VerticalResolution);

            Graphics g1 = Graphics.FromImage(FinalImage);
            g1.SmoothingMode = SmoothingMode.AntiAlias;
            g1.InterpolationMode = InterpolationMode.HighQualityBicubic;
            g1.PixelOffsetMode = PixelOffsetMode.HighQuality;

            //Performs the cut if enabled automatic cutting and there is need to cut
            if ((autoCrop) && ((x1 > 0) || (y1 > 0) || (x2 < NewImage.Height) || (y2 < NewImage.Height)))
            {
                Rectangle cropRect = new Rectangle(x1, y1, x2 - x1, y2 - y1);
                g1.DrawImage(NewImage, new Rectangle(-1, -1, FinalImage.Width + 1, FinalImage.Height + 1), cropRect.X, cropRect.Y, cropRect.Width, cropRect.Height, GraphicsUnit.Pixel);
            }
            else
            {
                g1.DrawImage(NewImage, new Rectangle(-1, -1, FinalImage.Width + 1, FinalImage.Height + 1), 0, 0, NewImage.Width, NewImage.Height, GraphicsUnit.Pixel);
            }

            g1.Dispose();
            g1 = null;

            NewImage = null;
            FinalImage.Save("output.jpg");
            FinalImage.Dispose();
        }
        finally
        {
            srcBitmapData = null;
            srcPixels = null;
            dstPixels = null;
            dstBitmapData = null;
        }

Answer 1

Such a distortion as a symmetry of revolution.

In polar coordinates, with the pole at the center of the image, it is expressed as

r' = f(r)
Θ' = Θ

where the quote indicates the distorted coordinates. The function f is unknown and should be measured empirically, by calibration (looking at a regular target).

To correct the image, you need to invert the function f and apply the reverse transform to the image. In fact, it is easier to measure g directly by calibration. As a starting approximation, a simple model like

r = r' + a.r'³ 

can do.

Most probably you don't have a picture of a grid taken with the same lens. Your last resort is to implement the undistortion function with adjustable parameters, and optimize these by trial and error.

It should also be possible to derive the calibration curve by looking at the deformation of straight lines, but this is more "technical".

---

In Cartesian coordinates, you can express the correction transform as

x = g(r').x'/r'
y = g(r').y'/r'

where r' = √x'²+y'².

Answer 2

Use the algorithm from here:

http://www.helviojunior.com.br/fotografia/barrel-and-pincushion-distortion/

It worked for me

Answer 3

I've made some revamp to the HelvioJunior's library (that was linked by @Tarek.Mh), I think this may suit your need:

Below, the code:

using System.Drawing;
using System.Drawing.Imaging;
using System.Linq;
using System.Runtime.InteropServices;
using static System.Math;

namespace HelvioJunior
{
    //https://www.helviojunior.com.br/fotografia/barrel-and-pincushion-distortion/
    public class Program
    {
        private static void Main(string[] args)
        {
            Bitmap source = (Bitmap)Image.FromFile(@"JpwX0.png");
            Bitmap bmp = BarrelDistortion(source, 4/10f, true);
            bmp.Save(@"test.png");
            bmp.Dispose();
            source.Dispose();
        }

        static public Bitmap BarrelDistortion(Bitmap sourceImage, double factor = 0, bool autoCrop = true, uint previewRectangleWidth = 0, Color? fillerColor = null)
        {
            int sourceRight = sourceImage.Width - 1, sourceBottom = sourceImage.Height - 1;

            // Vertical amplitude is half the height times factor
            // Horizontal amplitude is missing ; vertical amplitude's applied to both directions
            double amp = sourceBottom / 2f * factor;

            // Inner shrinking area points
            RePoint[] lPts;
            bool inverse = factor < 0;

            // Shrinking area coordinates (center point is considered always available)
            double x1 = sourceRight / 2f,
                y1 = sourceBottom / 2f,
                x2 = sourceRight / 2f,
                y2 = sourceBottom / 2f;

            if (inverse)
            {
                lPts = new RePoint[]
                {
                    new RePoint(0, 0),
                    new RePoint(0, sourceBottom),
                    new RePoint(sourceRight, sourceBottom),
                    new RePoint(sourceRight, 0)
                };
            }
            else
            {
                lPts = new RePoint[]
                {
                    new RePoint(sourceRight * 1 / 2f, 0),
                    new RePoint(0, sourceBottom * 1 / 2f),
                    new RePoint(sourceRight, sourceBottom * 1 / 2f),
                    new RePoint(sourceRight * 1 / 2f, sourceBottom)
                };
            }

            foreach (var pN in lPts.Select(pt => NewPoint(pt, sourceImage.Width, sourceImage.Height, amp, inverse)))
            {
                if (pN.Y < y1) y1 = pN.Y;
                if (pN.Y > y2) y2 = pN.Y;
                if (pN.X < x1) x1 = pN.X;
                if (pN.X > x2) x2 = pN.X;
            }

            // Bytes per color from bit per pixel (bpp) format
            int bpcCount = Image.GetPixelFormatSize(sourceImage.PixelFormat) / 8;

            Rectangle sourceRectangle = new Rectangle(0, 0, sourceImage.Width, sourceImage.Height);
            int srcLength = sourceImage.Width * sourceImage.Height * bpcCount;

            // Gets sourceImage byte array as srcpixels
            BitmapData srcBitmapData = sourceImage.LockBits(sourceRectangle, ImageLockMode.ReadOnly, sourceImage.PixelFormat);
            byte[] srcPixels = new byte[srcLength];
            Marshal.Copy(srcBitmapData.Scan0, srcPixels, 0, srcLength);
            sourceImage.UnlockBits(srcBitmapData);
            srcBitmapData = null;

            // Destination byte array preparation as dstPixels
            byte[] dstPixels = new byte[srcLength];
            int dstIndex = 0;

            // Filler color preparation
            Color fillColor = fillerColor ?? Color.Transparent;
            if (!autoCrop)
            {
                if (bpcCount <= 4) // Depth > 32bpp may not work as expected, filler color's not applied for bit safety reason
                    do
                    {
                        dstPixels[dstIndex++] = fillColor.B;
                        if (bpcCount > 1)
                        {
                            dstPixels[dstIndex++] = fillColor.G;
                            dstPixels[dstIndex++] = fillColor.R;
                            if (bpcCount > 3)
                                dstPixels[dstIndex++] = fillColor.A; // a
                        }
                    } while (dstIndex < srcLength);
            }

            // Byte-to-byte copy (incl. Point transformation)
            int index = 0, srcBpcLength = srcLength - bpcCount;
            do
            {
                int comp = index / bpcCount; // comp yields the current "pixel" position
                int y = comp / sourceImage.Width; // Each line is sourceImage.Width bytes wide
                int x = comp - (y * sourceImage.Width); // Remaining (comp - lines) bytes is target column (ranges from 0 to width - 1)

                // Destination "pixel"
                RePoint pt = NewPoint(new RePoint(x, y), sourceImage.Width, sourceImage.Height, amp, inverse);
                dstIndex = (((int)pt.Y * sourceImage.Width) + (int)pt.X) * bpcCount; // dstIndex++ overflows when |amp| >= 2

                if (dstIndex >= 0 && dstIndex <= srcBpcLength)
                    for (int i = 0; i++ < bpcCount;)
                        dstPixels[dstIndex++] = srcPixels[index++];
                else
                    index += bpcCount;
            } while (index < srcLength);
            srcPixels = null;

            // Destination bytes application
            BitmapData dstBitmapData = sourceImage.LockBits(sourceRectangle, ImageLockMode.WriteOnly, sourceImage.PixelFormat);
            Marshal.Copy(dstPixels, 0, dstBitmapData.Scan0, srcLength);
            sourceImage.UnlockBits(dstBitmapData);
            dstBitmapData = null;
            dstPixels = null;

            // Final Image area
            Rectangle cropRect = new Rectangle((int)Ceiling(x1), (int)Ceiling(y1), (int)Ceiling(x2 - x1), (int)Ceiling(y2 - y1));
            Rectangle destRectangle = autoCrop ? cropRect : sourceRectangle;

            // Final image preparation
            Bitmap FinalImage = new Bitmap(destRectangle.Width, destRectangle.Height, sourceImage.PixelFormat);
            FinalImage.SetResolution(sourceImage.HorizontalResolution, sourceImage.VerticalResolution);

            Graphics g1 = Graphics.FromImage(FinalImage);
            g1.DrawImage(sourceImage, -destRectangle.X, -destRectangle.Y);

            // Previsualization rectangle
            if (previewRectangleWidth > 0)
                g1.DrawRectangle(new Pen(Color.Red, previewRectangleWidth), cropRect.X - 1, cropRect.Y - 1, cropRect.Width + previewRectangleWidth, cropRect.Height + previewRectangleWidth);

            g1.Dispose();
            g1 = null;

            return FinalImage;
        }

        private static RePoint NewPoint(RePoint aP, double Width, double Height, double Amplitude, bool inverse)
        {
            double h = aP.Y / (Height - 1);
            double w = aP.X / (Width - 1);

            // Works ok for [0/2] to [1/2]
            // Floating point error(s) here, in the range of ]1/2] to [2/2] (No workaround found)
            double sinX = Round(Sin(PI * w), 15); // Range of [0] to [1] * PI ; result ranges from 0 (far from center) to 1 (at center)
            double sinY = Round(Sin(PI * h), 15);

            double caX = Amplitude * (1 - 2 * w);
            double caY = Amplitude * (1 - 2 * h);

            double aY = 0, aX = 0;
            if (inverse)
            {
                aX = -caX;
                aY = -caY;
            }


            double pY = aP.Y + aY + caY * sinX;
            double pX = aP.X + aX + caX * sinY;

            return new RePoint(pX, pY);
        }

        private struct RePoint
        {
            public double X;
            public double Y;

            public RePoint(double x, double y)
            {
                X = x;
                Y = y;
            }
        }
    }
}