OpenCV matrix weird behavior of addition, multiplication with float point and 8bit value

Question

I'm using OpenCV version 3.4.5. I tried to understand the behavior of cv::Mat in OpenCV. I multiple some 8bit matrices and add them up. The result seems inconsistent.

Does OpenCV cast the float point result immediately back to 8bit after each multiplication, or after each addition, or maybe some "Fused Multiply Add" happened here ?

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

int main() {
    for (int value = 1; value < 10; value++) {
        cv::Mat im(1, 1, CV_8U, value);
        cv::Mat result[10];
        result[0] = im*0.1f;
        result[1] = im*0.1f + im*0.1f;
        result[2] = im*0.1f + im*0.1f + im*0.1f;
        result[3] = im*0.1f + im*0.1f + im*0.1f + im*0.1f;
        result[4] = im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f;
        result[5] = im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f;
        result[6] = im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f;
        result[7] = im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f;
        result[8] = im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f;
        result[9] = im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f + im*0.1f;
        std::cout << "base value: " << value << std::endl;
        for (int i = 0; i < 10; i++) {
            std::cout << i + 1 << ": " << int(result[i].at<uint8_t>(0, 0)) << "\t";
        }
        std::cout << std::endl;
    }
    return 0;
}

The outputs are:

base value: 1
1: 0    2: 0    3: 0    4: 0    5: 0    6: 0    7: 0    8: 0    9: 0    10: 0
base value: 2
1: 0    2: 0    3: 0    4: 0    5: 0    6: 0    7: 0    8: 0    9: 0    10: 0
base value: 3
1: 0    2: 1    3: 1    4: 1    5: 1    6: 1    7: 1    8: 1    9: 1    10: 1
base value: 4
1: 0    2: 1    3: 1    4: 1    5: 1    6: 1    7: 1    8: 1    9: 1    10: 1
base value: 5
1: 0    2: 1    3: 2    4: 2    5: 2    6: 2    7: 2    8: 2    9: 2    10: 2
base value: 6
1: 1    2: 1    3: 2    4: 3    5: 4    6: 5    7: 6    8: 7    9: 8    10: 9
base value: 7
1: 1    2: 1    3: 2    4: 3    5: 4    6: 5    7: 6    8: 7    9: 8    10: 9
base value: 8
1: 1    2: 2    3: 3    4: 4    5: 5    6: 6    7: 7    8: 8    9: 9    10: 10
base value: 9
1: 1    2: 2    3: 3    4: 4    5: 5    6: 6    7: 7    8: 8    9: 9    10: 10

Answer 1

In OpenCV, the operations between cv::Mat is not evaluated right away for performance reason. It uses delayed evaluation. If you do A = B*0.1 + C*0.1 (A,B,C is cv::Mat), a cv::MatExpr that holds B,C and their scale value 0.1 and 0.1 is created. And it is evaluated only when the operations are more than 3 or assigned to an actual cv::Mat.

If you don't assign B*0.1 + C*0.1 to an actual cv::Mat, the evaluation will not happen.

If you do A = B*0.1 + C*0.1 + D*0.1, a cv::MatExpr that holds B*0.1 and C*0.1is created, and then another cv::MatExpr is created which holds last cv::MatExpr's evaluation and D*0.1. This is why operations more than 3 is absurd. During evaluation, since the type is not floating-point, the 0.x is lost.

To avoid this problem, just set type to CV_32F or CV_64F.

Your im*0.1f + im*0.1f + im*0.1f + ... is actually working like this.

result = im*0.1f + im*0.1f + im*0.1f + im*0.1f +   ...
  ∧         │         │         │         │
  │      MatExpr   MatExpr   MatExpr   MatExpr     ...
  │         │         │         │         │
  │         └────┬────┘         │         │
  │           MatExpr           │         │
  │eval          │              │         │
  │              └>────────┬───<┘         │
  │              eval   MatExpr           │
  │                        │              │
  │                        └>────────┬───<┘
  │                        eval   MatExpr
  │                                 ...            ...
  │      
  │                                  └>────────┬───<┘
  │                                  eval   MatExpr
  │                                            │
  └<──────────────────────────────────────────<┘

matop.cpp

cv::Mat*double

class MatOp_AddEx : public MatOp
{
public:

    void assign(const MatExpr& expr, Mat& m, int type=-1) const;
    void add(const MatExpr& e1, const Scalar& s, MatExpr& res) const;
    static void makeExpr(MatExpr& res, const Mat& a, const Mat& b, double alpha, double beta, const Scalar& s=Scalar());

    // erased others for clarity
};

static MatOp_AddEx g_MatOp_AddEx; // only visible in matop.cpp

MatExpr operator * (const Mat& a, double s)
{
    MatExpr e;
    MatOp_AddEx::makeExpr(e, a, Mat(), s, 0); // meaning a*s + NULL-mat*0
    return e;
}

inline void MatOp_AddEx::makeExpr(MatExpr& res, const Mat& a, const Mat& b, double alpha, double beta, const Scalar& s)
{
    res = MatExpr(&g_MatOp_AddEx, 0, a, b, Mat(), alpha, beta, s); // MatExpr constructor
}

cv::MatExpr + cv::MatExpr

MatExpr operator + (const MatExpr& e1, const MatExpr& e2)
{
    MatExpr en;
    e1.op->add(e1, e2, en); // MatOp_AddEx inherits MatOp
    return en;
}

void MatOp::add(const MatExpr& e1, const MatExpr& e2, MatExpr& res) const
{
    CV_INSTRUMENT_REGION()

    if( this == e2.op )
    {
        double alpha = 1, beta = 1;
        Scalar s;
        Mat m1, m2;
        if( isAddEx(e1) && (!e1.b.data || e1.beta == 0) )
        {
            m1 = e1.a;
            alpha = e1.alpha;
            s = e1.s;
        }
        else
            e1.op->assign(e1, m1);  // <- Evaluation. Remember that type is set to auto(-1) by default

        if( isAddEx(e2) && (!e2.b.data || e2.beta == 0) )
        {
            m2 = e2.a;
            beta = e2.alpha;
            s += e2.s;
        }
        else
            e2.op->assign(e2, m2);  // <- Evaluation
        MatOp_AddEx::makeExpr(res, m1, m2, alpha, beta, s);
    }
    else
        e2.op->add(e1, e2, res);  // <- Evaluation

}

void MatOp_AddEx::assign(const MatExpr& e, Mat& m, int _type) const
{
    Mat temp, &dst = _type == -1 || e.a.type() == _type ? m : temp;
    if( e.b.data )
    {
        if( e.s == Scalar() || !e.s.isReal() )
        {
            if( e.alpha == 1 )
            {
                if( e.beta == 1 )
                    cv::add(e.a, e.b, dst);
                else if( e.beta == -1 )
                    cv::subtract(e.a, e.b, dst);
                else
                    cv::scaleAdd(e.b, e.beta, e.a, dst);
            }
            else if( e.beta == 1 )
            {
                if( e.alpha == -1 )
                    cv::subtract(e.b, e.a, dst);
                else
                    cv::scaleAdd(e.a, e.alpha, e.b, dst);
            }
            else
                cv::addWeighted(e.a, e.alpha, e.b, e.beta, 0, dst);

            if( !e.s.isReal() )
                cv::add(dst, e.s, dst);
        }
        else
            cv::addWeighted(e.a, e.alpha, e.b, e.beta, e.s[0], dst);
    }
    else if( e.s.isReal() && (dst.data != m.data || fabs(e.alpha) != 1))
    {
        e.a.convertTo(m, _type, e.alpha, e.s[0]);
        return;
    }
    else if( e.alpha == 1 )
        cv::add(e.a, e.s, dst);
    else if( e.alpha == -1 )
        cv::subtract(e.s, e.a, dst);
    else
    {
        e.a.convertTo(dst, e.a.type(), e.alpha);
        cv::add(dst, e.s, dst);
    }

    if( dst.data != m.data )
        dst.convertTo(m, m.type());
}