數位影像處理基本演算法實現（2）–section3.3長條圖處理

      從上一個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 imageMat & 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 pixelsfor (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; //bigint heightS = 200;//smallint heightS3 = 600;//3*smallif(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 backgroundm_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 verticalif(hist_graph[j].x < hist_graph[j+1].x){pix = saturate_cast<int>(hist_graph[j].y+k*(i-hist_graph[j].x));}else{//verticalpix = hist_graph[j].y;}hist[i] = pix;total_pix += pix;}//shift the section and checkif (++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 averagevector<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 sigma2float 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 offsetsint *box = new int[rgSize];for (int i=0; i<rgSize; ++i){box[i] = -rgSize/2+i;}//prepare varsint 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 assignfor (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！