快速中值滤波算法之黄氏算法

#include <stdio.h>
#include <iostream>
#include <immintrin.h>
#include <opencv2/opencv.hpp>
using namespace cv;
using namespace std;

//计算中值
int getMediaValue(const int hist[], int thresh) {
  int sum = 0;
  for (int i = 0; i < 256; i++) {
    sum += hist[i];
    if (sum >= thresh) {
      return i;
    }
  }
  return 255;
}
//快速中值滤波，灰度图
Mat fastMedianBlur(Mat src, int diameter) {
  int row = src.rows;
  int col = src.cols;
  Mat dst(row, col, CV_8UC1);
  int Hist[256] = { 0 };
  int radius = (diameter - 1) / 2;
  int windowSize = diameter * diameter;
  int threshold = windowSize / 2 + 1;
  uchar *srcData = src.data;
  uchar *dstData = dst.data;
  int right = col - radius;
  int bot = row - radius;
  for (int j = radius; j < bot; j++) {
    for (int i = radius; i < right; i++) {
      //每一行第一个待滤波元素建立直方图
      if (i == radius) {
        memset(Hist, 0, sizeof(Hist));
        for (int y = j - radius; y <= min(j + radius, row); y++) {
          for (int x = i - radius; x <= min(i + radius, col); x++) {
            uchar val = srcData[y * col + x];
            Hist[val]++;
          }
        }
      }
      else {
        int L = i - radius - 1;
        int R = i + radius;
        for (int y = j - radius; y <= min(j + radius, row); y++) {
          //更新左边一列
          Hist[srcData[y * col + L]]--;
          //更新右边一列
          Hist[srcData[y * col + R]]++;
        }
      }
      uchar medianVal = getMediaValue(Hist, threshold);
      dstData[j * col + i] = medianVal;
    }
  }
  //边界直接赋值
  for (int i = 0; i < col; i++) {
    for (int j = 0; j < radius; j++) {
      int id1 = j * col + i;
      int id2 = (row - j - 1) * col + i;
      dstData[id1] = srcData[id1];
      dstData[id2] = srcData[id2];
    }
  }
  for (int i = radius; i < row - radius - 1; i++) {
    for (int j = 0; j < radius; j++) {
      int id1 = i * col + j;
      int id2 = i * col + col - j - 1;
      dstData[id1] = srcData[id1];
      dstData[id2] = srcData[id2];
    }
  }
  return dst;
}

void meanFilter(int height, int width, unsigned char * __restrict src, unsigned char * __restrict dst) {
  for (int i = 1; i < height - 1; i++) {
    for (int j = 1; j < width - 1; j++) {
      unsigned char a[9];
      a[0] = src[i * width + j];
      a[1] = src[i * width + j + 1];
      a[2] = src[i * width + j - 1];

      a[3] = src[(i + 1) * width + j];
      a[4] = src[(i + 1) * width + j + 1];
      a[5] = src[(i + 1) * width + j - 1];

      a[6] = src[(i - 1) * width + j];
      a[7] = src[(i - 1) * width + j + 1];
      a[8] = src[(i - 1) * width + j - 1];
      for (int ii = 0; ii < 5; ii++) {
        for (int jj = ii + 1; jj < 9; jj++) {
          if (a[ii] > a[jj]) {
            unsigned char temp = a[ii];
            a[ii] = a[jj];
            a[jj] = temp;
          }
        }
      }
      dst[i * width + j] = a[4];
    }
  }
  for (int i = 0; i < width; i++) {
    dst[i] = src[i];
    dst[(height - 1) * width + i] = src[(height - 1) * width + i];
  }
  for (int i = 0; i < height; i++) {
    dst[i * width] = src[i * width];
    dst[i * width + width - 1] = src[i * width + width - 1];
  }
}
Mat speed_MedianFilter(Mat src) {
  int row = src.rows;
  int col = src.cols;
  unsigned char * data = (unsigned char *)src.data;
  unsigned char *dst = new unsigned char[row * col];
  //for (int i = 0; i < row * col; i++) {
  //  printf("%d\n", data[i]);
  //}
  meanFilter(row, col, data, dst);
  Mat res(row, col, CV_8UC1, dst);
  return res;
}

Mat speed_rgb2gray(Mat src) {
  Mat dst(src.rows, src.cols, CV_8UC1);
#pragma omp parallel for num_threads(4)
  for (int i = 0; i < src.rows; i++) {
    for (int j = 0; j < src.cols; j++) {
      dst.at<uchar>(i, j) = ((src.at<Vec3b>(i, j)[0] << 18) + (src.at<Vec3b>(i, j)[0] << 15) + (src.at<Vec3b>(i, j)[0] << 14) +
        (src.at<Vec3b>(i, j)[0] << 11) + (src.at<Vec3b>(i, j)[0] << 7) + (src.at<Vec3b>(i, j)[0] << 7) + (src.at<Vec3b>(i, j)[0] << 5) +
        (src.at<Vec3b>(i, j)[0] << 4) + (src.at<Vec3b>(i, j)[0] << 2) +
        (src.at<Vec3b>(i, j)[1] << 19) + (src.at<Vec3b>(i, j)[1] << 16) + (src.at<Vec3b>(i, j)[1] << 14) + (src.at<Vec3b>(i, j)[1] << 13) +
        (src.at<Vec3b>(i, j)[1] << 10) + (src.at<Vec3b>(i, j)[1] << 8) + (src.at<Vec3b>(i, j)[1] << 4) + (src.at<Vec3b>(i, j)[1] << 3) + (src.at<Vec3b>(i, j)[1] << 1) +
        (src.at<Vec3b>(i, j)[2] << 16) + (src.at<Vec3b>(i, j)[2] << 15) + (src.at<Vec3b>(i, j)[2] << 14) + (src.at<Vec3b>(i, j)[2] << 12) +
        (src.at<Vec3b>(i, j)[2] << 9) + (src.at<Vec3b>(i, j)[2] << 7) + (src.at<Vec3b>(i, j)[2] << 6) + (src.at<Vec3b>(i, j)[2] << 5) + (src.at<Vec3b>(i, j)[2] << 4) + (src.at<Vec3b>(i, j)[2] << 1) >> 20);
    }
  }
  return dst;
}

int main() {
  Mat src = cv::imread("F:\\1.jpg");
  src = speed_rgb2gray(src);
  Mat dst1 = fastMedianBlur(src, 3);
  Mat dst2 = speed_MedianFilter(src);
  cv::imshow("dst1", dst1);
  cv::imshow("dst2", dst2);
  int row = src.rows;
  int col = src.cols;
  for (int i = 0; i < row; i++) {
    for (int j = 0; j < col; j++) {
      if (dst1.at<uchar>(i, j) != dst1.at<uchar>(i, j)) {
        printf("Failture\n");
      }
    }
  }
  cv::waitKey(0);
  return 0;
}