opencv中熵的求法

Question

我需要一个像matlab中的entropyfilt()这样的函数，它不存在于opencv中。

在matlab中，J = entropyfilt(I)返回数组J，其中每个输出像素包含输入图像I中相应像素周围9乘9邻域的熵值。

我在c++中编写了一个函数来实现它，像素的熵如下所示：

1. 使用带有适当设置的掩码参数的cvCalHist获取图像ROI (即9*9矩形)。
2. 将直方图规范化，使其回收箱之和等于1。
3. 使用(Shannon)熵的公式。

我在下面列出了C++代码：

GetLocalEntroyImage( const IplImage*gray_src,IplImage*entopy_image){
    int hist_size[]={256};
    float gray_range[]={0,255};
    float* ranges[] = { gray_range};
    CvHistogram * hist = cvCreateHist( 1, hist_size, CV_HIST_SPARSE, ranges,1);
    for(int i=0;i<gray_src.width;i++){
            for(int j=0;j<gray_src.height;j++){
                //calculate entropy for pixel(i,j) 
                //1.set roi rect(9*9),handle edge pixel
                CvRect roi;
                int threshold=Max(0,i-4);
                roi.x=threshold;
                threshold=Max(0,j-4);
                roi.y=threshold;
                roi.width=(i-Max(0,i-4))+1+(Min(gray_src->width-1,i+4)-i);
                roi.height=(j-Max(0,j-4))+1+(Min(gray_src->height-1,j+4)-j);
                cvSetImageROI(const_cast<IplImage*>(gray_src),roi);
                IplImage*gray_src_non_const=const_cast<IplImage*>(gray_src);                            

                //2.calHist,here I chose CV_HIST_SPARSE to speed up
                cvCalcHist( &gray_src_non_const, hist, 0, 0 );*/
                cvNormalizeHist(hist,1.0);
                float total=0;
                float entroy=0;

               //3.get entroy
                CvSparseMatIterator it;
                for(CvSparseNode*node=cvInitSparseMatIterator((CvSparseMat*)hist-   >bins,&it);node!=0;node=cvGetNextSparseNode(&it)){
                float gray_frequency=*(float*)CV_NODE_VAL((CvSparseMat*)hist->bins,node);
                entroy=entroy-gray_frequency*(log(gray_frequency)/log(2.0f));//*(log(gray_frequency)/log(2.0))
                }
                ((float*)(local_entroy_image->imageData + j*local_entroy_image->widthStep))[i]=entroy;
                cvReleaseHist(&hist);
            }
        }
        cvResetImageROI(const_cast<IplImage*>(gray_src));
    }

然而，的代码太慢了。我在一张600*1200的图像中测试了它，它的成本为120 s，而matlab中的entroyfilt只需要5s.。

是否有人知道如何加快它的速度，或者知道任何其他好的实现

Answer 1

在您的代码中，最大的慢速是：log(gray_frequency)/log(2.0f))。

您不应该调用cvNormalizeHist()。你知道回收箱要加到81，所以从计算出的熵中减去81 * log(81)/log(2) (当然，这是一个常数，不是每次循环计算的)。如果您不对组图进行规范化，则其条目将是整数，您可以使用它们访问查找表。

由于您有一个9x9内核，所以gray_frequency的最大值是81 (只要您不对直方图进行规范化)，您可以轻松地用一个预先计算的表的一个查找来替换对log()的两个调用。这将产生巨大的变化。您可以初始化这样的表：

    double entropy_table[82]; // 0 .. 81
    const double log2 = log(2.0);
    entropy_table[0] = 0.0;
    for(int i = 1; i < 82; i ++)
    {
        entropy_table[i] = i * log(double(i)) / log2;
    }

后来，它只是：

entroy -= entropy_table[gray_frequency];

此外，您可能会发现，实现您自己的组织代码是一个胜利。例如，如果您有一个小内核，您可以跟踪您将要使用的回收箱，并且只清除这些。但既然你用的是81/256个垃圾箱，这可能是不值得的。

另一个地方，你可以得到一个速度是在边缘像素处理。您正在检查每个像素。但是，如果对边界像素和内部像素有单独的循环，则可以避免对max和min的多次调用。

如果这还不够快，您可以考虑在条纹上使用parallel_for。作为如何这样做的一个很好的例子，请看OpenCV的形态过滤器的源代码。

Answer 2

我检查了entropyfilt的源代码，它位于"entropyfilt.m“中。

它首先垫子src垫，然后调用entropyfiltmex.。

我们知道entropyfiltmex是用C++代码编写的(参考MEX文件file)，并且可以在Matlab文件中找到这些C++源代码文件。

我检查了entroyfiltemex.cpp，主要逻辑是：

void local_entropy(_T *inBuf, double *outBuf){
  ......
  for (p = 0; p < numElements; p++)
        {           
            nhSetWalkerLocation(walker,p);

            // Get Idx into image
            while (nhGetNextInboundsNeighbor(walker, &n, NULL))
            {
                histCountPtr[(int) inBuf[n]]++;
            }

            // Calculate Entropy based on normalized histogram counts
            // (sum should equal one).
            for (k = 0; k < numBins;k++)
            {
                if (histCountPtr[k] != 0)
                {
                    temp = (double) histCountPtr[k] / numNeighbors;

                    // log base 2 (temp) = log(temp) / log(2)
                    entropy = temp * (log(temp)/log((double) 2));
                    outBuf[p] -= entropy;

                    //re-initialize for next neighborhood
                    histCountPtr[k] = 0;
                }
            }
        }
......
}

在这里，nhSetWalkerLocation和nhGetNextInboundsNeighbor是Matlab的邻居操作。

根据Matlab的源代码，非常感谢@B.，我实现了一个改进了这些方面的新版本:

1. 先拍拍图像
2. 避免调用opencv cvCalHist() func，使用hist256获取直方图。
3. 重用matlab邻域操作，使点数学快速。
4. 使用entropy_table保存log()结果，这确实产生了很大的差异(40多个到3s)。

下面是代码：

    void ImageProcess::GetLocalEntroyImage( const IplImage*gray_src,CvRect roi_rect,IplImage*local_entroy_image,IplImage*mask){
        using namespace cv;
        clock_t func_begin,func_end;
        func_begin=clock();
        //1.define nerghbood model,here it's 9*9
        int neighbood_dim=2;
        int neighbood_size[]={9,9};

        //2.Pad gray_src
        Mat gray_src_mat(gray_src);
        Mat pad_mat;
        int left=(neighbood_size[0]-1)/2;
        int right=left;
        int top=(neighbood_size[1]-1)/2;
        int bottom=top;
        copyMakeBorder(gray_src_mat,pad_mat,top,bottom,left,right,BORDER_REPLICATE,0);
        IplImage*pad_src=&IplImage(pad_mat);
        roi_rect=cvRect(roi_rect.x+top,roi_rect.y+left,roi_rect.width,roi_rect.height);

        //3.initial neighbood object,reference to Matlab build-in neighbood object system
        int element_num=roi_rect.width*roi_rect.height;
        //here,implement a histogram by ourself ,each bin calcalate gray value frequence
        int hist_count[256]={0};
        int neighbood_num=1;
        for(int i=0;i<neighbood_dim;i++)
            neighbood_num*=neighbood_size[i];
        //neighbood_corrds_array is a neighbors_num-by-neighbood_dim array containing relative offsets
        int*neighbood_corrds_array=(int*)malloc(sizeof(int)*neighbood_num*neighbood_dim);
        //Contains the cumulative product of the image_size array;used in the sub_to_ind and ind_to_sub calculations.
        int *cumprod;
        cumprod = (int *)malloc(neighbood_dim * sizeof(*cumprod));
        cumprod[0]=1;
        for(int i=1;i<neighbood_dim;i++){
            cumprod[i]=cumprod[i-1]*neighbood_size[i-1];
        }
        int*image_cumprod=(int*)malloc(2*sizeof(*image_cumprod));
        image_cumprod[0]=1;
        image_cumprod[1]=pad_src->width;
        //initialize neighbood_corrds_array
        int p;
        int q;
        int*coords;
        for(p=0;p<neighbood_num;p++){
            coords=neighbood_corrds_array+p*neighbood_dim;
            ind_to_sub(p, neighbood_dim, neighbood_size, cumprod, coords);
            for (q = 0; q < neighbood_dim; q++)
            {
                coords[q] -= (neighbood_size[q] - 1) / 2;
            }
        }
        //initlalize neighbood_offset in use of neighbood_corrds_array
        int*neighbood_offset=(int*)malloc(sizeof(int)*neighbood_num);
        int*elem;
        for(int i=0;i<neighbood_num;i++){
            elem=neighbood_corrds_array+i*neighbood_dim;
            neighbood_offset[i]=sub_to_ind(elem, image_cumprod,2);
        }

        //4.calculate entroy for pixel
        uchar*array=(uchar*)pad_src->imageData;
        //here,use entroy_table to avoid frequency log function which cost losts of time
        float entroy_table[82];
        const float log2=log(2.0f);
        entroy_table[0]=0.0;
        float frequency=0;
        for(int i=1;i<82;i++){
            frequency=(float)i/81;
            entroy_table[i]=frequency*(log(frequency)/log2);
        }
        int neighbood_index;
        int max_index=pad_src->width*pad_src->height;
        float temp;
        float entropy;
        int current_index=0;
        int current_index_in_origin=0;
        for(int y=roi_rect.y;y<roi_rect.height;y++){
            current_index=y*pad_src->width;
            current_index_in_origin=(y-4)*gray_src->width;
            for(int x=roi_rect.x;x<roi_rect.width;x++,current_index++,current_index_in_origin++){
                for(int j=0;j<neighbood_num;j++){
                    int offset=neighbood_offset[j];
                    neighbood_index=current_index+neighbood_offset[j];
                    hist_count[array[neighbood_index]]++;
                }
                //get entroy
                entropy=0;
                for(int k=0;k<256;k++){
                    if(hist_count[k]!=0){
                        int frequency=hist_count[k];
                        entropy -= entroy_table[hist_count[k]];
                        hist_count[k]=0;
                    }
                }
                ((float*)local_entroy_image->imageData)[current_index_in_origin]=entropy;
            }
        }
        func_end=clock();
        double func_time=(double)(func_end-func_begin)/CLOCKS_PER_SEC;
        cout<<"func time"<<func_time<<endl;
    }

新版本的速度要快得多，在相同的图像上只需要3s。

注意：

1. Matlab中的邻接对象确实很花哨。实际上，我们可以更改这个函数接口，以允许不同的内核大小。现在没有时间了，所以这只是一个快速的reuse.aha

参考资料: 1ftp://196.203.130.15/pub/logiciels/matlab2007/toolbox/images/images/private/entropyfiltmex.h 2ftp://196.203.130.15/pub/logiciels/matlab2007/toolbox/images/images/private/neighborhood.cpp

Answer 3

漂亮的(已经投了票)。以下是一些有助于使用它的更改和注意事项。通常，我修复了内存泄漏以及一些将其转换为c++ opencv的内容(尽管还可以做更多改进)。在ios上也很好。

void getLocalEntropyImage(cv::Mat &gray, cv::Rect &roi, cv::Mat &entropy)
{
        using namespace cv;
        clock_t func_begin, func_end;
        func_begin = clock();
        //1.define nerghbood model,here it's 9*9
        int neighbood_dim = 2;
        int neighbood_size[] = {9, 9};

        //2.Pad gray_src
        Mat gray_src_mat(gray);
        Mat pad_mat;
        int left = (neighbood_size[0] - 1) / 2;
        int right = left;
        int top = (neighbood_size[1] - 1) / 2;
        int bottom = top;
        copyMakeBorder(gray_src_mat, pad_mat, top, bottom, left, right, BORDER_REPLICATE, 0);
        Mat *pad_src = &pad_mat;
        roi = cv::Rect(roi.x + top, roi.y + left, roi.width, roi.height);

        //3.initial neighbood object,reference to Matlab build-in neighbood object system
        //        int element_num = roi_rect.area();
        //here,implement a histogram by ourself ,each bin calcalate gray value frequence
        int hist_count[256] = {0};
        int neighbood_num = 1;
        for (int i = 0; i < neighbood_dim; i++)
                neighbood_num *= neighbood_size[i];

        //neighbood_corrds_array is a neighbors_num-by-neighbood_dim array containing relative offsets
        int *neighbood_corrds_array = (int *)malloc(sizeof(int)*neighbood_num * neighbood_dim);
        //Contains the cumulative product of the image_size array;used in the sub_to_ind and ind_to_sub calculations.
        int *cumprod = (int *)malloc(neighbood_dim * sizeof(*cumprod));
        cumprod[0] = 1;
        for (int i = 1; i < neighbood_dim; i++)
                cumprod[i] = cumprod[i - 1] * neighbood_size[i - 1];
        int *image_cumprod=(int*)malloc(2 * sizeof(*image_cumprod));
        image_cumprod[0] = 1;
        image_cumprod[1]= pad_src->cols;
        //initialize neighbood_corrds_array
        int p;
        int q;
        int *coords;
        for (p = 0; p < neighbood_num; p++){
                coords = neighbood_corrds_array+p * neighbood_dim;
                ind_to_sub(p, neighbood_dim, neighbood_size, cumprod, coords);
                for (q = 0; q < neighbood_dim; q++)
                        coords[q] -= (neighbood_size[q] - 1) / 2;
        }
        //initlalize neighbood_offset in use of neighbood_corrds_array
        int *neighbood_offset = (int *)malloc(sizeof(int) * neighbood_num);
        int *elem;
        for (int i = 0; i < neighbood_num; i++){
                elem = neighbood_corrds_array + i * neighbood_dim;
                neighbood_offset[i] = sub_to_ind(elem, image_cumprod, 2);
        }

        //4.calculate entroy for pixel
        uchar *array=(uchar *)pad_src->data;
        //here,use entroy_table to avoid frequency log function which cost losts of time
        float entroy_table[82];
        const float log2 = log(2.0f);
        entroy_table[0] = 0.0;
        float frequency = 0;
        for (int i = 1; i < 82; i++){
                frequency = (float)i / 81;
                entroy_table[i] = frequency * (log(frequency) / log2);
        }
        int neighbood_index;
        //        int max_index=pad_src->cols*pad_src->rows;
        float e;
        int current_index = 0;
        int current_index_in_origin = 0;
        for (int y = roi.y; y < roi.height; y++){
                current_index = y * pad_src->cols;
                current_index_in_origin = (y - 4) * gray.cols;
                for (int x = roi.x; x < roi.width; x++, current_index++, current_index_in_origin++) {
                        for (int j=0;j<neighbood_num;j++) {
                                neighbood_index = current_index+neighbood_offset[j];
                                hist_count[array[neighbood_index]]++;
                        }
                        //get entropy
                        e = 0;
                        for (int k = 0; k < 256; k++){
                                if (hist_count[k] != 0){
                                        //                                        int frequency=hist_count[k];
                                        e -= entroy_table[hist_count[k]];
                                        hist_count[k] = 0;
                                }
                        }
                        ((float *)entropy.data)[current_index_in_origin] = e;
                }
        }
        free(neighbood_offset);
        free(image_cumprod);
        free(cumprod);
        free(neighbood_corrds_array);

        func_end = clock();
        double func_time = (double)(func_end - func_begin) / CLOCKS_PER_SEC;
        std::cout << "func time" << func_time << std::endl;
}

这里还有一些遗漏的功能。

static int32_t sub_to_ind(int32_t *coords, int32_t *cumprod, int32_t num_dims)
{
        int index = 0;
        int k;

        assert(coords != NULL);
        assert(cumprod != NULL);
        assert(num_dims > 0);

        for (k = 0; k < num_dims; k++)
        {
                index += coords[k] * cumprod[k];
        }

        return index;
}

static void ind_to_sub(int p, int num_dims, const int size[],
                       int *cumprod, int *coords)
{
        int j;

        assert(num_dims > 0);
        assert(coords != NULL);
        assert(cumprod != NULL);

        for (j = num_dims-1; j >= 0; j--)
        {
                coords[j] = p / cumprod[j];
                p = p % cumprod[j];
        }
}

最后，这里是如何使用它，以了解它的外观(示例)。

            cv::Rect roi(0, 0, gray.cols, gray.rows);
            cv::Mat dst(gray.rows, gray.cols, CV_32F);
            getLocalEntropyImage(gray, roi, dst);
            cv::normalize(dst, dst, 0, 255, cv::NORM_MINMAX);
            cv::Mat entropy;
            dst.convertTo(entropy, CV_8U);

这里@熵是你要展示的形象。

​

​