kedoumy的专栏

每天进步一点点...

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

摘要: 关于空间域图像处理算法框架,直方图处理,空间域滤波器算法框架的编程心得,使用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

 
 
阅读更多
个人分类: 图像算法收集
想对作者说点什么? 我来说一句

没有更多推荐了,返回首页

加入CSDN,享受更精准的内容推荐,与500万程序员共同成长!
关闭
关闭