从上一个section3.2到现在又过去一段时间了，断断续续终于把section3.3的直方图那些事弄出来了，还是那样：不求速度与健壮性，但求自己取实现一下，尽量自己去写代码，不用现成的函数。
《数字图像处理 Digital Image Processing-3E冈萨雷斯》(后面简称《DIP》)的section3.3是直方图处理有直方图均衡，直方图规则化，局部直方图处理以及用局部直直方图做图像增强，主要是这几个，所以这次就把这几个算法实现了一下，附带着写了一个简单的显示直方图的函数，能够以曲线或者条状显示直方图，1个通道或者3个通道。当然使用了一些OpenCV的函数。
算法内容都是按照书上来的，所以废话不说，先上声明：
enum HistType{Line_hist, Bin_hist};
class section_3_3
{
public:
//calculate image histogram;
bool getImgHist(Mat &src, vector< vector<int> > &hist, int chanels_use);
//equalize given image
Mat & equalizeImgHist(Mat &src);
//show given histogram--note! this change the histogram!
bool showHist(string wndname, vector< vector<int> > &hist, HistType histType, int thickness=1,bool bSave=false, string savename=string());
//normalize histogram to given range;
void normalizeHist(vector< vector<int> > &hist, int maxLeval);
//hist matching or hist normalization;
//use the point given histogram curve;
Mat& histMatching(Mat &src, vector<cv::Point> &hist_graph);
//use the given histogram;
Mat& histMatching(Mat &src, vector<int> &dest_hist);
//equalize use local hist;
Mat& localEqualize(Mat &src, int rgSize=3);
//image enhancement use local hist;
Mat& enhanceImg(Mat &src, int rgSize=3, float E=4.0f, float k0=0.4f, float k1=0.02f, float k2=0.4f);
private:
Mat m_destMat;
Mat m_showMat;
};
主要的就是前面说的那几个函数，采用了STL的vector来存放直方图，其中需要注意的是normalizeHist函数，其实就是将直方图的bin的最大值限制到maxLeval，方便直方图的显示，histMatching是直方图规则化或者叫直方图匹配，其中给出了两种实现，一个用vector点集指定直方图曲线的几个拐点，另一个是直接输入一个目标直方图，将源图像的直方图变换为与目标直方图一样(大体一样)。下面是函数实现--有点长
：
bool section_3_3::getImgHist(Mat &src, vector< vector<int> > &hist, int chanels_use)
{
if(src.channels() < chanels_use || !(chanels_use==1||chanels_use==3))
return false;
Mat temp_gray;
if(chanels_use ==1 && src.channels()==3)
cvtColor(src, temp_gray, COLOR_BGR2GRAY);
//create hist vectors;
vector<int> hh(256, 0);
for(int i=0; i<chanels_use; ++i)
{
hist.push_back(hh);
}
//calculate pixels
for (int i=0; i<src.rows; ++i)
{
uchar *src_ptr = src.ptr<uchar>(i);
for (int j=0; j<src.cols; ++j)
{
if(chanels_use==1)
{
++hist[0][src_ptr[j]];
}
else
{
int j3 = j*3;
++hist[0][src_ptr[j3]];
++hist[1][src_ptr[j3+1]];
++hist[2][src_ptr[j3+2]];
}
}
}
return true;
}
//normalize hist to the given range, maxLeval gives the upper bound, lower bound is 0;
void section_3_3::normalizeHist(vector< vector<int> > &hist, int maxLeval)
{
for (size_t i=0; i<hist.size(); ++i)
{
int maxvalue = 0;
for (size_t j=0; j<hist[i].size(); ++j)
{
if(hist[i][j]>maxvalue)
maxvalue = hist[i][j];
}
//if the maxvalue is less than maxLeval, no need to change;
if(maxvalue <= maxLeval)
continue;
//normalize the hist---note!-this change the given histogram;
for (size_t j=0; j<hist[i].size(); ++j)
{
hist[i][j] = (int)(hist[i][j]*maxLeval/(float)(maxvalue));
}
}
}
//show histogram in image, can choose hist type:bin or line, save it or not;
bool section_3_3::showHist(string wndname, vector< vector<int> > &hist, HistType histType, int thickness, bool bSave, string savename)
{
int width = 512+20;
int heightB = 400; //big
int heightS = 200; //small
int heightS3 = 600; //3*small
if(thickness > 2)
thickness = 2;
//curve hist ;
if(histType == Line_hist)
{
m_showMat.create(heightB+20, width, CV_8UC3);
normalizeHist(hist, heightB-10);
if(thickness==-1)
thickness = 1;
}
else //bin hist;
{
if(hist.size() == 3)
{
m_showMat.create(heightS3+20, width, CV_8UC3);
normalizeHist(hist, heightS-10);
}
else
{
m_showMat.create(heightB+20, width, CV_8UC3);
normalizeHist(hist, heightB-10);
}
}
//black background
m_showMat = Scalar::all(0);
//if bin hist and 3 channels, up to down=channel 0 to 2;
//if line hist all in one;
int chanel = hist.size();
for(int i=1; i<256;)
{
if (histType == Line_hist)
{
if(chanel == 1)
{
line(m_showMat, cv::Point(10+2*(i-1), heightB-hist[0][i-1]), cv::Point(10+2*i, heightB-hist[0][i]), Scalar::all(255), thickness);
}
else
{
line(m_showMat, cv::Point(10+2*(i-1), heightB-hist[0][i-1]), cv::Point(10+2*i, heightB-hist[0][i]), Scalar(255,0,0), thickness);
line(m_showMat, cv::Point(10+2*(i-1), heightB-hist[1][i-1]), cv::Point(10+2*i, heightB-hist[1][i]), Scalar(0,255,0), thickness);
line(m_showMat, cv::Point(10+2*(i-1), heightB-hist[2][i-1]), cv::Point(10+2*i, heightB-hist[2][i]), Scalar(0,0,255), thickness);
}
++i;
}
else
{
if (chanel == 1)
{
rectangle(m_showMat, cv::Point(10+2*(i-1), heightB),cv::Point(10+2*i, heightB-hist[0][i-1]),Scalar::all(255), thickness);
rectangle(m_showMat, cv::Point(10+2*i, heightB),cv::Point(10+2*(i+1), heightB-hist[0][i]),Scalar::all(255), thickness);
}
else
{
rectangle(m_showMat, cv::Point(10+2*(i-1), heightS),cv::Point(10+2*i, heightS-hist[0][i-1]),Scalar(255,0,0), thickness);
rectangle(m_showMat, cv::Point(10+2*i, heightS),cv::Point(10+2*(i+1), heightS-hist[0][i]),Scalar(255,0,0), thickness);
rectangle(m_showMat, cv::Point(10+2*(i-1), heightB),cv::Point(10+2*i, heightB-hist[1][i-1]),Scalar(0,255,0), thickness);
rectangle(m_showMat, cv::Point(10+2*i, heightB),cv::Point(10+2*(i+1), heightB-hist[1][i]),Scalar(0,255,0), thickness);
rectangle(m_showMat, cv::Point(10+2*(i-1), heightS3),cv::Point(10+2*i, heightS3-hist[2][i-1]),Scalar(0,0,255), thickness);
rectangle(m_showMat, cv::Point(10+2*i, heightS3),cv::Point(10+2*(i+1), heightS3-hist[2][i]),Scalar(0,0,255), thickness);
}
i+=2;
}
}
imshow(wndname, m_showMat);
if(bSave)
{
if(savename.empty())
savename = "default.bmp";
imwrite(savename, m_showMat);
}
return true;
}
//equalize image hist and give the changed image;
Mat & section_3_3::equalizeImgHist(Mat &src)
{
src.copyTo(m_destMat);
int channel = src.channels();
float k = 255.f/src.total();
//calculate hist;
vector< vector<int> > hist;
getImgHist(src, hist, channel);
//create map table;
vector<uchar> pixtable(256, 0);
vector< vector<uchar> > ttable;
for(int i=0; i<channel; ++i)
{
size_t accum = 0;
for (size_t j=0; j<pixtable.size(); ++j)
{
accum += hist[i][j];
pixtable[j] = static_cast<uchar>(k*accum);
}
ttable.push_back(pixtable);
}
//equalize src image;
for (int i=0; i<src.rows; ++i)
{
uchar *src_ptr = src.ptr<uchar>(i);
uchar *dest_ptr = m_destMat.ptr<uchar>(i);
for (int j=0; j<src.cols; ++j)
{
if(channel==1)
{
dest_ptr[j] = ttable[0][src_ptr[j]];
}
else
{
int j3 = j*3;
dest_ptr[j3] = ttable[0][src_ptr[j3]];
dest_ptr[j3+1] = ttable[1][src_ptr[j3+1]];
dest_ptr[j3+2] = ttable[2][src_ptr[j3+2]];
}
}
}
return m_destMat;
}
//hist matching---see <<digital image processing-3E-chs>>-section3.3.2,page:77-83;
//hist_graph--histogram curve given by multiplie points;
Mat& section_3_3::histMatching(Mat &src, vector<cv::Point> &hist_graph)
{
if(src.channels()!=1 || hist_graph.size()<2 || hist_graph[0].x!=0||hist_graph[hist_graph.size()-1].x!=255)
return m_destMat;
//check each point--every point's x coord should not less than the last one;
for (size_t i=1; i<hist_graph.size(); ++i)
{
if(hist_graph[i].x<hist_graph[i-1].x)
return m_destMat;
}
//create dest mat;
src.copyTo(m_destMat);
//equalize src image;
float k = 255.f/src.total();
//calculate hist;
vector< vector<int> > src_hist;
getImgHist(src, src_hist, 1);
//create src equalize table--index is r, value is s;
uchar r_s_table[256] = {0};
size_t accum = 0;
for (size_t i=0; i<256; ++i)
{
accum += src_hist[0][i];
r_s_table[i] = static_cast<uchar>(k*accum);
}
//create the required hist;
int hist[256]={0};
int uplimit = min(256, hist_graph[1].x);
size_t j=0;
size_t total_pix = 0;
for(int i=0;;)
{
float k = (hist_graph[j+1].y-hist_graph[j].y)/(float)(hist_graph[j+1].x-hist_graph[j].x);
for(; i<=uplimit; ++i)
{
int pix = 0;
//not vertical
if(hist_graph[j].x < hist_graph[j+1].x)
{
pix = saturate_cast<int>(hist_graph[j].y+k*(i-hist_graph[j].x));
}
else
{//vertical
pix = hist_graph[j].y;
}
hist[i] = pix;
total_pix += pix;
}
//shift the section and check
if (++j < hist_graph.size()-1)
{
uplimit = hist_graph[j+1].x;
}
else
{
for (;i<256; ++i)
{
hist[i] = 0;
}
break;
}
}
//create dest equalize table--index is z, value is s;
k = 255.f/total_pix;
uchar z_s_table[256] = {0};
accum = 0;
for (size_t i=0; i<256; ++i)
{
accum += hist[i];
z_s_table[i] = static_cast<uchar>(k*accum);
}
// cout<<"total="<<total_pix<<",accum="<<accum<<endl;
//create s_z table,init to -1 --index is s, value is z;
int s_z_table[256];
for(int t=0; t<256; ++t)
s_z_table[t] = -1;
//from big to small--if multiple z mapped to same s, the big overlayed by small z;
for (int i=255; i>-1; --i)
{
s_z_table[z_s_table[i]] = i;
// cout<<i<<"-"<<(int)z_s_table[i]<<endl;
}
//create the final transform table--index is r, value is z;
uchar r_z_table[256]={0};
uchar last_z = 0;
for (int i=0; i<256; ++i)
{
if (s_z_table[i] != -1)
{
last_z = r_z_table[i] = r_s_table[s_z_table[i]];
}
else
{
r_z_table[i] = last_z;
}
// cout<<i<<"-"<<(int)r_z_table[i]<<endl;
}
//tranform the src image pixels value;
for (int i=0; i<src.rows; ++i)
{
uchar *src_ptr = src.ptr<uchar>(i);
uchar *dest_ptr = m_destMat.ptr<uchar>(i);
for (int j=0; j<src.cols; ++j)
{
dest_ptr[j] = r_z_table[src_ptr[j]];
}
}
return m_destMat;
}
//match the image hist to the given hist, most codes are same with the one above;
Mat& section_3_3::histMatching(Mat &src, vector<int> &dest_hist)
{
if(src.channels()!=1 || dest_hist.size()!=256)
return m_destMat;
//create dest mat;
src.copyTo(m_destMat);
//equalize src image;
float k = 255.f/src.total();
//calculate hist;
vector< vector<int> > src_hist;
getImgHist(src, src_hist, 1);
//create src equalize table--index is r, value is s;
uchar r_s_table[256] = {0};
size_t accum = 0;
for (size_t i=0; i<256; ++i)
{
accum += src_hist[0][i];
r_s_table[i] = static_cast<uchar>(k*accum);
}
//create dest equalize table--index is z, value is s;
size_t total = 0;
for(int i=0; i<256; ++i)
total += dest_hist[i];
k = 255.f/total;
uchar z_s_table[256] = {0};
accum = 0;
for (size_t i=0; i<256; ++i)
{
accum += dest_hist[i];
z_s_table[i] = static_cast<uchar>(k*accum);
}
//create s_z table,init to -1 --index is s, value is z;
int s_z_table[256];
for(int t=0; t<256; ++t)
s_z_table[t] = -1;
//from big to small--if multiple z mapped to same s, the big overlayed by small z;
for (int i=255; i>-1; --i)
{
s_z_table[z_s_table[i]] = i;
// cout<<i<<"-"<<(int)z_s_table[i]<<endl;
}
//create the final transform table--index is r, value is z;
uchar r_z_table[256]={0};
uchar last_z = 0;
for (int i=0; i<256; ++i)
{
if (s_z_table[i] != -1)
{
last_z = r_z_table[i] = r_s_table[s_z_table[i]];
}
else
{
r_z_table[i] = last_z;
}
// cout<<i<<"-"<<(int)r_z_table[i]<<endl;
}
//tranform the src image pixels value;
for (int i=0; i<src.rows; ++i)
{
uchar *src_ptr = src.ptr<uchar>(i);
uchar *dest_ptr = m_destMat.ptr<uchar>(i);
for (int j=0; j<src.cols; ++j)
{
dest_ptr[j] = r_z_table[src_ptr[j]];
}
}
return m_destMat;
}
//use local hist to equalize image hist;--see<<digital image processing>>(call it DIP from now on),section 3.3.3;
//note! rgSize should less than 17, and this function is slow if rgSize is big!!
Mat& section_3_3::localEqualize(Mat &src, int rgSize)
{
//only deals odds that less than 17;
//if rgSize>16 then transtable will bigger than 256, this method will have no attraction;
if (src.channels()!=1 || (rgSize%2 != 1) || rgSize>15)
return m_destMat;
/*
the equalize formula: Sk=((L-1)/M*N)*sigma_i=0~k(Ri);
to this M*N=rgSize*rgSize, L=256, and i at most can be numbers,others are 0, so this formula equals(take 3*3 as example):
Sk = a*i; a=255/9, i=1~9; so we can create a table to accelerate the speed;
*/
int tabSize = rgSize*rgSize;
uchar *trans_table = new uchar[tabSize];
memset(trans_table, 0, tabSize);
float k = 255.f/tabSize;
for (int i=0; i<tabSize; ++i)
{
trans_table[i] = static_cast<uchar>(k*(i+1));
// cout<<i<<"="<<int(trans_table[i])<<endl;
}
//create the dest mat;
src.copyTo(m_destMat);
//
int *box = new int[rgSize];
for (int i=0; i<rgSize; ++i)
{
box[i] = -rgSize/2+i;
}
int step = src.step[0];
int offset_i=0,offset_j=0;
int count=0;
for (int i=0; i<src.rows; ++i)
{
for (int j=0; j<src.cols; ++j)
{
uchar pix = *(src.data+i*step+j);
count = 0;
for (int a=0; a<rgSize; ++a)
{
for (int b=0;b<rgSize; ++b)
{
offset_i = std::min(std::max(0,i+box[a]), src.rows-1);
offset_j = std::min(std::max(0,j+box[b]), src.cols-1);
if(*(src.data+step*offset_i+offset_j) <= pix)
++count;
}
}
*(m_destMat.data+i*step+j) = trans_table[count];
}
}
return m_destMat;
}
//image enhancement use local hist, <DIP>section 3.3.4;
Mat& section_3_3::enhanceImg(Mat &src, int rgSize, float E, float k0, float k1, float k2)
{
if (src.channels()!=1 || (rgSize !=3 && rgSize != 5))
return m_destMat;
//get histogram of total image;
vector< vector<int> > hist;
getImgHist(src, hist, 1);
//calculate average
vector<float> pers(256,0.f);
float m_g = 0.f;
float totals = src.total();
for (int i=0; i<256; ++i)
{
pers[i] = hist[0][i]/totals;
m_g += i*pers[i];
}
//calculate sigma2
float sigma2_g = 0.f;
for (int i=0; i<256; ++i)
{
sigma2_g += (i-m_g)*(i-m_g)*pers[i];
}
cout<<"mg="<<m_g<<",sigma2_g="<<sigma2_g<<endl;
//create the dest mat;
src.copyTo(m_destMat);
//roi offsets
int *box = new int[rgSize];
for (int i=0; i<rgSize; ++i)
{
box[i] = -rgSize/2+i;
}
//prepare vars
int roiSize = rgSize*rgSize;
int step = src.step[0];
int offset_i=0,offset_j=0;
float count = 0;
float k0_mg = k0*m_g;
float k1_sigma2_g = k1*sigma2_g;
float k2_sigma2_g = k2*sigma2_g;
//caluc and assign
for (int i=0; i<src.rows; ++i)
{
for (int j=0; j<src.cols; ++j)
{
uchar pix = *(src.data+i*step+j);
count = 0.f;
//calculate local avg;
for (int a=0; a<rgSize; ++a)
{
for (int b=0;b<rgSize; ++b)
{
offset_i = std::min(std::max(0,i+box[a]), src.rows-1);
offset_j = std::min(std::max(0,j+box[b]), src.cols-1);
count += *(src.data+step*offset_i+offset_j);
}
}
float m = count/roiSize;
float sigma2 = 0.f;
//calculate local sigma2;
for (int a=0; a<rgSize; ++a)
{
for (int b=0;b<rgSize; ++b)
{
offset_i = std::min(std::max(0,i+box[a]), src.rows-1);
offset_j = std::min(std::max(0,j+box[b]), src.cols-1);
uchar vl = *(src.data+step*offset_i+offset_j);
sigma2 += (vl-m)*(vl-m);
}
}
sigma2 /=roiSize;
//check conditions;
if(m <= k0_mg && sigma2<=k2_sigma2_g && sigma2>=k1_sigma2_g)
*(m_destMat.data+i*step+j) = MIN(saturate_cast<uchar>(E*pix), 255);
else
*(m_destMat.data+i*step+j) = pix;
if(0 && i>230&&j>190&&i<270&&j<230)
cout<<j<<","<<i<<",m="<<m<<",sigma="<<sigma2<<endl;
}
}
return m_destMat;
}
对应声明，按顺序来的。贴几个处理的图像：
原始的dog图像，以lena直方图为目标直方图变换后的图像。
对应的变换前，变换后，以及lena的直方图--从直方图可以看出，变换后的直方图确实与lena的很相似：
然后是局部直方图均衡化，对lena图像使用7*7邻域；
最后是直方图用于图像增强，对钨丝图像，采用的参数为邻域为3*3，E=4,k0=0.35,k1=0.0001,k2=0.05,下面是变换前后的图
虽然黑白交界处出现了轮廓，但是对右边很暗的那部分，确实出现了与书中一样的效果。
继续努力，下一个section！
扫一扫
198
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
