Pixel-wise median of sequence of cv::Mat's

Question

Note: I am NOT asking about Median Filter.

I have a sequence of images let us say:

std::array<cv::Mat,N> sequence;

I want to blend all those images in one. This one image should satisfies:

Each pixel of the new image is the median of its corresponding pixels from the sequence. In other words:

Result(i,j)=median(sequence[0](i,j), sequence[1](i,j), ..., sequence[N](i,j));

Is there built-in function for doing that? What would be the fastest way?

What I tried till now: To iterate over each pixel from all the sequence and sort then take the median then store it in the result. However, it is so overkill.

Answer 1

You can use the following technique if the number of images in your sequence is odd.

• Prepare an N (which is odd) channel image from the sequence of images that you have
• Reshape this image into a 1-channel column vector
• Now apply a median filter of size N to this column vector. As the image is a 1-channel column vector, the filter will calculate the median of the channels(of course, there will be some additional calculations that are not useful to us)
• Reshape this filtered image into its original form having N channels and the original number of rows and columns
• Pick the middle channel of this N-channel image, which contains the median image of the sequence

Below I've illustrated the above items with a simple code and its output.

individual channels

channel 0:
[1, 1, 1;
  1, 1, 1;
  1, 1, 1]
channel 1:
[2, 2, 2;
  2, 2, 2;
  2, 2, 2]
channel 2:
[3, 3, 3;
  3, 3, 3;
  3, 3, 3]
channel 3:
[4, 4, 4;
  4, 4, 4;
  4, 4, 4]
channel 4:
[5, 5, 5;
  5, 5, 5;
  5, 5, 5]

output for N = 3

3-channel image data:
[1, 2, 3, 1, 2, 3, 1, 2, 3;
  1, 2, 3, 1, 2, 3, 1, 2, 3;
  1, 2, 3, 1, 2, 3, 1, 2, 3]
1-channel column vector image data:
[1; 2; 3; 1; 2; 3; 1; 2; 3; 1; 2; 3; 1; 2; 3; 1; 2; 3; 1; 2; 3; 1; 2; 3; 1; 2; 3]
median of the 1-channel column vector image data:
[1; 2; 2; 2; 2; 2; 2; 2; 2; 2; 2; 2; 2; 2; 2; 2; 2; 2; 2; 2; 2; 2; 2; 2; 2; 2; 3]
reshaped filtered image data:
[1, 2, 2, 2, 2, 2, 2, 2, 2;
  2, 2, 2, 2, 2, 2, 2, 2, 2;
  2, 2, 2, 2, 2, 2, 2, 2, 3]
median image data:
[2, 2, 2;
  2, 2, 2;
  2, 2, 2]

output for N = 5

5-channel image data:
[1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5;
  1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5;
  1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5]
1-channel column vector image data:
[1; 2; 3; 4; 5; 1; 2; 3; 4; 5; 1; 2; 3; 4; 5; 1; 2; 3; 4; 5; 1; 2; 3; 4; 5; 1; 2
; 3; 4; 5; 1; 2; 3; 4; 5; 1; 2; 3; 4; 5; 1; 2; 3; 4; 5]
median of the 1-channel column vector image data:
[1; 2; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3
; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 3; 4; 5]
reshaped filtered image data:
[1, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3;
  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3;
  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 5]
median image data:
[3, 3, 3;
  3, 3, 3;
  3, 3, 3]

Code:

// number of channels (= number of images in the sequence)
// N MUST BE ODD
const int N = 5;    
// channel data
uchar ch0[] = {1, 1, 1, 1, 1, 1, 1, 1, 1};
uchar ch1[] = {2, 2, 2, 2, 2, 2, 2, 2, 2};
uchar ch2[] = {3, 3, 3, 3, 3, 3, 3, 3, 3};
uchar ch3[] = {4, 4, 4, 4, 4, 4, 4, 4, 4};
uchar ch4[] = {5, 5, 5, 5, 5, 5, 5, 5, 5};
// images
Mat m0 = Mat(3, 3, CV_8U, ch0);
Mat m1 = Mat(3, 3, CV_8U, ch1);
Mat m2 = Mat(3, 3, CV_8U, ch2);
Mat m3 = Mat(3, 3, CV_8U, ch3);
Mat m4 = Mat(3, 3, CV_8U, ch4);
// prepare image sequence
Mat channels[] = {m0, m1, m2, m3, m4};
// put the images into channels of matrix m
Mat m;
merge(channels, N, m);
// reshape data so that we have a single channel column vector as the image
Mat n = m.reshape(1, m.rows * m.cols * m.channels());
// apply median filter to the 1-channel column vector image. filter size must be the number of channels
Mat med;
medianBlur(n, med, N);

cout << N << "-channel image data:" << endl;
cout << m << endl;

cout << "1-channel column vector image data:" << endl;
cout << n << endl;

cout << "median of the 1-channel column vector image data:" << endl;
cout << med << endl;

// reshape the filtered 1-channel column vector image into its original form having N channels
med = med.reshape(N, m.rows);

cout << "reshaped filtered image data:" << endl;
cout << med << endl;

// split the image
split(med, channels);

// extract the middle channel which contains the median image of the sequence
cout << "median image data:" << endl;
cout << channels[(N+1)/2 - 1] << endl;

Answer 2

You can compute the sequential median for each position using histograms.

Assuming that you're using Mat1b images, each histogram would have 256 values. You need to store the histogram, as well as the sum of all bins:

struct Hist {
    vector<short> h;
    int count;
    Hist() : h(256, 0), count(0) {};
};

The median value is the index in the histogram that corresponds to half of the pixels count / 2. Snippet from Rosetta Code:

int i;
int n = hist.count / 2; // 'hist' is the Hist struct at a given pixel location
for (i = 0; i < 256 && ((n -= hist.h[i]) >= 0); ++i);
// 'i' is the median value

When you add or remove an image, you update the histogram for each pixel location, and recompute the median value. This operation is quite fast because you don't need to sort.

There are some drawback to this:

1. This will work only for uchar values, otherwise the length of each histogram will be too large
2. This approach will use a lot of RAM since it needs rows x cols histograms. It may not work for large images.
3. The following implementation works for single channel images, but it can be easily extended to 3 channels.

4. You can use an approach based on two heaps, or approximate methods. You can find details here:

   • Find running median from a stream of integers
   • Rolling median algorithm in C
   • Find median value from a growing set

Code:

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

using namespace std;
using namespace cv;

struct Hist {
    vector<short> h;
    int count;
    Hist() : h(256, 0), count(0) {};
};

void addImage(vector<Mat1b>& images, Mat1b& img, vector<vector<Hist>>& M, Mat1b& med)
{
    assert(img.rows == med.rows);
    assert(img.cols == med.cols);

    for (int r = 0; r < img.rows; ++r) {
        for (int c = 0; c < img.cols; ++c){

            // Add pixel to histogram
            Hist& hist = M[r][c];
            ++hist.h[img(r, c)];
            ++hist.count;

            // Compute median
            int i;
            int n = hist.count / 2;
            for (i = 0; i < 256 && ((n -= hist.h[i]) >= 0); ++i);

            // 'i' is the median value
            med(r,c) = uchar(i);
        }
    }

    // Add image to my list
    images.push_back(img.clone());
}

void remImage(vector<Mat1b>& images, int idx, vector<vector<Hist>>& M, Mat1b& med)
{
    assert(idx >= 0 && idx < images.size());

    Mat1b& img = images[idx];
    for (int r = 0; r < img.rows; ++r) {
        for (int c = 0; c < img.cols; ++c){

            // Remove pixel from histogram
            Hist& hist = M[r][c];
            --hist.h[img(r, c)];
            --hist.count;

            // Compute median
            int i;
            int n = hist.count / 2;
            for (i = 0; i < 256 && ((n -= hist.h[i]) >= 0); ++i);

            // 'i' is the median value
            med(r, c) = uchar(i);
        }
    }

    // Remove image from list
    images.erase(images.begin() + idx);
}

void init(vector<vector<Hist>>& M, Mat1b& med, int rows, int cols)
{
    med = Mat1b(rows, cols, uchar(0));
    M.resize(rows);
    for (int i = 0; i < rows; ++i) {
        M[i].resize(cols);
    }
}


int main()
{
    // Your images... be sure that they have the same size
    Mat1b img0 = imread("path_to_image", IMREAD_GRAYSCALE);
    Mat1b img1 = imread("path_to_image", IMREAD_GRAYSCALE);
    Mat1b img2 = imread("path_to_image", IMREAD_GRAYSCALE);
    resize(img0, img0, Size(800, 600));
    resize(img1, img1, Size(800, 600));
    resize(img2, img2, Size(800, 600));



    int rows = img0.rows;
    int cols = img0.cols;

    vector<Mat1b> images;   // All your images, needed only if you need to remove an image
    vector<vector<Hist>> M; // histograms
    Mat1b med;              // median image

    // Init data strutctures
    init(M, med, rows, cols);

    // Add images. 'med' will be the median image and will be updated each time
    addImage(images, img0, M, med);
    addImage(images, img1, M, med);
    addImage(images, img2, M, med);

    // You can also remove an image from the median computation
    remImage(images, 2, M, med);

    // Hey, same median as img0 and img1 ;D

    return 0;
}