Pentadiagonal Matrix in OpenCV

Question

I'm trying to implement the active contour models algorithm with c++ and opencv in VisualStudio. I'm having trouble creating a pentadiagonal matrix, that is a matrix with the five main diagonals with values and zeros in the rest of the matrix. I've already implemented this in MatLab, that came out like this:

K = diag(repmat(a,1,n));
K = K + diag(repmat(b,1,n-1),1) + diag(repmat(b,1,n-1),-1)...
    + diag(b,n-1) + diag(b,-n+1);
K = K + diag(repmat(c,1,n-2),2) + diag(repmat(c,1,n-2),-2)...
+ diag([c c],n-2) + diag([c c],-n+2);

How do I do the same thing in opencv? Thanks!

Answer 1

You can wrap the Matlab functions repmat and diag as in the code below.

Then you can write the OpenCV code as in Matlab. The only difference is that you cannot concatenate 2 matrices with [c c], but you should use repmat(c, 2, 1), or repmat(c,1,2) for [c; c].

Code:

#include <opencv2\opencv.hpp>
using namespace cv;

Mat repmat(const Mat& src, int nx, int ny)
{
    return repeat(src, nx, ny);
}

Mat diag(const Mat& src, int k=0)
{
    // src must be a row or column matrix
    //CV_Assert(src.rows == 1 || src.cols == 2);
    int n = src.rows * src.cols;

    Mat1d _src;
    src.convertTo(_src, CV_64F);

    // Create output matrix of correct dimension
    int dim = n + abs(k);
    Mat1d _dst(dim, dim, double(0.0));


    // Select the ranges where to put src data
    Range rrows;
    Range rcols;
    if (k >= 0)
    {
        rrows = Range(0, n);
        rcols = Range(k, k+n);
    }
    else
    {
        rrows = Range(-k, -k + n);
        rcols = Range(0, n);
    }

    // Create square n x n submatrix
    Mat1d sub(_dst(rrows, rcols));

    // Put data on the diagonal of the submatrix
    for (int i = 0; i < n; ++i)
    {
        sub(i, i) = _src(i);
    }

    Mat dst;
    _dst.convertTo(dst, src.type());
    return dst;
}

int main()
{
    Mat a;
    Mat b;
    Mat c;
    int n;

    // ... init a, b, c, n

    Mat K;
    K = diag(repmat(a, 1, n));
    K = K + diag(repmat(b, 1, n - 1), 1) + diag(repmat(b, 1, n - 1), -1) + diag(b, n - 1) + diag(b, -n + 1);
    K = K + diag(repmat(c, 1, n - 2), 2) + diag(repmat(c, 1, n - 2), -2) + diag(repmat(b, 2, 1), n - 2) + diag(repmat(b, 2, 1), -n + 2);

    return 0;
}