一维熵局部熵C++实现

基于OpenCV的一维熵与局部熵C++实现

由于要将MATLAB代码都转换为C++代码，因此开始了艰苦卓绝的码代码过程。这其中又遇到了很多的坑，以及爬坑过程。

我的环境的是Visual studio2017+ OpenCV3

本文内容：

• 一维熵与局部熵的原理
• 详解OpenCV存在的坑
• 局部熵entropyfilt()源代码实现

一维熵与局部熵的原理

熵这个东西，最先是由一个叫香农的人提出来的。所以你在谷歌上搜一维熵或者局部熵，也会出来shannon entropy 相关信息（比如局部熵 ，你搜local shannon entropy也对）。熵的相关计算在MATLAB上都有比较完备的计算函数，比如一维熵的函数是entropy（图片计算出来的结果是一个浮点型数值），局部熵的函数是entropyfilt（计算出来的结果是一张和原图一样大小的熵图）.因此将原来的MATLAB转换为C++时，一个比较方便的是，你可以各种对比计算得到的数值。

1. 一维熵

   图像熵表示的是一种特征的统计形式，它反映了图像中的平均信息量。图像一维熵表示图像中灰度分布的聚集特征所包含的信息量，令 $P_i$ 表示图像中灰度值为$i$ 的像素值所占的比例，则定义灰度图像的一维熵为：

   $P_{i,j}=f(i,j)/N^2$

   其中$P_i$是某个灰度在图片中出现的概率，可以由灰度直方图获得。

2. 局部熵

   局部熵反应的是图片中某一区域所含的信息的多少。对于图像中任意一点$(i,j)$ ，选取局部$M×N$ 窗口（一般是3×3）,定义局部熵为：

   $F_{i,j}=-\sum^M_{i=1}\sum_{j=1}^Np_{i,j}log_2p_{i,j}$

   其中$P_{i,j]}$是$（i.j）$ 处的概率。（这段内容出自《基于局部熵的量子衍生医学超声图像去斑 》）具体原理大家查询这篇论文。

详解OpenCV的坑

在说坑之前我先贴出局部熵的MATLAB的实现代码。

clear all;close all;clc;
img = imread('271.jpg');
[m1 n1]=size(img);
%% 将整个图片外围扩大一圈
x1=zeros(3,424);
img=[x1;img];
img=[img;x1];
y1=zeros(93,3);
img=[y1,img];
img=[img,y1];
w=3;%模板半径
imgn=zeros(m1+6,n1+6);%m*n维矩阵
[m,n]=size(img);
for i=1+w:m-w
    for j=1+w:n-w

        Hist=zeros(1,256);%1*256维
        for p=i-w:i+w%核大小
            for q=j-w:j+w
                Hist(img(p,q)+1)=Hist(img(p,q)+1)+1;%统计局部直方图
            end
        end

        Hist=Hist/sum(Hist);

        for k=1:256
            if Hist(k)~=0
               imgn(i,j)=imgn(i,j)+Hist(k)*log(1/Hist(k));%局部熵
            end
        end

    end
end

imshow(imgn,[])

imgn=entropyfilt(img); %系统的局部熵函数entropyfilt

figure;

imshow(imgn,[]);

这段代码网上有很多博客都贴出来了，实现的效果也很不错。其实代码的逻辑也比较简单。就是先算出来局部的直方图，然后根据熵的公式计算局部熵。这段代码是用在灰度图上的，可以直接运行，没任何问题。

好，我当时看见这段代码，我就乐了。直接按照这个代码写一段自己的C++代码应该也不是很难。然后终于开始了掉坑之旅了。有几个比较坑的地方：

• 图片的格式选择
• C++矩阵与MATLAB矩阵存储的区别
• 插件的坑
• 0的坑

---

1. 图片的格式选择

   由于MATLAB上对数字的类型没有很多的限制，你想用哪种类型用哪种。但是到了OpenCV上，就需要格外注意了。建议大家看这篇《阿洲的程式教學 》以及这篇《OpenCV 之 Mat 类》这两篇博客都比较详细，系统地讲解了OpenCV Mat类型对应的数据类型：

   CV_8U  - 8-bit 无符号整数 ( 0..255 )
   CV_8S  - 8-bit 有符号整数 ( -128..127 )
   CV_16U - 16-bit 无符号整数 ( 0..65535 )
   CV_16S - 16-bit 有符号整数 ( -32768..32767 )
   CV_32S - 32-bit 有符号整数 ( -2147483648..2147483647 )
   CV_32F - 32-bit 浮点数 ( -FLT_MAX..FLT_MAX, INF, NAN )
   CV_64F - 64-bit 浮点数 ( -DBL_MAX..DBL_MAX, INF, NAN )

这很重要，比如你想建立一个Mat Hist 来存储直方图的相关信息，你就需要用CV_8U来存储。如果你需要将最终的局部熵特征值存储到一个图片上，你就需要用到CV_32F。用来存储相关的浮点值。有些博客上说读取灰度图需要采用image.at<uchar>(i,j)其实这个是有些问题的。如果你直接这么读取了一个像素是输出不了一个值的。你需要在前面加上一个int()强制转换才行。因此在整个编码的过程中，被这些数值类型整的精疲力尽。但是也不能全部吐槽，讲点小收获吧：

一个就是像素的读取方式前面已经说了用at。还有一个就是OpenCV的矩阵存储形式Mat 基本就是一个万能的数据结构。用在矩阵，用在向量上，都很好。还有一个就是图像的剪切：Rect rect(top_point,left_point,length,height) img=img(rect)。当然还有好多小的技巧，都可以在源代码里面读出来。

1. C++矩阵与MATLAB矩阵存储的区别

   我相信每一个做过转换的人都会有这个槽点可以吐槽。就是C++矩阵是从0开始的，然而MATLAB是从1开始的。

2. 插件的坑

   其实也不算特别大的坑，只能说是一个比较不合理的地方。如果用OpenCV+vs，一般我推荐大家使用一下Image Watch。这个是OpenCV在vs推出来的图片查看插件，可以看到具体图片的相关信息，具体的下载和使用，可以看这篇博客《VS2015中OpenCV编程插件Image Watch安装和使用介绍 》，我说它坑主要是因为在不调试的情况下，这个插件是没有任何作用的。

3. 0的坑

   这个坑呢，就是个人编程习惯的问题，在C++中由于受精度的影响，C++中最好不要用！=表示不相等，最好使用差值<0.00001表示相等。这个在大量存在浮点运算的程序里显得尤为重要。

局部熵entropyfilt()源代码实现

这段代码最终实现效果还是差强人意，不算太好。但是比网上的大部分都要好，网上的C++局部熵代码，我找了好久每一个靠谱的。我可以给大家推荐几个还算能看的。一个是这个Entropy Filter in OpenCV similar to entropyfilt() function in matlab 在stackoverflow上，这个老外依照另一个回答的代码，写出了自己的代码，但是实现出的效果不好是一块白板。还有一个是GitHub上的代码HuangSuqi/GaoDSH 他写的这段代码主题思想和我之前给出的MATLAB代码差不多，但是我没试。

最终我打算还是用自己的。

#include <iostream>
#include <stdio.h>
#include"opencv2/opencv.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <math.h>
using namespace std;
using namespace cv;
Mat entropyfilt(Mat img,Size size) {
    //模板大小 半径是3  模板大小9*9
    int w = 3;
    int m1 = size.height;
    int n1 = size.width;
    Mat img_padding;
    //先对图片的四周进行补零  
    Mat Hist1 = cv::Mat::zeros(1, 256, CV_32F);
    copyMakeBorder(img, img_padding, 3, 3, 3, 3, BORDER_CONSTANT, cv::Scalar(0, 0, 0,0));


    Mat imgn = Mat::zeros(m1+6, n1+6,CV_64F);

    int m = m1 + 6;
    int n = n1 + 6;

    for (int i = w; i < m - w; i++) {
        for (int j = w; j < n - w;j++) {
            Mat Hist = Mat::zeros(1, 800,CV_32F);
            int m_loop = i;
            int n_loop = j;
            for (int p = m_loop - w; p <= m_loop + w;p++) {
                for (int q = n_loop - w; q <= n_loop + w;q++) {
                    if(int(img_padding.at<uchar>(p, q))<= 800){
                        //cout << "  " << int(img_padding.at<float>(p, q));
                        Hist.at<float>(int(img_padding.at<uchar>(p, q))) =int( Hist.at<float>(int(img_padding.at<uchar>(p, q))))+1;  //这个地方有问题
                        //cout << int(Hist.at<float>(int(img_padding.at<float>(p, q)))) << " ";
                    }
                    else {
                        cout << p << " " << q<<" ";
                        cout << int(img_padding.at<uchar>(p, q));
                        cout << "-----------------------------------------";
                    }
                }

            }
            Hist = Hist / float(sum(Hist));    // 将直方图归一化 
            //cout << "???:" << Hist.at<float>(0, 0);
            //normalize(Hist, Hist,1.0000, 0.0000, NORM_MINMAX);

            /*for (int i = 0; i < 256; i++) {
                if (Hist.at<float>(0,i) != 0) {
                    cout << Hist.at<float>(0, i) << " ";
                }
            }*/
            //cout << double(imgn.at<float>(i,j)) << endl;
            //cout << "什么数值：" << Hist.at<float>(20) << " ";
            for (int k = 0; k <256;k++) {
                if (float(Hist.at<float>(k)) > 0.0001 ) {

                    imgn.at<double>(i, j) = float(imgn.at<double>(i, j)) + float(Hist.at<float>(k))* (float(log2(float(1.0 / Hist.at<float>(k))))); // 局部熵

                }

            }
        }
    }
    Rect rect(3, 3, 424, 87);
    imgn = imgn(rect);
    return imgn;
}
int main(int ,char** argv) {
    Mat src, hist;
    //load image
    src = imread("271.jpg",cv::IMREAD_GRAYSCALE);            //输入图片的方式  
    if (src.empty())
    {
        return -1;
    }
    namedWindow("Image", WINDOW_AUTOSIZE);
    imshow("Image", src);
    Mat img_temp = entropyfilt(src, src.size());
    //进行反色处理
    normalize(img_temp, img_temp, 1.0000, 0.0000, NORM_MINMAX);
    namedWindow("img_temp", WINDOW_AUTOSIZE);
    imshow("img_temp", img_temp);
    waitKey(0);
    return 0;

}