c++ - How can I rewrite this warp-affine using OpenCV?

Question

I'm trying to optimize this code, in particular:

bool interpolate(const Mat &im, float ofsx, float ofsy, float a11, float a12, float a21, float a22, Mat &res)
{         
   bool ret = false;
   // input size (-1 for the safe bilinear interpolation)
   const int width = im.cols-1;
   const int height = im.rows-1;
   // output size
   const int halfWidth  = res.cols >> 1;
   const int halfHeight = res.rows >> 1;
   float *out = res.ptr<float>(0);
   for (int j=-halfHeight; j<=halfHeight; ++j)
   {
      const float rx = ofsx + j * a12;
      const float ry = ofsy + j * a22;
      for(int i=-halfWidth; i<=halfWidth; ++i)
      {
         float wx = rx + i * a11;
         float wy = ry + i * a21;
         const int x = (int) floor(wx);
         const int y = (int) floor(wy);
         if (x >= 0 && y >= 0 && x < width && y < height)
         {
            // compute weights
            wx -= x; wy -= y;
            // bilinear interpolation
            *out++ = 
               (1.0f - wy) * ((1.0f - wx) * im.at<float>(y,x)   + wx * im.at<float>(y,x+1)) +
               (       wy) * ((1.0f - wx) * im.at<float>(y+1,x) + wx * im.at<float>(y+1,x+1));
         } else {
            *out++ = 0;
            ret =  true; // touching boundary of the input            
         }
      }
   }
   return ret;
}

According to Intel Advisor, this is a very time consuming function. In this question I asked how I could optimize this, and someone made me notice that this is warp-affine transformation.

Now, since I'm not the image processing guy, I had to read this article to understand what a warp-affine transformation is.

To my understanding, given a point p=(x,y), you apply a transformation A (a 2x2 matrix) and then translate it by a vector b. So the obtained point after the transformation p' can be expressed as p' = A*p+b. So far so good.

However, I'm a little bit confused on how to apply cv::warpAffine() to this case. First of all, from the function above interpolate() I can see only the 4 A components (a11, a12, a21, a22), while I can't see the 2 b components...Are they ofsx and ofy?

In addition notice that this function returns a bool value, which is not returned by warpAffine (this boolean value is used here at line 126), so I don't know I could this with the OpenCV function.

But most of all I'm so confused by for (int j=-halfHeight; j<=halfHeight; ++j) and for(int i=-halfWidth; i<=halfWidth; ++i) and all the crap that happens inside.

I understand that:

        // bilinear interpolation
        *out++ = 
           (1.0f - wy) * ((1.0f - wx) * im.at<float>(y,x)   + wx * im.at<float>(y,x+1)) +
           (       wy) * ((1.0f - wx) * im.at<float>(y+1,x) + wx * im.at<float>(y+1,x+1));

Is what INTER_LINEAR does, but apart from that I'm totally lost.

So, to test my approach, I tried to do the equivalent of line 131 of this as:

     bool touchesBoundary = interpolate(smoothed, (float)(patchImageSize>>1), (float)(patchImageSize>>1), imageToPatchScale, 0, 0, imageToPatchScale, patch);
     Mat warp_mat( 2, 3, CV_32FC1 );
     float a_11 = imageToPatchScale;
     float a_12 = 0;
     float a_21 = 0;
     float a_22 = imageToPatchScale;
     float ofx = (float)(patchImageSize>>1);
     float ofy = (float)(patchImageSize>>1);
     float ofx_new = ofx - a12*halfHeight - a11*halfWidth;
     float ofy_new = ofy - a22*halfHeight - a21*halfWidth;
     warp_mat.at<float>(0,0) = imageToPatchScale;
     warp_mat.at<float>(0,1) = 0;
     warp_mat.at<float>(0,2) = ofx_new;
     warp_mat.at<float>(1,0) = 0;
     warp_mat.at<float>(1,1) = imageToPatchScale;
     warp_mat.at<float>(1,2) = ofy_new;
     cv::Mat myPatch;
     std::cout<<"Applying warpAffine"<<std::endl;
     warpAffine(smoothed, myPatch, warp_mat, patch.size());
     std::cout<<"WarpAffineApplied patch size="<<patch.size()<<" myPatch size="<<myPatch.size()<<std::endl;
     cv::Mat diff = patch!=myPatch;
     if(cv::countNonZero(diff) != 0){
         throw std::runtime_error("Warp affine doesn't work!");
     }
     else{
         std::cout<<"It's working!"<<std::endl;
     }

And of course at the first time the this is executed, the exception is thrown (so the two methods are not equivalent)...How can I solve this?

Can someone help me please?

As I already written in the comments, the resulting matrix by using the code above is a zero matrix. While this is maPatch obtained by using ofx and ofy instead of ofx_new and ofy_new, while patch has all the values different from zero:

[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 229.78679, 229.5752, 229.11732, 229.09612, 229.84615, 230.28633, 230.35257, 230.70955, 230.99368, 231.00777, 231.20511, 231.63196, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 230.60367, 230.16417, 230.07034, 230.06793, 230.02016, 230.14925, 230.60413, 230.84822, 230.92368, 231.02249, 230.99162, 230.9149, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 232.76547, 231.39716, 231.26674, 231.34512, 230.746, 230.25253, 229.65276, 227.83998, 225.43642, 229.57695, 230.31363, 230.16011, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 234.01663, 232.88118, 232.15475, 231.40129, 223.21553, 208.22626, 205.58975, 214.53882, 220.32681, 228.11552, 229.31509, 228.86545, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 234.04565, 233.00443, 231.9902, 230.14912, 198.0849, 114.86175, 97.901344, 160.0218, 217.38528, 231.07045, 231.13109, 231.10185, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 233.293, 232.69095, 217.03873, 190.56714, 167.61592, 94.968391, 81.302032, 150.72263, 194.79535, 215.15564, 230.01717, 232.37894, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 231.70988, 227.81319, 207.59377, 173.35149, 113.88276, 73.171112, 71.523285, 103.05875, 160.05588, 194.65132, 226.4287, 231.45871, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 231.93924, 224.24269, 199.1693, 150.65695, 103.33984, 79.489555, 77.509094, 87.893059, 122.01918, 168.37506, 219.22086, 231.05161, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 232.2706, 232.12926, 206.97635, 127.69308, 92.714355, 81.512207, 74.89402, 75.968353, 84.518105, 157.07962, 223.18773, 229.92766, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 232.64882, 222.16704, 161.95021, 92.577881, 83.757164, 76.764214, 67.041054, 66.195595, 71.112335, 131.66878, 188.27278, 217.6635, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 234.77202, 231.75511, 178.64326, 104.27015, 95.664223, 82.791382, 67.68969, 72.78054, 72.355469, 104.77696, 172.32361, 204.92691, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 236.49684, 235.5802, 185.34337, 115.96995, 106.85963, 82.980408, 61.703068, 69.540627, 76.200562, 82.429321, 101.46993, 119.75877, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

Answer 1

I'm more familiar with this warpAffine, whose basic statement is

cv::warpAffine (InputArray src,     // input mat
                OutputArray dst,    // output mat
                InputArray M,       // affine transformation mat
                Size dsize)         // size of the output mat

where M is the matrix

      a11 a12 ofx
      a21 a22 ofy

In your term, the first two columns is the linear transformation matrix A, the last is the translation vector b.

The cv::hal::warpAffine() is just the same, where double M[6] corresponds to the above affine transformation matrix, but I'm not sure in which order it is flatten (most likely, [a11,a12,ofx,a21,a22,ofy]).

---

In OpenCV, the origin (0,0) is the top-left conner as usual, while in Intel's code, the origin (0,0) is in the middle of the image. That's what the part

for (int j=-halfHeight; j<=halfHeight; ++j)
{
   for(int i=-halfWidth; i<=halfWidth; ++i)
      {
         const int y = (int) floor(wy);
         //...
   }
}

does: (i,j) is the coordinate in res, j from -halfHeight to halfHeight and i from -halfHeight to halfHeight. So in this case (0,0) is in the center of the res image.

In the provided code, if you want to map src onto res (i guess), you would need to do:

    bool touchesBoundary = interpolate(smoothed, (float)(imageSize>>1), (float)(imageSize>>1), imageToPatchScale, 0, 0, imageToPatchScale, patch);

Notice here imageSize>>1 instead of patchImageSize>>1. Why? You want the center of the res (i=0,j=0) maps to the center of src, i.e. the value src.at<float>(src.cols/2, src.rows/2) (why?)

Now to make that work in your example, the equivalent of cv::warpedAffine() would be

warpAffine(smoothed, myPatch, warp_mat, patch.size(),WARP_INVERSE_MAP);

where the warp_mat has ofsx=0,ofsy=0.

Finally, here's an illustration of what I tried:

where diff = mypatch - patch >5 and smoothed is scaled up by OS. Notice the black border in patch, it is because the restrictions x < width and y<height in the code.