opencv实现otsu阈值分割

#include <cv.h> #include <cxcore.h> #include <highgui.h> using namespace cv; #include <string> using std::string; #include <cstring> #include <iostream> using std::cout; using std::cin; using std::cerr; using std::endl; using std::ios; using std::right; using std::left; #include <fstream> using std::ofstream; #include <locale> using std::locale;//对该类还不太明白，有待学习 #include <iomanip> using std::setw; //得到图像直方图 Mat getHist(const Mat & img) { //Mat *arrays=&imgGray; const Mat *arrays=&img; int narrays=1; int channels[]={0};//单通道图像计算直方图 MatND hist; int dims=1; int histSize[]={256}; float hranges[] = { 0, 255 }; const float *ranges[] = { hranges }; //调用 calcHist 计算直方图, 结果存放在 hist 中 calcHist(arrays, narrays, channels, Mat(), hist, dims, histSize, ranges); normalize(hist,hist,1,0,cv::NORM_L1);//直方图标准化,归一化到[0,1] return hist; } double getThreshold_otus(const Mat &hist)//1维直方图,hist类型为32FC1 { const int n_level=hist.rows; double threshold=0;//要得到的otus阈值 double mean[256]={0},pro[256]={0}; double max_sigma=-1000,sigmaTemp=0; for(int k=0;k<n_level;k++) { double mean_sum=0; for(int i=0;i<k;i++) { pro[k]+=hist.at<float>(i,0); mean[k]+=i*hist.at<float>(i,0); } for(int i=0;i<n_level;i++) mean_sum+=i*hist.at<float>(i,0); sigmaTemp=((mean_sum*pro[k]-mean[k])*(mean_sum*pro[k]-mean[k]))/(pro[k]*(1-pro[k])); if(sigmaTemp>max_sigma) { max_sigma=sigmaTemp; threshold=k; } } return threshold; } int main() { char file[40]="JPCNN001.bmp"; const int imgFileNum=154; for(int i=1;i<=imgFileNum;i++) { sprintf(file,"F:\\DATA\\JPCNN1024\\JPCNN%03d.bmp",i);//批量读入文件名"F:\\DATA\\JPCNN1024\\JPCNN%03d.bmp" cout<<file<<endl; string filename(file);//将字符数组转换为string类型字符串 int location=filename.find_last_of("\\"); string imgName=filename.substr(location+1,strlen("JPCNN001.bmp"));//从路径中获得图片名称 Mat img=imread(filename);//读入图像,flag=1(默认)为3通道,8UC3 //Mat img; //simg.convertTo(img, CV_32F);//转换图像类型 if(img.empty()) { cout<<"load image error!"<<endl; fflush(stdin); getchar(); exit(0); } Mat imgGray;//定义单通道图像，用来将彩色图像转换为单通道灰度图像 if(3==img.channels()) cvtColor(img,imgGray,CV_BGR2GRAY);//三通道RGB图像转换为单通道灰度图像,8UC3-8UC1 else imgGray=img; Mat imgGray_inv; imgGray_inv=-imgGray+255;//图像反转 /*for(int i=100 ; i!=108 ; ++i) { cout<<endl; for( int j=100 ; j!= 108 ; ++j) cout<<(int)imgGray_inv.at<uchar>(i,j)<<" "; }*/ /******************图像直方图均衡****************/ Mat imgEqualize; equalizeHist(imgGray, imgEqualize); //输出图像区域 //cv::namedWindow(imgName,0); //cvResizeWindow(imgName.c_str(),512,512); //cv::imshow(imgName,imgEqualize); //cv::waitKey(0); /***********************************图像直方图操作*********************************/ Mat imgHist=getHist(imgEqualize); /****************************************otsu阈值分割********************************/ double threshold=0; threshold=getThreshold_otus(imgHist); //图像二值化 Mat imgThresh; //cv::threshold(imgGray,imgThresh,threshold,255,CV_THRESH_TOZERO_INV); cv::threshold(imgEqualize,imgThresh,threshold,255,CV_THRESH_TOZERO_INV); /*************************利用水平垂直投影找出肺部所在矩形区域**************************/ int rightLungLeft=0,rightLungRight=0; int leftLungLeft=0,leftLungRight=0; int top=0,bottom=0,middle=0; /***************************水平投影找肺左右边界*************************/ //Mat tempImg; //imgGray.convertTo(tempImg, CV_32F);//转换图像类型 //Mat_<double> imgGray2(imgGray.size()); const int rows=imgGray_inv.rows; const int cols=imgGray_inv.cols; Mat_<double> imgHpro(1,cols);//存储水平投影数值 for(int col=0;col<cols;col++) { cv::Scalar s=cv::sum(imgGray_inv.col(col)); imgHpro(0,col)=s[0]/cols; //cout<<s[0]<<" "; //if(0==col%10) cout<<endl; } Mat maskLpro( 1 , cols , CV_8UC1 , cv::Scalar(0) );//用模板限制找最小值区间,左边界 Mat maskRpro( 1 , cols , CV_8UC1 , cv::Scalar(0) );//用模板限制找最小值区间,右边界 Mat maskMpro( 1 , cols , CV_8UC1 , cv::Scalar(0) );//用模板限制找最小值区间,中线 for(int col = (int)((1.0/20)*cols) ; col < (int)((1.0/4)*cols) ; col++) maskLpro.at<uchar>(0,col)=1; for(int col = (int)((7.0/10)*cols) ; col<(int)((49.0/50)*cols) ; col++) maskRpro.at<uchar>(0,col)=1; for(int col = (int)((9.0/20)*cols) ; col < (int)((13.0/20)*cols) ; col++) maskMpro.at<uchar>(0,col)=1; double minValue_pro=0; Point minLoc_Lpro(0,0); Point minLoc_Rpro(0,0); Point minLoc_Mpro(0,0); minMaxLoc(imgHpro, &minValue_pro, 0, &minLoc_Lpro, 0, maskLpro);//找左边界 minMaxLoc(imgHpro, &minValue_pro, 0, &minLoc_Rpro, 0, maskRpro);//找右边界 minMaxLoc(imgHpro, &minValue_pro, 0, &minLoc_Mpro, 0, maskMpro);//找中线 rightLungLeft=minLoc_Lpro.x; leftLungRight=minLoc_Rpro.x; middle=minLoc_Mpro.x; rightLungRight=middle;//-(leftLungRight-rightLungLeft)/12;//经验值 leftLungLeft=middle;//+(leftLungRight-rightLungLeft)/12;//经验值 //cout<<endl<<middle<<endl; /*****************************输出肺部区域*****************************/ top=1,bottom=rows;//先不确定顶部和底部 //输出整个肺部区域 cv::Rect lung(0,0,0,0); lung.x=rightLungLeft; lung.y=top; lung.width=leftLungRight-rightLungLeft; lung.height=bottom-lung.y; if( (lung.width+lung.x) >=cols ) lung.width=cols-lung.x; if((lung.height+lung.y) >= rows ) lung.height=rows-lung.y; //输出右肺 cv::Rect rightLung(0,0,0,0); rightLung.x=rightLungLeft;//列 rightLung.y=top;//行 rightLung.width=rightLungRight-rightLungLeft; rightLung.height=bottom-rightLung.y; if((rightLung.width+rightLung.x)>=cols) rightLung.width=cols-rightLung.x; if((rightLung.height+rightLung.y)>=rows) rightLung.height=rows-rightLung.y; //输出左肺 cv::Rect leftLung(0,0,0,0); leftLung.x=leftLungLeft;//列 leftLung.y=top;//行 leftLung.width=leftLungRight-leftLungLeft; leftLung.height=bottom-leftLung.y; if((leftLung.width+leftLung.x)>=cols) leftLung.width=cols-leftLung.x; if((leftLung.height+leftLung.y)>=rows) leftLung.height=rows-leftLung.y; Mat imgRightLung=imgThresh(rightLung);//右肺 Mat imgLeftLung=imgThresh(rightLung);//左肺 Mat outImgLung1=imgGray_inv(rightLung); string outImg1="F:\\C++\\opencv\\otus阈值分割\\otus_segmentation\\otus_segmentation\\outRightLung\\"+imgName; cv::threshold(outImgLung1 , outImgLung1 , 250 , 0 , CV_THRESH_TOZERO_INV); outImgLung1=-outImgLung1+255; Mat lungEqualize; equalizeHist(outImgLung1, lungEqualize); imwrite( outImg1 , lungEqualize ); //Mat outImgLung2=imgGray_inv(rightLung); //cv::threshold(outImgLung2 , outImgLung2 , 250 , 0 , CV_THRESH_TOZERO_INV); //string outImg2="F:\\C++\\opencv\\otus阈值分割\\otus_segmentation\\otus_segmentation\\outRightLung1\\"+imgName; //imwrite( outImg2 , outImgLung2 ); const int rightLungRows=imgRightLung.rows; const int leftLungRows=imgLeftLung.rows; const int rightLungCols=imgRightLung.cols; const int leftLungCols=imgLeftLung.cols; /***********************对otsu阈值分割后的左右肺分别垂直投影找肺底部边界******************/ Mat_<double> imgVpro(1,rightLungRows);//存储垂直投影数值 int bottomRight=0; int bottomLeft=0; //右肺投影找底部 for(int row=0;row<rightLungRows;row++) { cv::Scalar s=cv::sum(imgRightLung.row(row)); imgVpro(0,row)=s[0]/rightLungCols; } Mat maskBpro(1,rightLungRows,CV_8UC1,cv::Scalar(0));//用模板限制找最大值区间,底部边界 minValue_pro=0; Point minLoc_Tpro(0,0); for(int row=(int)(0.5*rightLungRows);row<(int)((19.0/20)*rightLungRows);row++) maskBpro.at<uchar>(0,row)=1; minMaxLoc(imgVpro, &minValue_pro, 0, &minLo