10

I am implementing a Kuwahara filter in C++, with OpenCV to help opening and displaying images. The idea is quite straight forward but somehow I got weird result from it. Here' the cose:

#include "opencv2/opencv.hpp"
#include <iostream>
#include <iomanip>
#include <cmath>

using namespace std;
using namespace cv;

//This class is essentially a struct of 4 Kuwahara regions surrounding a pixel, along with each one's mean, sum and variance.
class Regions{
    int* Area[4];
    int Size[4];
    unsigned long long Sum[4];
    double Var[4];
    int kernel;
public:
    Regions(int _kernel) : kernel(_kernel) {
        for (int i = 0; i<4; i++) {
            Area[i] = new int[kernel*kernel];
            Size[i] = 0;
            Sum[i] = 0;
            Var[i] = 0.0;
        }
    }

    //Update data, increase the size of the area, update the sum
    void sendData(int area, int data){
        Area[area][Size[area]] = data;
        Sum[area] += data;
        Size[area]++;
    }
    //Calculate the variance of each area
    double var(int area) {
        int __mean = Sum[area]/Size[area];
        double temp = 0;
        for (int i = 0; i<Size[area]; i++) {
            temp+= (Area[area][i] - __mean) * (Area[area][i] - __mean);
        }
        if (Size[area]==1) return 1.7e38; //If there is only one pixel inside the region then return the maximum of double
                                           //So that with this big number, the region will never be considered in the below minVar()
        return sqrt(temp/(Size[area]-1));
    }
    //Call the above function to calc the variances of all 4 areas
    void calcVar() {
        for (int i = 0; i<4; i++) {
            Var[i] = var(i);
        }
    }
    //Find out which regions has the least variance
    int minVar() {
        calcVar();
        int i = 0;
        double __var = Var[0];
        if (__var > Var[1]) {__var = Var[1]; i = 1;}
        if (__var > Var[2]) {__var = Var[2]; i = 2;}
        if (__var > Var[3]) {__var = Var[3]; i = 3;}
        return i;
    }

    //Return the mean of that regions
    uchar result(){
        int i = minVar();
        return saturate_cast<uchar> ((double) (Sum[i] *1.0 / Size[i]));
    }
};

class Kuwahara{
private:
    int wid, hei, pad, kernel;
    Mat image;
public:
    Regions getRegions(int x, int y){
        Regions regions(kernel);

        uchar *data = image.data;

        //Update data for each region, pixels that are outside the image's boundary will be ignored.

        //Area 1 (upper left)
        for (int j = (y-pad >=0)? y-pad : 0; j>= 0 && j<=y && j<hei; j++)
            for (int i = ((x-pad >=0) ? x-pad : 0); i>= 0 && i<=x && i<wid; i++) {
                regions.sendData(1,data[(j*wid)+i]);
            }
        //Area 2 (upper right)
        for (int j = (y-pad >=0)? y-pad : 0; j<=y && j<hei; j++)
            for (int i = x; i<=x+pad && i<wid; i++) {
                regions.sendData(2,data[(j*wid)+i]);
            }
        //Area 3 (bottom left)
        for (int j = y; j<=y+pad && j<hei; j++)
            for (int i = ((x-pad >=0) ? x-pad : 0); i<=x && i<wid; i++) {
                regions.sendData(3,data[(j*wid)+i]);
            }
        //Area 0 (bottom right)
        for (int j = y; j<=y+pad && j<hei; j++)
            for (int i = x; i<=x+pad && i<wid; i++) {
                regions.sendData(0,data[(j*wid)+i]);
            }
        return regions;
    }

    //Constructor
    Kuwahara(const Mat& _image, int _kernel) : kernel(_kernel) {
        image = _image.clone();
        wid = image.cols; hei = image.rows;
        pad = kernel-1;
    }

    //Create new image and replace its pixels by the results of Kuwahara filter on the original pixels
    Mat apply(){
        Mat temp;
        temp.create(image.size(), CV_8U);
        uchar* data = temp.data;

        for (int j= 0; j<hei; j++) {
            for (int i = 0; i<wid; i++)
                data[j*wid+i] = getRegions(i,j).result();
        }
        return temp;
    }
};

int main() {
    Mat img = imread("limes.tif", 1);
    Mat gray, dest;
    int kernel = 15;
    gray.create(img.size(), CV_8U);
    cvtColor(img, gray, CV_BGR2GRAY);

    Kuwahara filter(gray, kernel);

    dest = filter.apply();

    imshow("Result", dest);
    imwrite("result.jpg", dest);
    waitKey();
}

And here's the result: enter image description here

As you can see it's different from the correct result, the borders of those limes seem to be duplicated and moved upward. If I apply a 15x15 filter, it gives me a complete mess like this:

enter image description here

I've spent my whole day to debug, but so far nothing is found. I even did the calculation on small images by hand and compare with the result and see no differences. Could anyone help me find out what did I do wrong? Many many thanks.

Max
  • 3,824
  • 8
  • 41
  • 62
  • 1
    There's a lot of code to check there, but one immediate thing I noticed is in your `var()` function you're using integer arithmetic in places and this may lead to some truncation of values. I'm not sure if it will be significant with the range expected but it's worth doing everything at double precision to see if that makes any difference. – Roger Rowland Mar 25 '13 at 06:04
  • @roger_rowland Thank you for taking time reading my code and pointing out that mistake. However, the result doesn't change when I use `double` instead of `int` and forcing all the arithmetic calculations to deal with `double`. – Max Mar 25 '13 at 06:22
  • 3
    This won't solve you problem, but it would be good pratice to use `std::numeric_limits::max()` instead of manually writing the maximum value of `double`. Also, you should use constants to name your areas, say `UPPER_LEFT` instead of `1`. As I said, those details won't solve your problem, but people will be more likely to help if your code is easy to read and self-documented :) – Morwenn Mar 25 '13 at 15:06
  • @Morwenn Cool, thank you. I tried searching for that `numeric_limits` function before I had to hard code the max xalue. – Max Mar 25 '13 at 18:38
  • @crazyfffan: Which implementation did you use to create the groundtruth (the correct result)? – DCS Mar 27 '13 at 07:16
  • @DCS The correct one was given by my teacher. I also tried with lots of examples of kuwahara filter online, and compare with their outcome. My implementation of Kuwahara filter never produces anything close to those example. :(( – Max Mar 27 '13 at 07:31
  • If you need some base code to build from, I found the initial implementation of the Kuwahara filter by Jan Eric Kyprianidis to be a good starting point: http://code.google.com/p/gpuakf/source/browse/glsl/kuwahara.glsl (from this project: http://code.google.com/p/gpuakf/ and this publication: http://www.kyprianidis.com/p/pg2009/ ). It's GL shader code, but is C-like enough that you should be able to adapt it to your situation without too much effort. I used something based on this in a shader here: http://stackoverflow.com/a/9402041/19679 – Brad Larson Mar 27 '13 at 21:42

1 Answers1

7

It turns out that there's nothing wrong with my code, but the way I defined a kernel was the source of problem. My kernel is actually one of four small kuwahara sections, while the correct definition of a kernel is the whole area where data is calculated for each pixel, therefore the area that contains all four sections is actually the kernel. So when talked about a 7x7 "kernel", I actually applied a 15x15 one, and the horrible result came not from a 15x15 kernel as I thought, but from a 31x31. At that size, Kuwahara filter simply doesn't make sense and bizarre results are inevitable.

Max
  • 3,824
  • 8
  • 41
  • 62
  • Old post, but I couldn't reach you directly. I'm researching some edge-preserving filters, and thought that I'd try Kuwahara. Did you get anywhere with this? Any chance that you added any variants to the OpenCV contrib library? – StringTheory Mar 27 '21 at 01:51