Revision history [back]

This is the function. Is not well written, is not optimized, but it seems to work. Thanks to berak who helped me understand better the opencv functions.

/// code start

    double MahalanobisDistance(cv::Mat& model, cv::Mat& scene) {

std::vector<int> mean(3);
mean[0] = 0;
mean[1] = 0;
mean[2] = 0;

//cv::Mat covar = cv::Mat::zeros(3, 3, CV_32F);
std::vector<double> covar, invCovar(9);

int numberOfPixel = model.cols * model.rows;

for (int iter = 0; iter < 9; iter++)
    covar.push_back(0.0);

//mean: Sum for each channel
for (int x = 0; x < model.cols; x++)
    for (int y = 0; y < model.rows; y++)
    {
        //cout << "\n Pixel: " << int(model.at<Vec3b>(y, x)[1]) << " \n";
        mean[0] += int(model.at<Vec3b>(y, x)[0]);
        mean[1] += int(model.at<Vec3b>(y, x)[1]);
        mean[2] += int(model.at<Vec3b>(y, x)[2]);
    }

mean[0] = cvRound(double(mean[0]) / double(numberOfPixel));
mean[1] = cvRound(double(mean[1]) / double(numberOfPixel));
mean[2] = cvRound(double(mean[2]) / double(numberOfPixel));


for (int x = 0; x < model.cols; x++)
    for (int y = 0; y < model.rows; y++)
    {
        for (int i = 0; i < 3; i++)
            for (int j = 0; j <3; j++)
            {
                covar[i * 3 + j] += (int(model.at<Vec3b>(y, x)[i]) - mean[i]) * (int(model.at<Vec3b>(y, x)[j]) - mean[j]);
            }

    }

for (int iter = 0; iter < 9; iter++)
    covar[iter] = cvRound(double(covar[iter]) / double(numberOfPixel));

//invert covar
//compute det of A.
double detA = (covar[0] * covar[4] * covar[8]) + (covar[1] * covar[5] * covar[6]) + (covar[2] * covar[3] * covar[7])
- (covar[2] * covar[4] * covar[6]) - (covar[0] * covar[5] * covar[7]) - (covar[8] * covar[3] * covar[1]);

//compute det of minors
std::vector<double> minors(9);

minors[0] = (covar[4] * covar[8]) - (covar[5] * covar[7]);
minors[1] = -((covar[3] * covar[8]) - (covar[5] * covar[6]));
minors[2] = (covar[3] * covar[7]) - (covar[4] * covar[6]);

minors[3] = -((covar[1] * covar[8]) - (covar[2] * covar[7]));
minors[4] = (covar[0] * covar[8]) - (covar[2] * covar[6]);
minors[5] = -((covar[0] * covar[7]) - (covar[1] * covar[6]));

minors[6] = (covar[1] * covar[5]) - (covar[4] * covar[2]);
minors[7] = -((covar[0] * covar[5]) - (covar[2] * covar[3]));
minors[8] = (covar[0] * covar[4]) - (covar[3] * covar[1]);

//transpose and detA
invCovar[0] = minors[0] / detA;
invCovar[1] = minors[3] / detA;
invCovar[2] = minors[6] / detA;

invCovar[3] = minors[1] / detA;
invCovar[4] = minors[4] / detA;
invCovar[5] = minors[7] / detA;

invCovar[6] = minors[2] / detA;
invCovar[7] = minors[5] / detA;
invCovar[8] = minors[8] / detA;


//test
/*std::vector<double> test(9);

test[0] = (covar[0] * invCovar[0]) + (covar[1] * invCovar[3]) + (covar[2] * invCovar[6]);
test[1] = (covar[0] * invCovar[1]) + (covar[1] * invCovar[4]) + (covar[2] * invCovar[7]);
test[2] = (covar[0] * invCovar[2]) + (covar[1] * invCovar[5]) + (covar[2] * invCovar[8]);

test[3] = (covar[3] * invCovar[0]) + (covar[4] * invCovar[3]) + (covar[5] * invCovar[6]);
test[4] = (covar[3] * invCovar[1]) + (covar[4] * invCovar[4]) + (covar[5] * invCovar[7]);
test[5] = (covar[3] * invCovar[2]) + (covar[4] * invCovar[5]) + (covar[5] * invCovar[8]);

test[6] = (covar[6] * invCovar[0]) + (covar[7] * invCovar[3]) + (covar[8] * invCovar[6]);
test[7] = (covar[6] * invCovar[1]) + (covar[7] * invCovar[4]) + (covar[8] * invCovar[7]);
test[8] = (covar[6] * invCovar[2]) + (covar[7] * invCovar[5]) + (covar[8] * invCovar[8]);

for (int iter = 0; iter < 9; iter++)
{
    cout << "covar: " << covar[iter] << " , inverse: " << invCovar[iter] <<  " , mul: " << test[iter] << "\n";
}
*/

double MahabDistance = 0.0;

for (int x = 0; x < scene.cols; x++)
    for (int y = 0; y < scene.rows; y++)
    {
        std::vector<int>  Pixel(3);
        Pixel[0] = int(scene.at<Vec3b>(y, x)[0]);
        Pixel[1] = int(scene.at<Vec3b>(y, x)[1]);
        Pixel[2] = int(scene.at<Vec3b>(y, x)[2]);

        //cout << "\n Mean: " << mean[0] << " , " << mean[1] << " , " << mean[2] << " \n";
        //cout << "Pixel: " << Pixel[0] << " , " << Pixel[1] << " , " << Pixel[2] << " \n";

        for (int i = 0; i < 3; i++)
            for (int j = 0; j < 3; j++)
            {
                //cout << double(Pixel[j] - mean[j]) << " \n";
                //cout << invCovar[i * 3 + j] << " \n";
                //cout << double(Pixel[i] - mean[i]) << " \n";

                MahabDistance += double(Pixel[j] - mean[j]) * invCovar[i * 3 + j] * double(Pixel[i] - mean[i]);
                //cout << "M dist (x, y) " << x << y << ": = " << MahabDistance << "\n\n";
            }
    }

MahabDistance = sqrt(MahabDistance);
cout << "M dist: " << " = " << MahabDistance << "\n";

return MahabDistance;
    }

/// code end