数字图像处理算法实现 编程心得

摘要: 关于空间域图像处理算法框架，直方图处理，空间域滤波器算法框架的编程心得，使用GDI+(C++)

一，图像文件的读取

    初学数字图像处理时，图像文件的读取往往是一件麻烦的事情，我们要面对各种各样的图像文件格式，如果仅用C++的fstream库那就必须了解各种图像编码格式，这对于初学图像处理是不太现实的，需要一个能帮助轻松读取各类图像文件的库。在Win32平台上GDI+(C++)是不错的选择，不光使用上相对于Win32 GDI要容易得多，而且也容易移植到.Net平台上的GDI+。

    Gdiplus::Bitmap类为我们提供了读取各类图像文件的接口，Bitmap::LockBits方法产生的BitmapData类也为我们提供了高速访问图像文件流的途径。这样我们就可以将精力集中于图像处理算法的实现，而不用关心各种图像编码。具体使用方式请参考MSDN中GDI+文档中关于Bitmap类和BitmapData类的说明。另外GDI+仅在Windows XP/2003上获得直接支持，对于Windows 2000必须安装相关DLL，或者安装有Office 2003，Visual Studio 2003 .Net等软件。

二，空间域图像处理算法框架

 (1) 在空间域图像处理中，对于一个图像我们往往需要对其逐个像素的进行处理，对每个像素的处理使用相同的算法（或者是图像中的某个矩形部分）。即，对于图像f(x,y),其中0≤ x≤M,0≤y≤N,图像为M*N大小，使用算法algo，则f(x,y) = algo(f(x,y))。事先实现一个算法框架，然后再以函数指针或函数对象(functor,即实现operator()的对象)传入算法，可以减轻编程的工作量。

    如下代码便是一例 :

#ifndef PROCESSALGO_H #define PROCESSALGO_H   #include <windows.h> #include <Gdiplus.h>     namespace nsimgtk {          template <typename pixelType, Gdiplus::PixelFormat pixelFormat, class Processor>     bool ProcessPixelsOneByOne(Gdiplus::Bitmap* const p_bitmap, Processor processor, unsigned int x, unsigned int y,                                                            unsigned int width, unsigned int height)     {                    if (p_bitmap == NULL)                    {                             return false;                    }                      if ((width + x > p_bitmap->GetWidth()) || (height + y >p_bitmap->GetHeight()))                    {                             return false;                    }           Gdiplus::BitmapData bitmapData;                    Gdiplus::Rect rect(x, y, width,height);                 if (p_bitmap->LockBits(&rect, Gdiplus::ImageLockModeWrite, pixelFormat, &bitmapData) != Gdiplus::Ok)             {                             return false;                    }                      pixelType *pixels = (pixelType*)bitmapData.Scan0;                              for (unsigned int row=0; row<height; ++row)                    {                             for (unsigned int col=0; col<width; ++col)                             {                                      processor(&pixels[col+row*bitmapData.Stride/sizeof(pixelType)]);                                  }                    }                      if (p_bitmap->UnlockBits(&bitmapData) != Gdiplus::Ok)                    {                             return false;                    }                            return true;          } }   #endif

ProcessPixelsOneByOne函数可以对图像中从(x,y)位置起始，width*height大小的区域进行处理。模板参数pixelType用于指定像素大小，例如在Win32平台上传入unsigned char即为8位，用于8阶灰度图。模板参数Processor为图像处理算法实现，可以定义类实现void operator(pixelType *)函数，或者传入同样接口的函数指针。

    如下便是一些算法示例(说明见具体注释):

#ifndef SPATIALDOMAIN_H #define SPATIALDOMAIN_H #include <cmath> #include <string>   namespace nsimgtk {     // 8 阶灰度图的灰度反转算法            class NegativeGray8          {          public:                    void operator()(unsigned char *const p_value)                    {                             *p_value ^= 0xff;                    }          };         // 8 阶灰度图的 Gamma 校正算法          class GammaCorrectGray8          {          private:                    unsigned char d_s[256];          public:                    GammaCorrectGray8::GammaCorrectGray8(double c, double gamma);                      void operator()(unsigned char*const p_value)                    {                             *p_value = d_s[*p_value];                    }          };       // 8 阶灰度图的饱和度拉伸算法          class ContrastStretchingGray8          {          private:                    unsigned char d_s[256];          public:                    ContrastStretchingGray8::ContrastStretchingGray8(double a1, double b1, unsigned int x1,                             double a2, double b2, unsigned int x2, double a3, double b3);                      void operator()(unsigned char*const p_value)                    {                             *p_value = d_s[*p_value];                    }          };         // 8 阶灰度图的位平面分割，构造函数指定位平面号          class BitPlaneSliceGray8          {          private:                          unsigned char d_s[256];          public:                    BitPlaneSliceGray8(unsigned char bitPlaneNum);                      void operator()(unsigned char* const p_value)                    {                             *p_value = d_s[*p_value];                    }          }; }   #endif   // 上述类中各构造函数的实现代码，应该分在另一个文件中，此处为说明方便，一并列出 #include "SpatialDomain/spatialDomain.h"   namespace nsimgtk {          GammaCorrectGray8::GammaCorrectGray8(double c, double gamma)          {                    double temp;                    for (unsigned int i=0; i<256; ++i)                    {                             temp = ceil(c * 255.0 * pow(double(i)/255.0, gamma));                             d_s[i] = unsigned char(temp);                    }          }            ContrastStretchingGray8::ContrastStretchingGray8(double a1, double b1, unsigned int x1,                             double a2, double b2, unsigned int x2, double a3, double b3)          {                    if (x1 > 255 || x2 > 255 || x1 > x1)                    {                             for (unsigned int i=0; i<256; ++i)                                      d_s[i] = i;                    }                    else                    {                             double tmp;                             for (unsigned int i=0; i<x1; ++i)                             {                                      tmp = ceil(a1*double(i)+b1);                                      d_s[i] = (unsigned char)tmp;                             }                               for (unsigned int i=x1; i<x2; ++i)                             {                                      tmp = ceil(a2*double(i)+b2);                                      d_s[i] = (unsigned char)tmp;                             }                               for (unsigned int i=x2; i<256; ++i)                             {                                      tmp = ceil(a3*double(i)+b3);                                      d_s[i] = (unsigned char)tmp;                             }                    }          }            BitPlaneSliceGray8::BitPlaneSliceGray8(unsigned char bitPlaneNum)          {                   unsigned char bitMaskArray[8] =                    {                             0x01, 0x02, 0x04, 0x08,                             0x10, 0x20, 0x40, 0x80                    };                      for (unsigned int i=0; i<256; ++i)                    {                             unsigned char tmp = i;                             tmp &= bitMaskArray[bitPlaneNum];                             tmp = (tmp >> bitPlaneNum) * 255;                             d_s[i] = tmp;                   }          } }

(2) 直方图在GDI+1.0中没有获得支持，我们必须自行实现。直方图相关的处理在数字图像处理中占有重要地位，可以通过它获取图像灰度级的统计信息，且可以通过直方图进行一些重要的图像增强技术，如直方图均衡化，直方图规定化，基本全局门限等。

下面是获取8阶图像直方图的算法实现：

namespace nsimgtk {          bool GetHistogramNormalizeGray8(Gdiplus::Bitmap * const p_bitmap, float *histogramArray)          {                    if (p_bitmap == NULL || histogramArray == NULL)                    {                             return false;                    }                      Gdiplus::BitmapData bitmapData;                    Gdiplus::Rect rect(0, 0, p_bitmap->GetWidth(), p_bitmap->GetHeight());                      if (p_bitmap->LockBits(&rect, Gdiplus::ImageLockModeRead, PixelFormat8bppIndexed, &bitmapData) != Gdiplus::Ok)             {                             return false;                    }                      unsigned char *pixels = (unsigned char*)bitmapData.Scan0;         unsigned int histogram[256];                    for (int i=0; i<256; ++i)                    {                             histogram[i] = 0;                    }                      for (unsigned int row=0; row<p_bitmap->GetWidth(); ++row)                    {                             for (unsigned int col=0; col<p_bitmap->GetHeight(); ++col)                             {                                      ++histogram[pixels[col+row*bitmapData.Stride]];                             }                    }                      const unsigned int totalPixels = p_bitmap->GetWidth() * p_bitmap->GetHeight();                    for (int i=0; i<256; ++i)                    {                             histogramArray[i] = float(histogram[i]) / float(totalPixels);                    }                      if (p_bitmap->UnlockBits(&bitmapData) != Gdiplus::Ok)                    {                             return false;                    }                      return true;          } }

在获取直方图后(即上面算法的第二个参数)，再将其作为参数传入下面的对象的构造函数，然后以该对象为仿函数传入ProcessPixelsOneByOne即可实现8阶图像直方图均衡化:

#ifndef SPATIALDOMAIN_H #define SPATIALDOMAIN_H   #include <cmath> #include <string>   namespace nsimgtk {     // 8 阶灰度图的直方图均衡化          class HistogramEqualizationGray8          {          private:                    unsigned char d_s[256];          public:                    HistogramEqualizationGray8(const float *const histogramArray);                                       void operator()(unsigned char *const p_value)                    {                             *p_value = d_s[*p_value];                    }          }; }   #endif            // #include "SpatialDomain/spatialDomain.h"   namespace nsimgtk {          HistogramEqualizationGray8::HistogramEqualizationGray8(const float *const histogramArray)          {                    if (histogramArray != NULL)                    {                             float sum = 0.0;                             for (int i=0; i<256; ++i)                             {                                      sum += histogramArray[i];                                      d_s[i] = unsigned char(sum * 255);                             }                    }          } }

(3)空间域滤波器，滤波器是一个m*n大小的掩模，其中m,n均为大于1的奇数。滤波器逐像素地通过图像的全部或部分矩形区域，然后逐像素地对掩模覆盖下的像素使用滤波器算法获得响应，将响应赋值于当前像素即掩模中心像素，另外滤波器算法使用中将会涉及到图像边缘的问题，这可以对边缘部分掩模使用补零法补齐掩模下无像素值的区域，或者掩模的移动范围以不越出图像边缘的方式移动，当然这些处理方法都会给图像边缘部分带来不良效果，但是一般情况下，图像边缘部分往往不是我们关注的部分或者没有重要的信息。

下面的滤波器算法框架SpatialFilterAlgo即以补零法(zero-padding)实现:

#ifndef SPATIALFILTER_H #define SPATIALFILTER_H   #include <vector> #include <numeric> #include <algorithm> #include <gdiplus.h> #include <fstream> #include <cmath>   namespace nsimgtk {     template <typename pixelType, Gdiplus::PixelFormat pixelFormat, class FilterMask>     bool SpatialFilterAlgo(Gdiplus::Bitmap* const p_bitmap, FilterMask filterMask, unsigned int x, unsigned int y,                                                            unsigned int width, unsigned int height)     {                    if (p_bitmap == NULL)                    {                             return false;                    }                      if ((width + x > p_bitmap->GetWidth()) || (height + y >p_bitmap->GetHeight()))                    {                             return false;                    }           Gdiplus::BitmapData bitmapData;                    Gdiplus::Rect rect(x, y, width,height);                 if (p_bitmap->LockBits(&rect, Gdiplus::ImageLockModeWrite, pixelFormat, &bitmapData) != Gdiplus::Ok)             {                             return false;                    }                      pixelType *pixels = (pixelType*)bitmapData.Scan0;                            const unsigned int m = filterMask.d_m;                                         // mask's width         const unsigned int n = filterMask.d_n;                                          // mask's height         std::vector<pixelType> tmpImage((m-1+width)*(n-1+height));   // extend image to use zero-padding                 // copy original bitmap to extended image with zero-padding method         for (unsigned int row=0; row<height; ++row)                    {                             for (unsigned int col=0; col<width; ++col)                             {                                      tmpImage[(col+m/2)+(row+n/2)*(bitmapData.Stride/sizeof(pixelType)+m-1)] =                     pixels[col+row*bitmapData.Stride/sizeof(pixelType)];                                 }                    }                 // process every pixel with filterMask         for (unsigned int row=0; row<height; ++row)                    {                             for (unsigned int col=0; col<width; ++col)                             {                 // fill the "m*n" mask with the current pixel's neighborhood                 for (unsigned int i=0; i<n; ++i)                 {                     for (unsigned int j=0; j<m; ++j)                     {                         filterMask.d_mask[i*m+j] = tmpImage[(col+j)+(row+i)*(bitmapData.Stride/sizeof(pixelType)+m-1)];                     }                 }                               // replace the current pixel with filter mask's response                                      pixels[col+row*bitmapData.Stride/sizeof(pixelType)] = filterMask.response();                                   }                    }           if (p_bitmap->UnlockBits(&bitmapData) != Gdiplus::Ok)                    {                             return false;                    }                            return true;          } }   #endif

其中模板参数FilterMask即为滤波掩模算法。通常的滤波算法有均值滤波器，可以模糊化图像，去除图形中的细节部分，使得我们可以关注图像中较为明显的部分，均值滤波器用于周期性噪声。中值滤波器用于图像中存在椒盐噪声也即脉冲噪声的情况下。另外有基于一阶微分的Sobel梯度算子和基于两阶微分的拉普拉斯算子，它们往往被用于边缘检测中。

下面是一些滤波器算法的具体实现，所以滤波器算法都应该实现pixelType response ()函数以及有d_mask,d_m,d_n成员，这可以通过继承__filteMask类获得(不需要付出虚函数代价)。

#ifndef SPATIALFILTER_H #define SPATIALFILTER_H   #include <vector> #include <numeric> #include <algorithm> #include <gdiplus.h> #include <fstream> #include <cmath>   namespace nsimgtk {     // 滤波器掩模的基类，提供掩模大小 d_m,d_n ，掩模覆盖下的 m*n 个像素值 d_mask     // others filterMask should inherit it     template <typename pixelType>     struct __filterMask     {         const unsigned int d_m;         const unsigned int d_n;           // image's pixels under the m*n filter mask         std::vector<pixelType> d_mask;           // filter mask's width and heigh must be a odd, if not, it will plus one for the width or the height         __filterMask(unsigned int m, unsigned int n)             : d_m(m%2 ? m:m+1), d_n(n%2 ? n:n+1), d_mask(d_m*d_n)         {         }     };       // 掩模权值为全 1 的均值滤波器     template <typename pixelType>     class averagingFilterMaskSp         : public __filterMask<pixelType>     {     public:         averagingFilterMaskSp(unsigned int m, unsigned int n)             : __filterMask<pixelType>(m, n)         { }           pixelType response()         {             return std::accumulate(d_mask.begin(), d_mask.end(), 0) / (d_m * d_n);         }     };       // 可自定义掩模权值的均值滤波器     template <typename pixelType>     class averagingFilterMask         : public __filterMask<pixelType>     {     private:         std::vector<pixelType> d_weight;                 // weights' vector(m*n)         int d_weight_sum;                                        // all weights' sum       public:         averagingFilterMask(unsigned int m, unsigned int n, const std::vector<pixelType>& weightVec)             : __filterMask<pixelType>(m, n), d_weight(weightVec)         {             if (weightVec.size() != d_mask.size())             {                 // if weight's size isn't equal to mask's size, it will change filter mask as a special filter mask                 d_weight.resize(d_mask.size(), 1);             }               d_weight_sum = std::accumulate(d_weight.begin(), d_weight.end(), 0);         }           pixelType response()         {             return std::inner_product(d_mask.begin(), d_mask.end(), d_weight.begin(), 0) / d_weight_sum;         }     };       // 中值滤波器     template <typename pixelType>     class medianFilterMask         : public __filterMask<pixelType>     {     public:         medianFilterMask(unsigned int m, unsigned int n)             : __filterMask<pixelType>(m, n)         { }           pixelType response()         {             std::sort(d_mask.begin(), d_mask.end());             return d_mask[d_mask.size()/2];         }     };       // 3*3 拉普拉斯滤波器     // the mask is: [0 1 0           [0 -1 0     //             1 -5 1     or    -1 5 -1     //             0 1 0]          0 -1 0]     // if pixel's brightness is less than min, set it to min     // if pixel's brightness is larger than max, set it to max     template <typename pixelType, pixelType min, pixelType max>     class laplacianFilter         : public __filterMask<pixelType>     {     public:         laplacianFilter()             : __filterMask<pixelType>(3, 3)         { }           pixelType response()         {             int ret = (int)(5*(int)d_mask[4]) - ((int)d_mask[5]+d_mask[3]+d_mask[1]+d_mask[7]);             if (ret < min)                 ret = min;             if (ret > max)                 ret = max;             return ret;         }     };       // 3*3Sobel 滤波器     // the mask is: [-1 -2 -1            [-1 0 1     //            0 0 0    and       -2 0 2     //            1 2 1]             -1 0 1]     // if pixel's brightness is larger than max, set it to max     template <typename pixelType, pixelType max>     class sobelFilter         : public __filterMask<pixelType>     {     public:         sobelFilter()             : __filterMask<pixelType>(3, 3)         { }           pixelType response()         {             int ret = ::abs(d_mask[6]+2*d_mask[7]+d_mask[8]-d_mask[0]-2*d_mask[1]-d_mask[2])                 + ::abs(d_mask[2]+2*d_mask[5]+d_mask[8]-d_mask[0]-2*d_mask[3]-d_mask[6]);                         if (ret > max)                 ret = max;             return ret;         }     }; }   #endif

 

数字图像处理算法实现

                                 ------------ 编程心得 (1)

2001414 班 朱伟 20014123

摘要: 关于空间域图像处理算法框架，直方图处理，空间域滤波器算法框架的编程心得，使用GDI+(C++)

一，图像文件的读取

    初学数字图像处理时，图像文件的读取往往是一件麻烦的事情，我们要面对各种各样的图像文件格式，如果仅用C++的fstream库那就必须了解各种图像编码格式，这对于初学图像处理是不太现实的，需要一个能帮助轻松读取各类图像文件的库。在Win32平台上GDI+(C++)是不错的选择，不光使用上相对于Win32 GDI要容易得多，而且也容易移植到.Net平台上的GDI+。

    Gdiplus::Bitmap类为我们提供了读取各类图像文件的接口，Bitmap::LockBits方法产生的BitmapData类也为我们提供了高速访问图像文件流的途径。这样我们就可以将精力集中于图像处理算法的实现，而不用关心各种图像编码。具体使用方式请参考MSDN中GDI+文档中关于Bitmap类和BitmapData类的说明。另外GDI+仅在Windows XP/2003上获得直接支持，对于Windows 2000必须安装相关DLL，或者安装有Office 2003，Visual Studio 2003 .Net等软件。

二，空间域图像处理算法框架

 (1) 在空间域图像处理中，对于一个图像我们往往需要对其逐个像素的进行处理，对每个像素的处理使用相同的算法（或者是图像中的某个矩形部分）。即，对于图像f(x,y),其中0≤ x≤M,0≤y≤N,图像为M*N大小，使用算法algo，则f(x,y) = algo(f(x,y))。事先实现一个算法框架，然后再以函数指针或函数对象(functor,即实现operator()的对象)传入算法，可以减轻编程的工作量。

    如下代码便是一例 :

#ifndef PROCESSALGO_H #define PROCESSALGO_H   #include <windows.h> #include <Gdiplus.h>     namespace nsimgtk {          template <typename pixelType, Gdiplus::PixelFormat pixelFormat, class Processor>     bool ProcessPixelsOneByOne(Gdiplus::Bitmap* const p_bitmap, Processor processor, unsigned int x, unsigned int y,                                                            unsigned int width, unsigned int height)     {                    if (p_bitmap == NULL)                    {                             return false;                    }                      if ((width + x > p_bitmap->GetWidth()) || (height + y >p_bitmap->GetHeight()))                    {                             return false;                    }           Gdiplus::BitmapData bitmapData;                    Gdiplus::Rect rect(x, y, width,height);                 if (p_bitmap->LockBits(&rect, Gdiplus::ImageLockModeWrite, pixelFormat, &bitmapData) != Gdiplus::Ok)             {                             return false;                    }                      pixelType *pixels = (pixelType*)bitmapData.Scan0;                              for (unsigned int row=0; row<height; ++row)                    {                             for (unsigned int col=0; col<width; ++col)                             {                                      processor(&pixels[col+row*bitmapData.Stride/sizeof(pixelType)]);                                  }                    }                      if (p_bitmap->UnlockBits(&bitmapData) != Gdiplus::Ok)                    {                             return false;                    }                            return true;          } }   #endif

ProcessPixelsOneByOne函数可以对图像中从(x,y)位置起始，width*height大小的区域进行处理。模板参数pixelType用于指定像素大小，例如在Win32平台上传入unsigned char即为8位，用于8阶灰度图。模板参数Processor为图像处理算法实现，可以定义类实现void operator(pixelType *)函数，或者传入同样接口的函数指针。

    如下便是一些算法示例(说明见具体注释):

#ifndef SPATIALDOMAIN_H #define SPATIALDOMAIN_H #include <cmath> #include <string>   namespace nsimgtk {     // 8 阶灰度图的灰度反转算法            class NegativeGray8          {          public:                    void operator()(unsigned char *const p_value)                    {                             *p_value ^= 0xff;                    }          };         // 8 阶灰度图的 Gamma 校正算法          class GammaCorrectGray8          {          private:                    unsigned char d_s[256];          public:                    GammaCorrectGray8::GammaCorrectGray8(double c, double gamma);                      void operator()(unsigned char*const p_value)                    {                             *p_value = d_s[*p_value];                    }          };       // 8 阶灰度图的饱和度拉伸算法          class ContrastStretchingGray8          {          private:                    unsigned char d_s[256];          public:                    ContrastStretchingGray8::ContrastStretchingGray8(double a1, double b1, unsigned int x1,                             double a2, double b2, unsigned int x2, double a3, double b3);                      void operator()(unsigned char*const p_value)                    {                             *p_value = d_s[*p_value];                    }          };         // 8 阶灰度图的位平面分割，构造函数指定位平面号          class BitPlaneSliceGray8          {          private:                          unsigned char d_s[256];          public:                    BitPlaneSliceGray8(unsigned char bitPlaneNum);                      void operator()(unsigned char* const p_value)                    {                             *p_value = d_s[*p_value];                    }          }; }   #endif   // 上述类中各构造函数的实现代码，应该分在另一个文件中，此处为说明方便，一并列出 #include "SpatialDomain/spatialDomain.h"   namespace nsimgtk {          GammaCorrectGray8::GammaCorrectGray8(double c, double gamma)          {                    double temp;                    for (unsigned int i=0; i<256; ++i)                    {                             temp = ceil(c * 255.0 * pow(double(i)/255.0, gamma));                             d_s[i] = unsigned char(temp);                    }          }            ContrastStretchingGray8::ContrastStretchingGray8(double a1, double b1, unsigned int x1,                             double a2, double b2, unsigned int x2, double a3, double b3)          {                    if (x1 > 255 || x2 > 255 || x1 > x1)                    {                             for (unsigned int i=0; i<256; ++i)                                      d_s[i] = i;                    }                    else                    {                             double tmp;                             for (unsigned int i=0; i<x1; ++i)                             {                                      tmp = ceil(a1*double(i)+b1);                                      d_s[i] = (unsigned char)tmp;                             }                               for (unsigned int i=x1; i<x2; ++i)                             {                                      tmp = ceil(a2*double(i)+b2);                                      d_s[i] = (unsigned char)tmp;                             }                               for (unsigned int i=x2; i<256; ++i)                             {                                      tmp = ceil(a3*double(i)+b3);                                      d_s[i] = (unsigned char)tmp;                             }                    }          }            BitPlaneSliceGray8::BitPlaneSliceGray8(unsigned char bitPlaneNum)          {                   unsigned char bitMaskArray[8] =                    {                             0x01, 0x02, 0x04, 0x08,                             0x10, 0x20, 0x40, 0x80                    };                      for (unsigned int i=0; i<256; ++i)                    {                             unsigned char tmp = i;                             tmp &= bitMaskArray[bitPlaneNum];                             tmp = (tmp >> bitPlaneNum) * 255;                             d_s[i] = tmp;                   }          } }

(2) 直方图在GDI+1.0中没有获得支持，我们必须自行实现。直方图相关的处理在数字图像处理中占有重要地位，可以通过它获取图像灰度级的统计信息，且可以通过直方图进行一些重要的图像增强技术，如直方图均衡化，直方图规定化，基本全局门限等。

下面是获取8阶图像直方图的算法实现：

namespace nsimgtk {          bool GetHistogramNormalizeGray8(Gdiplus::Bitmap * const p_bitmap, float *histogramArray)          {                    if (p_bitmap == NULL || histogramArray == NULL)                    {                             return false;                    }                      Gdiplus::BitmapData bitmapData;                    Gdiplus::Rect rect(0, 0, p_bitmap->GetWidth(), p_bitmap->GetHeight());                      if (p_bitmap->LockBits(&rect, Gdiplus::ImageLockModeRead, PixelFormat8bppIndexed, &bitmapData) != Gdiplus::Ok)             {                             return false;                    }                      unsigned char *pixels = (unsigned char*)bitmapData.Scan0;         unsigned int histogram[256];                    for (int i=0; i<256; ++i)                    {                             histogram[i] = 0;                    }                      for (unsigned int row=0; row<p_bitmap->GetWidth(); ++row)                    {                             for (unsigned int col=0; col<p_bitmap->GetHeight(); ++col)                             {                                      ++histogram[pixels[col+row*bitmapData.Stride]];                             }                    }                      const unsigned int totalPixels = p_bitmap->GetWidth() * p_bitmap->GetHeight();                    for (int i=0; i<256; ++i)                    {                             histogramArray[i] = float(histogram[i]) / float(totalPixels);                    }                      if (p_bitmap->UnlockBits(&bitmapData) != Gdiplus::Ok)                    {                             return false;                    }                      return true;          } }

在获取直方图后(即上面算法的第二个参数)，再将其作为参数传入下面的对象的构造函数，然后以该对象为仿函数传入ProcessPixelsOneByOne即可实现8阶图像直方图均衡化:

#ifndef SPATIALDOMAIN_H #define SPATIALDOMAIN_H   #include <cmath> #include <string>   namespace nsimgtk {     // 8 阶灰度图的直方图均衡化          class HistogramEqualizationGray8          {          private:                    unsigned char d_s[256];          public:                    HistogramEqualizationGray8(const float *const histogramArray);                                       void operator()(unsigned char *const p_value)                    {                             *p_value = d_s[*p_value];                    }          }; }   #endif            // #include "SpatialDomain/spatialDomain.h"   namespace nsimgtk {          HistogramEqualizationGray8::HistogramEqualizationGray8(const float *const histogramArray)          {                    if (histogramArray != NULL)                    {                             float sum = 0.0;                             for (int i=0; i<256; ++i)                             {                                      sum += histogramArray[i];                                      d_s[i] = unsigned char(sum * 255);                             }                    }          } }

(3)空间域滤波器，滤波器是一个m*n大小的掩模，其中m,n均为大于1的奇数。滤波器逐像素地通过图像的全部或部分矩形区域，然后逐像素地对掩模覆盖下的像素使用滤波器算法获得响应，将响应赋值于当前像素即掩模中心像素，另外滤波器算法使用中将会涉及到图像边缘的问题，这可以对边缘部分掩模使用补零法补齐掩模下无像素值的区域，或者掩模的移动范围以不越出图像边缘的方式移动，当然这些处理方法都会给图像边缘部分带来不良效果，但是一般情况下，图像边缘部分往往不是我们关注的部分或者没有重要的信息。

下面的滤波器算法框架SpatialFilterAlgo即以补零法(zero-padding)实现:

#ifndef SPATIALFILTER_H #define SPATIALFILTER_H   #include <vector> #include <numeric> #include <algorithm> #include <gdiplus.h> #include <fstream> #include <cmath>   namespace nsimgtk {     template <typename pixelType, Gdiplus::PixelFormat pixelFormat, class FilterMask>     bool SpatialFilterAlgo(Gdiplus::Bitmap* const p_bitmap, FilterMask filterMask, unsigned int x, unsigned int y,                                                            unsigned int width, unsigned int height)     {                    if (p_bitmap == NULL)                    {                             return false;                    }                      if ((width + x > p_bitmap->GetWidth()) || (height + y >p_bitmap->GetHeight()))                    {                             return false;                    }           Gdiplus::BitmapData bitmapData;                    Gdiplus::Rect rect(x, y, width,height);                 if (p_bitmap->LockBits(&rect, Gdiplus::ImageLockModeWrite, pixelFormat, &bitmapData) != Gdiplus::Ok)             {                             return false;                    }                      pixelType *pixels = (pixelType*)bitmapData.Scan0;                            const unsigned int m = filterMask.d_m;                                         // mask's width         const unsigned int n = filterMask.d_n;                                          // mask's height         std::vector<pixelType> tmpImage((m-1+width)*(n-1+height));   // extend image to use zero-padding                 // copy original bitmap to extended image with zero-padding method         for (unsigned int row=0; row<height; ++row)                    {                             for (unsigned int col=0; col<width; ++col)                             {                                      tmpImage[(col+m/2)+(row+n/2)*(bitmapData.Stride/sizeof(pixelType)+m-1)] =                     pixels[col+row*bitmapData.Stride/sizeof(pixelType)];                                 }                    }                 // process every pixel with filterMask         for (unsigned int row=0; row<height; ++row)                    {                             for (unsigned int col=0; col<width; ++col)                             {                 // fill the "m*n" mask with the current pixel's neighborhood                 for (unsigned int i=0; i<n; ++i)                 {                     for (unsigned int j=0; j<m; ++j)                     {                         filterMask.d_mask[i*m+j] = tmpImage[(col+j)+(row+i)*(bitmapData.Stride/sizeof(pixelType)+m-1)];                     }                 }                               // replace the current pixel with filter mask's response                                      pixels[col+row*bitmapData.Stride/sizeof(pixelType)] = filterMask.response();                                   }                    }           if (p_bitmap->UnlockBits(&bitmapData) != Gdiplus::Ok)                    {                             return false;                    }                            return true;          } }   #endif

其中模板参数FilterMask即为滤波掩模算法。通常的滤波算法有均值滤波器，可以模糊化图像，去除图形中的细节部分，使得我们可以关注图像中较为明显的部分，均值滤波器用于周期性噪声。中值滤波器用于图像中存在椒盐噪声也即脉冲噪声的情况下。另外有基于一阶微分的Sobel梯度算子和基于两阶微分的拉普拉斯算子，它们往往被用于边缘检测中。

下面是一些滤波器算法的具体实现，所以滤波器算法都应该实现pixelType response ()函数以及有d_mask,d_m,d_n成员，这可以通过继承__filteMask类获得(不需要付出虚函数代价)。

#ifndef SPATIALFILTER_H #define SPATIALFILTER_H   #include <vector> #include <numeric> #include <algorithm> #include <gdiplus.h> #include <fstream> #include <cmath>   namespace nsimgtk {     // 滤波器掩模的基类，提供掩模大小 d_m,d_n ，掩模覆盖下的 m*n 个像素值 d_mask     // others filterMask should inherit it     template <typename pixelType>     struct __filterMask     {         const unsigned int d_m;         const unsigned int d_n;           // image's pixels under the m*n filter mask         std::vector<pixelType> d_mask;           // filter mask's width and heigh must be a odd, if not, it will plus one for the width or the height         __filterMask(unsigned int m, unsigned int n)             : d_m(m%2 ? m:m+1), d_n(n%2 ? n:n+1), d_mask(d_m*d_n)         {         }     };       // 掩模权值为全 1 的均值滤波器     template <typename pixelType>     class averagingFilterMaskSp         : public __filterMask<pixelType>     {     public:         averagingFilterMaskSp(unsigned int m, unsigned int n)             : __filterMask<pixelType>(m, n)         { }           pixelType response()         {             return std::accumulate(d_mask.begin(), d_mask.end(), 0) / (d_m * d_n);         }     };       // 可自定义掩模权值的均值滤波器     template <typename pixelType>     class averagingFilterMask         : public __filterMask<pixelType>     {     private:         std::vector<pixelType> d_weight;                 // weights' vector(m*n)         int d_weight_sum;                                        // all weights' sum       public:         averagingFilterMask(unsigned int m, unsigned int n, const std::vector<pixelType>& weightVec)             : __filterMask<pixelType>(m, n), d_weight(weightVec)         {             if (weightVec.size() != d_mask.size())             {                 // if weight's size isn't equal to mask's size, it will change filter mask as a special filter mask                 d_weight.resize(d_mask.size(), 1);             }               d_weight_sum = std::accumulate(d_weight.begin(), d_weight.end(), 0);         }           pixelType response()         {             return std::inner_product(d_mask.begin(), d_mask.end(), d_weight.begin(), 0) / d_weight_sum;         }     };       // 中值滤波器     template <typename pixelType>     class medianFilterMask         : public __filterMask<pixelType>     {     public:         medianFilterMask(unsigned int m, unsigned int n)             : __filterMask<pixelType>(m, n)         { }           pixelType response()         {             std::sort(d_mask.begin(), d_mask.end());             return d_mask[d_mask.size()/2];         }     };       // 3*3 拉普拉斯滤波器     // the mask is: [0 1 0           [0 -1 0     //             1 -5 1     or    -1 5 -1     //             0 1 0]          0 -1 0]     // if pixel's brightness is less than min, set it to min     // if pixel's brightness is larger than max, set it to max     template <typename pixelType, pixelType min, pixelType max>     class laplacianFilter         : public __filterMask<pixelType>     {     public:         laplacianFilter()             : __filterMask<pixelType>(3, 3)         { }           pixelType response()         {             int ret = (int)(5*(int)d_mask[4]) - ((int)d_mask[5]+d_mask[3]+d_mask[1]+d_mask[7]);             if (ret < min)                 ret = min;             if (ret > max)                 ret = max;             return ret;         }     };       // 3*3Sobel 滤波器     // the mask is: [-1 -2 -1            [-1 0 1     //            0 0 0    and       -2 0 2     //            1 2 1]             -1 0 1]     // if pixel's brightness is larger than max, set it to max     template <typename pixelType, pixelType max>     class sobelFilter         : public __filterMask<pixelType>     {     public:         sobelFilter()             : __filterMask<pixelType>(3, 3)         { }           pixelType response()         {             int ret = ::abs(d_mask[6]+2*d_mask[7]+d_mask[8]-d_mask[0]-2*d_mask[1]-d_mask[2])                 + ::abs(d_mask[2]+2*d_mask[5]+d_mask[8]-d_mask[0]-2*d_mask[3]-d_mask[6]);                         if (ret > max)                 ret = max;             return ret;         }     }; }   #endif