3.8 直方图规定化（匹配）

直方图均衡化算法可以自动确定灰度变换函数，从而获得具有均匀直方图的输出图像。它主要用于增强动态范围偏小的图像对比度，丰富图像的灰度级。这种方法的优点是操作简单，且结果可以预知，当图像需要自动增强时是一种不错的选择。

但读者有时希望可以对变换过程加以控制，如能够人为地修正直方图的形状，或者说是获得具有指定直方图的输出图像，这样就可以有选择地增强某个灰度范围内的对比度或使图像灰度值满足某种特定的分布。这种用于产生具有特定直方图的图像的方法叫作直方图规定化，或直方图匹配。

3.8.1 理论基础

直方图规定化是在运用均衡化原理的基础上，通过建立原始图像和期望图像（待匹配直方图的图像）之间的关系，使原始图像的直方图匹配特定的形状，从而弥补了直方图均衡化不具备交互作用的特性。

其匹配原理是先对原始的图像均衡化，转换公式如下。

同时对待匹配直方图的图像进行均衡化处理，公式如下。

由于都是均衡化，故可令s=v，则有如下关系。

于是可以按照如下的步骤由输入图像得到一个具有规定概率密度函数的图像。

（1）根据式（3-14）得到变换关系f(r)。

（2）根据式（3-15）得到变换关系g(z)。

（3）求得反变换函数g -1(s)。

（4）对输入图像所有像素应用式（3-16）中的变换，从而得到输出图像。

当然，在实际计算中利用的是上述公式的离散形式，这样就不必去关心函数f(r)、g(z)以及反变换函数g -1(s)具体的解析形式，而可以直接将它们作为映射表处理了。其中，f(r)为输入图像均衡化的离散灰度级映射关系，g(z)为标准图像均衡化的离散灰度级映射关系，而g -1(s)则是标准图像均衡化的逆映射关系，它给出了从经过均衡化处理的标准化图像到原标准图像的离散灰度映射，相当于均衡化处理的逆过程。

3.8.2 MATLAB编程实现

Histeq()函数不仅可以用于直方图均衡化，也可以用于直方图规定化，此时需要提供可选参数hgram。调用语法如下。

      [J, T] = histeq(I, hgram)

函数会将原始图像I处理成一幅以用户指定向量hgram作为直方图的图像。

参数hgram的分量数目即为直方图的收集箱数目。对于double型图像，hgram的元素取值范围是[0, 1]；对于uint8型图像取值范围为[0, 255]；对于uint16型图像取值范围则为[0, 65535]。

其他参数的意义与在直方图均衡化中的相同。

【例3.2】直方图的匹配。

下面的程序实现了从图像I分别到图像I1和I2的直方图匹配。

      I = imread('pout.tif'); %读入原图像
      I1 = imread('coins.png'); %读入要匹配直方图的图像
      I2 = imread('circuit.tif'); %读入要匹配直方图的图像

      % 计算直方图
      [hgram1, x] = imhist(I1);
      [hgram2, x] = imhist(I2);

      % 执行直方图均衡化
      J1=histeq(I, hgram1);
      J2=histeq(I, hgram2);

      % 绘图
      subplot(2,3,1);
      imshow(I); title(’原图’);
      subplot(2,3,2);
      imshow(I1); title(’标准图1');
      subplot(2,3,3);
      imshow(I2); title(’标准图2');
      subplot(2,3,5);
      imshow(J1); title(’规定化到1')
      subplot(2,3,6);
      imshow(J2); title(’规定化到2');

      % 绘直方图
      figure;

      subplot(2,3,1);
      imhist(I); title(’原图’);

      subplot(2,3,2);
      imhist(I1); title(’标准图1');

      subplot(2,3,3);
      imhist(I2); title(’标准图2');

      subplot(2,3,5);
      imhist(J1); title(’规定化到1')

      subplot(2,3,6);
      imhist(J2); title(’规定化到2');

上述程序的运行结果如图3.19和图3.20所示。读者可以看到，经过规定化处理，原图像的直方图与目标图像的直方图变得较为相似。

直方图规定化本质上是一种拟合过程，因此变换得到的直方图与标准目标图像的直方图并不会完全一致。然而即使只是相似的拟合，仍然使规定化的图像在亮度与对比度上具有类似标准图像的特性，这正是直方图规定化的目的所在。

3.8.3 Visual C++实现

根据前文提到的变换方法，可以首先计算原图像直方图均衡化的离散变换关系sk=f（rk）（这里rk=0,1,2,…,255），而后计算标准图像直方图均衡化的离散变换关系vk=g（zk）（这里zk=0,1, 2,…,255），从而进一步得到逆映射关系g-1(sk)或g-1(vk)，再利用关系vk=sk，可得zk=g-1(vk)=g-1(sk)，即对原图像的均衡化结果sk执行标准图像均衡化的反变换g-1(sk)或g-1(vk)。

下面Histst()的实现中首先记录下标准图像或其直方图进行直方图均衡化时每一灰度级别在均衡化后对应的灰度级，为了得到逆变换关系g-1，通过pdTran数组直接建立从均衡化后的灰度级vk到均衡化之前标准图像灰度级zk的逆映射。

      *(pdTran + (int)(255 * dTemp)) = i;

将均衡化后的原图像中每一灰度级按照这种逆映射关系进行变换，即可得到直方图规定化后的结果。

直方图规定化方法Histst( )的完整实现如下。

      /**************************************************
      BOOL CImgProcess::Histst(CImgProcess* pTo, double* pdStdHist)
      功能：      图像的灰度直方图规定化方法
      参数：      CImgProcess * pTo：输出CImgProcess对象的指针
                  double * pdStdHist：标准直方图数组（要求已经归一化的直方图）
      返回值：    BOOL类型，true为成功，false为失败
      ***************************************************/
      BOOL CImgProcess::Histst(CImgProcess* pTo, double* pdStdHist)
      {
            int i, j;

            // 首先检查图像是否是8位灰度图像
            if (m_pBMIH->biBitCount! =8) return false;

            BYTE gray;            // 临时变量，存储当前光标像素的灰度值
            int target;           // 临时变量，存储当前光标像素的目标值

            double pdHist[256];   // 临时变量，存储灰度直方图
            this->GenHist(pdHist);

            double dTemp;     // 临时变量，存储累加的直方图数据
            int pdTran[256];  // 临时变量，存储标准直方图均衡化的变换矩阵
            memset(pdTran, -1, sizeof(int)*256);

            // 求标准直方图的均衡化变换矩阵
            for (i=0; i<256; i++)
            {
                dTemp = 0;
                for (BYTE k=0; k<i; k++)
                {
                      dTemp+=*(pdStdHist + k);
                }
                *(pdTran + (int)(0.5+255 * dTemp)) = i;
            }

            // 去除均衡化变换矩阵中的间断点——插值
            {
                i=0, j=0;
                while(i<255)
                {
                      if(*(pdTran + i + 1)! =-1)
                      {
                          i++;
                          continue;
                      }
                      j = 1;
                      while((*(pdTran + i + j)==-1)&&((i + j)<=255))
                      {
                          *(pdTran + i + j)=*(pdTran + i);
                          j++;
                      }
                }
            }

            // 对原图像首先进行灰度均衡化后再进行规定化
            for (i=0; i<m_pBMIH->biHeight; i++)
            {
                for (j=0; j<m_pBMIH->biWidth; j++)
                {
                      dTemp = 0;
                      gray = GetGray(j, i);

                      for (BYTE k=0; k<gray; k++)
                      {
                          dTemp+=*(pdHist + k);
                      };

                      target = *(pdTran + (int)(255 * dTemp));

                      if (target < 0) target = 0;
                      if (target > 255) target = 255;

                      // 写入目标图像
                      pTo->SetPixel(j, i, RGB(target, target, target));
                }
            };

            return true;
      }

本书同时也提供了Histst()的一个重载方法，它接受CImg类型指针作为待匹配直方图的标准图像，实现如下。

      /**************************************************
      BOOL CImgProcess::Histst(CImg* pTo, CImg* pStd)
      功能：      图像的灰度直方图规定化方法
      参数  :   CImgProcess * pTo：输出CImg对象的指针
                CImgProcess* pStd：标准目标图像
      返回值：    BOOL类型，TRUE为成功，FALSE为失败
      ***************************************************/
      BOOL CImgProcess::Histst(CImg* pTo, CImg* pStd)
      {
        // 标准图像直方图
        double pdStdHist[256];

        pStd->GenHist(pdStdHist);

        return Histst(pTo, pdStdHist);
      }

利用Histst()函数实现直方图规定化的完整示例被封装在DIPDemo工程中的视图类函数CDIPDemoView::OnPointHistst ()，其中调用Histst()函数的代码片断如下所示。

        // 输出的临时对象
        CImgProcess imgOutput = imgInput;

        // 直方图规定化
        imgInput.Histst(&imgOutput, pdStdHist);
        // stdImage为待匹配直方图的标准图像

        // 将结果返回给文档类
        pDoc->m_Image = imgOutput;

上述程序运行时会弹出对话框，要求用户选择一幅标准图像。读者可以通过光盘中示例程序DIPDemo中的菜单命令“点运算→直方图规定化”来观察处理效果。