原创 数字图像处理算法实现

2008-5-25 23:48 3574 7 7 分类: 处理器与DSP



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


2001414 朱伟 20014123


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


一,图像文件的读取


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


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


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


 (1) 在空间域图像处理中,对于一个图像我们往往需要对其逐个像素的进行处理,对每个像素的处理使用相同的算法(或者是图像中的某个矩形部分)。即,对于图像f(x,y),其中0xM,0yN,图像为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类也为我们提供了高速访问图像文件流的途径。这样我们就可以将精力集中于图像处理算法的实现,而不用关心各种图像编码。具体使用方式请参考MSDNGDI+文档中关于Bitmap类和BitmapData类的说明。另外GDI+仅在Windows XP/2003上获得直接支持,对于Windows 2000必须安装相关DLL,或者安装有Office 2003Visual Studio 2003 .Net等软件。


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


 (1) 在空间域图像处理中,对于一个图像我们往往需要对其逐个像素的进行处理,对每个像素的处理使用相同的算法(或者是图像中的某个矩形部分)。即,对于图像f(x,y),其中0xM,0yN,图像为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


 

发表于 @ 2004年12月27日 09:33:00|评论(1)|编辑


旧一篇: 空间域滤波器(1)


PARTNER CONTENT

文章评论0条评论)

登录后参与讨论
EE直播间
更多
我要评论
0
7
关闭 站长推荐上一条 /3 下一条