Tiff图像直方图均衡化以及规定化算法

1.图像的截取

tiff格式图片大多为16位数的灰度片，往往色阶比较集中，肉眼难以识别，所以需要图片处理，考虑的方法为，截取集中的像素值段映射为8位图片后再处理，映射不应只考虑线性函数，可以根据具体图像加入γ因子，之后再考虑Tiff图像直方图均衡化以及规定化算法，这里不再赘述(使用的UI和库为opencv+qt)。

以下为部分示例代码：

    char g_GammaLUT[65536]; //全局数组
   
    Mat matTmp(Size(mimg.cols,mimg.rows), CV_8UC1);

    int win_min= m_histogramPara.level - m_histogramPara.wide/2;   //截取段最小位置
    if(win_min < 0) win_min =0;
    int win_max= m_histogramPara.level + m_histogramPara.wide/2;   //截取段最大位置
    if(win_max> 65535) win_max =65535;
    BuildTable(win_min,win_max,m_histogramPara.gama);
    for(int y = 0; y<qimg.height(); y++){                         //给图片的每个像素点用[0-255]范围重新赋值
        uchar * line = (uchar *)qimg.scanLine(y);
        for(int x = 0; x<qimg.width(); x++){
            if(ui->checkBox->isChecked())
                line[x]=255-g_GammaLUT[mimg.at<ushort>(y,x)];
            else
                line[x]=g_GammaLUT[mimg.at<ushort>(y,x)];
            matTmp.at<char>(y,x)=g_GammaLUT[mimg.at<ushort>(y,x)];
        }
    }
  

    scene->oneimg_theRect = QRect(0, 0, qimg.width(), qimg.height());
    scene->setBackgroundBrush(QBrush(QPixmap::fromImage(qimg)));

void BuildTable(int min ,int max ,float fPrecompensation )
{
    int i;
    float f;
    if(max> 65534 )max =65535;
    if(min >= max || min <0 || max >65535 )
    {
        qDebug()<<"error";
    }
    else
    {
        for( i=0;i<min;i++)
        {
            g_GammaLUT[i]=0;
        }
        for( i=max;i<65536;i++)
        {
            g_GammaLUT[i]=0xff;
        }
        for( i=min;i<max;i++)
        {
            int tmpi = i-min;
            int range =max-min;
            f=(tmpi+0.5F)/range;//归一化
            f=(float)pow(f,fPrecompensation); //伽马因子控制亮度，理论上只要值域和定义域都为[0-1]，
            //且为连续单调递增的函数都可以拿来做处理，来改变图像的对比度以及亮度。
            g_GammaLUT[i]=(char)(f*256-0.5F);//反归一化
        }
    }
}

2.直方图的均衡化

(1)数学原理

将图片的色阶图可以看成概率密度函数，由于原图分布不均匀，会导致对比度变差，那么如果让像素值在区域内均匀分布，那么识别度将会极大提高。根据概率论的定理：

设随机变量X具有概率密度函数 $f_{x} (x), - \infty < x < \infty$ , 又设函数 $g (x)$ 处处可导且恒有 $g^{'} (x) > 0$ (或且恒有 $g^{'} (x) < 0$ )，则 $Y = g (X)$ 是连续型随机变量，其概率密度为:
$f_{y} (y) = \binom{f_{x} [h (y)] | h^{'} (y) |, a < y < b}{0, 其他}$
其中 $a = m i n - \infty < x < + \infty$ , $b = m a x - \infty < + \infty$ , $h (y)$ 是 $g (x)$ 的反函数。

这个定理在书上很容易证明！我们现在将像素范围归一化处理，即所有像素位于 $[0 - 1]$ 的范围中，然后已之知我们的像素直方图为 $f_{x} (x)$ 和随机变量 $X$ ， $Y = g (X)$ 是我们想要得到均匀分布的随机变量，那么重点就在于如何求这个映射关系 $g (X)$ ，公式变形可得 : $f_{y} (y) = f_{x} (x) | h^{'} (y) |, a < y < b \to f_{y} (y) / | h^{'} (y) | = f_{x} (x)$

因为 $h (y) 和 g (x)$ 互为反函数那么根据定理：原函数的导数是反函数导数的倒数得知 $g^{'} (x) = 1 / h^{'} (y)$ ,即

$f_{y} (y) | g^{'} (x) | = f_{x} (x)$

对于这个等式对 $d x$ 进行积分可得

$f_{y} (y) g (x) = \int f_{x} (x) d x$ ,

因此如果确定 $f_{y} (y)$ 那么 $g (x)$ 就可以得出等式。由于 $f_{y} (y)$ 是均匀分布，由概率密度公式已知：

$f (x) = 1 / (b - a), a < X < b$ ,

我们已经做了归一化处理范围为 $(0, 1)$ ,所以 $f_{y} (y) = 1$ ,由此得知

$g (x) = \int f_{x} (x) d x$ 。

(2)加权均衡化算法

加权均衡化算法，需要加权矩阵，对每个像素点做加权矩阵范围的均衡化处理，作用是增加局部对比度，原理如下图所示：

(3)效果图及示范代码

   if(m_histogramPara.balanceHistogram){
    equalizeHist(matTmp, matTmp); //全局直方图均衡化，opencv默认的均衡化函数，直接调用
    }else{
        if(m_histogramPara.localBalanceHistogram){
            matTmp=calcLocalBalanceHistogram(matTmp,m_histogramPara.balanceLevel); //此处为加权均衡化函数，即将原始图像后做加权均衡化
        }
    }
    for(int y = 0; y<qimg.height(); y++){  //显示代码，QImage与mat的转换
        uchar * line = (uchar *)qimg.scanLine(y);
        for(int x = 0; x<qimg.width(); x++){
            line[x]=matTmp.at<uchar>(y,x);
            //                 matTmp.at<ushort>(y,x)=matTmp.at<ushort>(y,x);
        }
    }
    
    
    calcLocalBalanceHistogram(Mat &mat, int balanceLevel) //balanceLevel 均衡化块的大小
{

    mutex.lock();
    float valueCount[256];
    float cumulativeDistributionArray[256];//
    if(mat.cols<balanceLevel||mat.rows<balanceLevel){
        qDebug()<<"photo is too small!";
        return Mat();
    }
    qDebug()<<"photo is custom big!";
    int imageRows=mat.rows;
    int imageCols=mat.cols;
    Mat newMat=mat.clone();

    for(int row=balanceLevel/2;row<imageRows-1-balanceLevel/2;row++){

        for(int column=balanceLevel/2;column<imageCols-1-balanceLevel/2;column++){
            if(column==balanceLevel/2){

                Mat tmpMat=mat(Rect(column-balanceLevel/2,row-balanceLevel/2,balanceLevel,balanceLevel));
                  qDebug()<<mat.at<ushort>(column-balanceLevel/2,row-balanceLevel/2);

                calcValueCount(tmpMat,balanceLevel,valueCount);   //算得边角矩阵的像素频率值，这个只能硬算
                    
            }else{

                Mat oldColumn=mat.col(column-balanceLevel/2-1);
                Mat newColumn=mat.col(column+balanceLevel/2);
                for(int i=row-balanceLevel/2;i<row+balanceLevel/2;i++){
                    if(oldColumn.at<uchar>(i,0)!=newColumn.at<uchar>(i,0)){
                        valueCount[oldColumn.at<uchar>(i,0)]-=(float)1/row*row;  //矩阵位移1列后，仅仅是最左列和最右边发生变化
                        valueCount[newColumn.at<uchar>(i,0)]+=(float)1/row*row;  //判断某像素值的频率变化  
                    }
                        ;
                }

            }
            cumulativeDistributionArray[0]=valueCount[0];
            for(int i=1;i<256;i++){
                cumulativeDistributionArray[i]=cumulativeDistributionArray[i-1]+valueCount[i]; //积分运算，求得Y=g(x)的值
            }

            newMat.at<uchar>(row,column)=(float)(255)*cumulativeDistributionArray[mat.at<uchar>(row,column)];

        }
    }
    mutex.unlock();
    return newMat;
}

void calcValueCount(Mat &mat,int row,float *valueCount){
    qDebug()<<mat.rows<<mat.cols;
    for(int i=0;i<256;i++){
        valueCount[i]=0;
    }
    for(int i=0;i<row;i++){
        for(int j=0;j<row;j++){
            valueCount[mat.at<uchar>(i,j)]+=(float)1.0/(row*row); //计算频率出现的次数
        }
    }
}

3.直方图的规定化

(1)数学原理

如果想得到特定概率密度 $f_{z} (z)$ 的直方图，例如提高图像亮度,使概率分布位于偏1的位置： $f_{z} (z) = 2 z$ ,那么可以考虑将规定化的直方图也均衡化，由于都是均衡化为 $f_{y} (y)$ 这样可以直接算出映射关系,因为 $Y = g (x)$ 而且 $Y = k (z)$ ,那么

$z = k^{- 1} (g (x))$ ,

例如 $f_{z} (z) = 2 z$ ,由直方均化可算得:

$Y = k (z) = z^{2}$ ,

求反函数：

$k^{- 1} (y) = \sqrt{y}$ ,带入方程的最终解为

$z = \sqrt{g (x)}$

注意: $f (z)$ 的设定要满足之前的定理限制，且 $\int_{0}^{1}$ 上的积分为1，因为概率最大值为1。

(2)效果图及示范代码


  

    equalizeHist(matTmp, matTmp);
    matTmp=replaceHistogram(matTmp);
    
    replaceHistogram(Mat &mat)
{
    Mat newMat=mat.clone();
    int imageRows=mat.rows;
    int imageCols=mat.cols;
    for(int row=0;row<imageRows;row++){
        for(int column=0;column<imageCols;column++){
    double var=mat.at<uchar>(row,column); //直接拿现有的均衡化值，来计算，这里我取巧了，精度会比较模糊，最好还是自己先计算g(x).
    newMat.at<uchar>(row,column)=sqrt(var/256)*256;//注意先归一化，再运算，再反归一化，
        }
    }
    return newMat;
}


    for(int y = 0; y<qimg.height(); y++){  //显示代码，QImage与mat的转换
        uchar * line = (uchar *)qimg.scanLine(y);
        for(int x = 0; x<qimg.width(); x++){
            line[x]=matTmp.at<uchar>(y,x);
            //                 matTmp.at<ushort>(y,x)=matTmp.at<ushort>(y,x);
        }
    }

Tiff图像直方图均衡化以及规定化算法

1.图像的截取

2.直方图的均衡化

(1)数学原理

(2)加权均衡化算法

(3)效果图及示范代码

3.直方图的规定化

(1)数学原理

(2)效果图及示范代码

文章评论