Convert a fisheye image to panorama view

gistfile1.cpp

#include <cmath>
#include <iostream>
#include <opencv2/opencv.hpp>

using namespace cv;

int main(int argc, const char * argv[]) {
    Mat img = imread("lillestromfisheye.jpg");
    
    // assume the source image is square, and its width has even number of pixels
    int l = img.cols / 2;
    
    Mat destination = Mat(l, 4 * l, CV_8UC3);
    Mat bilinear = Mat(l, 4 * l, CV_8UC3);
    int i, j;
    int x, y;
    double radius, theta;
    
    // 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
    Vec3b 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;
    
    for (i = 0; i < destination.rows; ++i)
    {
        for (j = 0; j < destination.cols; ++j)
        {
            radius = (double)(l - i);
            theta = 2.0 * M_PI * (double)(4.0 * l - j) / (double)(4.0 * l);
            //theta = 2.0 * M_PI * (double)(-j) / (double)(4.0 * l);
            
            fTrueX = radius * cos(theta);
            fTrueY = radius * sin(theta);
            
            // "normal" mode
            x = (int)(round(fTrueX)) + l;
            y = l - (int)(round(fTrueY));
            // check bounds
            if (x >= 0 && x < (2 * l) && y >= 0 && y < (2 * l))
            {
                //bmDestination.SetPixel(j, i, bm.GetPixel(x, y));
                destination.at<Vec3b>(i, j) = img.at<Vec3b>(y, x);
            }
            
            // bilinear mode
            fTrueX = fTrueX + (double)l;
            fTrueY = (double)l - fTrueY;
            
            iFloorX = (int)(floor(fTrueX));
            iFloorY = (int)(floor(fTrueY));
            iCeilingX = (int)(ceil(fTrueX));
            iCeilingY = (int)(ceil(fTrueY));
            
            // check bounds
            if (iFloorX < 0 || iCeilingX < 0 ||
                iFloorX >= (2 * l) || iCeilingX >= (2 * l) ||
                iFloorY < 0 || iCeilingY < 0 ||
                iFloorY >= (2 * l) || iCeilingY >= (2 * l)) continue;
            
            fDeltaX = fTrueX - (double)iFloorX;
            fDeltaY = fTrueY - (double)iFloorY;
            
            clrTopLeft = img.at<Vec3b>(iFloorY, iFloorX);
            clrTopRight = img.at<Vec3b>(iFloorY, iCeilingX);
            clrBottomLeft = img.at<Vec3b>(iCeilingY, iFloorX);
            clrBottomRight = img.at<Vec3b>(iCeilingY, iCeilingX);
            
            // linearly interpolate horizontally between top neighbours
            fTopRed = (1 - fDeltaX) * clrTopLeft.val[2] + fDeltaX * clrTopRight.val[2];
            fTopGreen = (1 - fDeltaX) * clrTopLeft.val[1] + fDeltaX * clrTopRight.val[1];
            fTopBlue = (1 - fDeltaX) * clrTopLeft.val[0] + fDeltaX * clrTopRight.val[0];
            
            // linearly interpolate horizontally between bottom neighbours
            fBottomRed = (1 - fDeltaX) * clrBottomLeft.val[2] + fDeltaX * clrBottomRight.val[2];
            fBottomGreen = (1 - fDeltaX) * clrBottomLeft.val[1] + fDeltaX * clrBottomRight.val[1];
            fBottomBlue = (1 - fDeltaX) * clrBottomLeft.val[0] + fDeltaX * clrBottomRight.val[0];
            
            // 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;
            
            Vec3b color;
            color.val[0] = iBlue;
            color.val[1] = iGreen;
            color.val[2] = iRed;
            bilinear.at<Vec3b>(i, j) = color;
        }
    }
    
    imshow("Landscape", destination);
    imshow("Bilinear", bilinear);
    waitKey();
    
    img.release();
    destination.release();
    bilinear.release();
    
    return 0;
}