implement openCV method warpPerspective()

Question

I am trying to implement openCV method warpPerspective() from scratch, I made the code below, it can handle shifts in y and x but, when I pass homography matrix from findHomography() to the function I made it always gives blank image compared to warpPerspective() output.

I followed this definition to find the new locations of the pixels:

s*x' h1 h2 h3 x s*y' = h4 h5 h6 * y s h7 h8 1 1

my mapping works on simple shift like { 1, 0.5,-51,0,1,50,0,0,1}

but when matrix is like: [1.0340946, 0.032195676, -6.419126; 0.00302419, 1.0487343, -96.520393; 3.7013847e-06, 0.00010837225, 1] the output is like this:

my implementation: -Given H, and image A, -Find new locations of pixels in A and save them in TransArry.where the index of the array is the linearized index of A. -Remap pixels of A into tranImg.

    Mat transform(Mat A, Mat H)
{
// allocate array of all locations
int Numrows = A.rows;
int Numcols = A.cols;
int channels   = A.channels();
cout << "rows " << Numrows << "col " << Numcols << "channels " << channels <<endl;
int size = Numrows*Numcols;
int MaxX,MaxY = -1000;
int MinX,MinY =  1000;
int *TransArry = (int *)malloc(sizeof(int)*size);
int Idx;

 int homeX=Idx % Numcols;
 int homeY=Idx / Numcols;
 cout << H << endl;

 waitKey();         
for (Idx=0; Idx < size; ++Idx ){

        homeX=Idx % Numcols;
        homeY=Idx / Numcols;

        float x  = (H.at<float>(0,0) * (homeX)) +( H.at<float>(0,1) * (homeY)) + ( H.at<float>(0,2) * 1) ;
        float y  = (H.at<float>(1,0) * (homeX)) +( H.at<float>(1,1) * (homeY)) + ( H.at<float>(1,2) * 1) ;
        float s  = (H.at<float>(2,0) * (homeX)) +( H.at<float>(2,1) * (homeY)) + ( H.at<float>(2,2) * 1) ;

        cout << " x = " << x << " y= " << y << " s= " << s;
        x = (x/s);

        y = y/s;

        // for the first col in TransMatrix
        if (homeX ==0){
            if (x > MaxX) MaxX = x;
            if (x < MinX) MinX = x; 
        }

        //for thee first row in TransMatrix
        if (homeY ==0){
            if (y > MaxY) MaxY = y;
            if (y < MinY) MinY = y;
        }
        if((y)>=A.rows || (y)<0 || (x)>=A.cols || (x)<0){
            TransArry[Idx]  = -1;
            cout << "x= " << x << "y= "<< y << endl;
        }else{
            TransArry[Idx] = (y * Numcols + x); 
        }           

        //cout << Numcols << endl;
        cout <<     "New index of " << Idx << "is " << TransArry[Idx] << endl;
        }

         Mat   tranImg ;

         A.copyTo(tranImg);
         tranImg = tranImg - tranImg;
        cout <<     "Rows" << tranImg.rows << "cols" << tranImg.cols << "cha" <<  A.channels() << endl;


        waitKey();
        // Remap Image
        for (Idx=0; Idx < size; Idx ++ ){

            homeX=Idx % Numcols;
            homeY=Idx / Numcols;                
            //tranImg.at<uchar>(homeY, homeX) =0;
            if(TransArry[Idx] != -1){   
                //cout << "Index " << Idx << "Passed " << endl;
                int newhomeX=TransArry[Idx] % Numcols; // Col ID
                int newhomeY=TransArry[Idx] / Numcols;  // Row ID


                 cout << "Index is " << Idx << endl;
                 cout << "HomeX is " << homeX << " and HomeY is " << homeY << endl;
                 cout << "New Index is " << TransArry[Idx] << endl;
                 cout << "New HomeX is " << newhomeX << " and New HomeY is " << newhomeY << endl;   
                 cout << "*****************************************"<< endl; 
                // if (!(Idx%100)) sleep(20);  

                tranImg.at<uchar>(newhomeY, (newhomeX*channels)) = A.at<uchar>(homeY, homeX*channels);
                if(channels>1)
                    tranImg.at<uchar>(newhomeY, newhomeX*channels+1) = A.at<uchar>(homeY, homeX*channels+1);
                if(channels>2)
                    tranImg.at<uchar>(newhomeY, newhomeX*channels+2) = A.at<uchar>(homeY, homeX*channels+2);
                // if (!(Idx%100)){
                    // imshow("inside", tranImg);
                    // waitKey(1);
                    // }
                }
        }
        //cout << tranImg << endl;  

return tranImg;

}

H is calaculated and verified.

So, is there a problem in the way I access the matrices H and A?

Answer 1

I found the problem. The co-ordinates you are using are floating point values instead of integers. Typecast x and y to integers and floor the values x = floor(x/s) and y =floor(y/s).

Answer 2

See my answer here and it should answer both of your questions. This answer still uses warpPerspective() but it gives the full insight into how to calculate how far out of the image bounds your homography warps the image, so you can pad appropriately on each side of the image.

As for doing warpPerspective() manually, all you need to do is put all your image coordinates into a linearized homogeneous array and multiply by your homography, and then divide by the last coordinate to get back to Cartesian coordinates. Then, you can use remap() to do the interpolation. The syntax for remap() can be confusing, so you can reference my answer here to get an idea of how to use it. This answer shows a manual transformation and interpolation, it should give you (at least nearly) identical results to warpPerspective().