I have this code:
vector<Mat> blurs (par.numberOfScales+3, Mat());
vector<float> sigmas(par.numberOfScales+3);
//fill blurs and sigmas someway
...
Mat high = hessianResponse(blurs[i+1], sigmas[i+1]*sigmas[i+1]);
Mat comp = hessianResponse(blurs[i+1], sigmas[i+1]*sigmas[i+1]);
equalMats(high, comp);
...
Where:
void equalMats(const Mat &a, const Mat &b){
cv::Mat diff;
cv::bitwise_xor(a, b, diff);
// Equal if no elements disagree
assert( cv::countNonZero(diff) == 0);
}
And (taken from here):
Mat HessianDetector::hessianResponse(const Mat &inputImage, const float norm)
{
const int rows = inputImage.rows;
const int cols = inputImage.cols;
const int stride = cols;
// allocate output
Mat outputImage(rows, cols, CV_32FC1);
// setup input and output pointer to be centered at 1,0 and 1,1 resp.
const float *in = inputImage.ptr<float>(1);
float *out = outputImage.ptr<float>(1) + 1;
float norm2 = norm * norm;
/* move 3x3 window and convolve */
for (int r = 1; r < rows - 1; ++r)
{
float v11, v12, v21, v22, v31, v32;
/* fill in shift registers at the beginning of the row */
v11 = in[-stride]; v12 = in[1 - stride];
v21 = in[ 0]; v22 = in[1 ];
v31 = in[+stride]; v32 = in[1 + stride];
/* move input pointer to (1,2) of the 3x3 square */
in += 2;
for (int c = 1; c < cols - 1; ++c)
{
/* fetch remaining values (last column) */
const float v13 = in[-stride];
const float v23 = *in;
const float v33 = in[+stride];
// compute 3x3 Hessian values from symmetric differences.
float Lxx = (v21 - 2*v22 + v23);
float Lyy = (v12 - 2*v22 + v32);
float Lxy = (v13 - v11 + v31 - v33)/4.0f;
/* normalize and write out */
*out = (Lxx * Lyy - Lxy * Lxy)*norm2;
/* move window */
v11=v12; v12=v13;
v21=v22; v22=v23;
v31=v32; v32=v33;
/* move input/output pointers */
in++; out++;
}
out += 2;
}
return outputImage;
}
If you want to see full hessianResponse's implementation, see here.
Now, for some reason, equalMats fails (so high and comp are different). A subset of different values:
a(0,652)=11.7962 b(0,652)=0
a(0,653)=13.2172 b(0,653)=0
a(0,654)=14.4621 b(0,654)=0
a(0,655)=15.4783 b(0,655)=0
a(0,656)=16.2395 b(0,656)=0
a(0,657)=16.7433 b(0,657)=0
a(0,658)=17.0053 b(0,658)=0
a(0,659)=17.0513 b(0,659)=0
a(0,660)=16.9105 b(0,660)=0
a(0,661)=16.6106 b(0,661)=0
a(0,662)=16.1759 b(0,662)=0
a(0,663)=15.6269 b(0,663)=0
a(0,664)=14.9816 b(0,664)=0
a(0,665)=14.2563 b(0,665)=0
a(0,666)=13.4675 b(0,666)=0
a(0,667)=12.6334 b(0,667)=0
a(0,668)=11.7744 b(0,668)=0
a(0,669)=10.9132 b(0,669)=0
a(0,670)=10.0744 b(0,670)=0
a(0,671)=9.28306 b(0,671)=0
a(0,672)=8.56325 b(0,672)=0
a(0,673)=7.93605 b(0,673)=0
a(0,674)=7.41722 b(0,674)=0
a(0,675)=7.01438 b(0,675)=0
a(0,676)=6.72463 b(0,676)=0
a(1,676)=16.2824 b(1,676)=0
a(2,0)=16.5866 b(2,0)=0
a(2,676)=3.11132 b(2,676)=0
a(3,0)=2.95835 b(3,0)=0
a(3,676)=10.2835 b(3,676)=0
a(4,0)=9.54264 b(4,0)=0
a(4,676)=41.1547 b(4,676)=0
a(5,0)=40.6696 b(5,0)=0
a(5,676)=45.7742 b(5,676)=0
a(6,0)=45.757 b(6,0)=0
a(6,676)=45.5739 b(6,676)=0
a(7,0)=45.5451 b(7,0)=0
a(7,676)=43.6088 b(7,676)=0
a(8,0)=44.1428 b(8,0)=0
I don't really understand how this is possible, since hessianResponse arguments are const and there is no side effect in hessianResponse.
Notice that the algorithm works correctly if there is no check, so using only Mat high = hessianResponse(blurs[i+1], par.sigmas[i+1]*par.sigmas[i+1]); with no comp or equalsMats works fine and it produce the expected result.
Why this happen?
If you wonder why I'm doing this:
I'm trying to parallelize the algorithm of the linked code, and in order to parallelize the for (int i = 1; i < par.numberOfScales+2; i++) loop, I need to pre-compute all the variables. As you can see above, I successfully replaced sigma, blur and nextBlur with sigmas[i+1], blurs[i] and blurs[i+1], respectively, introducing this code before the for that I want to parallelize:
blurs[1] = firstLevel.clone();
for (int i = 1; i < par.numberOfScales+2; i++){
float sigma = par.sigmas[i]* sqrt(sigmaStep * sigmaStep - 1.0f);
blurs[i+1] = gaussianBlur(blurs[i], sigma);
std::cout<<"rows="<<blurs[i+1].rows<<" blurs.cols="<<blurs[i+1].cols<<std::endl;
//highs[i+1] = hessianResponse(blurs[i+1], par.sigmas[i+1]*par.sigmas[i+1]);
}
As you can see, the last line is commented because of the problem above (it came out that highs[i+1] != high) and so I introduced comp only for debugging purpose.