OTSU 算法学习，二值化，MFC，BMP
学习了一下OTSU算法，写成程序，加深理解。
这个算法是用来计算最优全局阈值的，所以对光照不均问题似乎不适用。按算法的原理，背景灰度有大的梯度时，会被分割成前景和背景。
对背景一致的图片的分割比较好，其实就是图片整体明暗变化时能自动调节阈值。好比白天到黑夜，阈值自动变。
程序如下：
// XBMPView.cpp : implementation of the CXBMPView class // #include "stdafx.h" // SHARED_HANDLERS can be defined in an ATL project implementing preview, thumbnail // and search filter handlers and allows sharing of document code with that project. #ifndef SHARED_HANDLERS #include "XBMP.h" #endif #include "XBMPDoc.h" #include "XBMPView.h" #ifdef _DEBUG #define new DEBUG_NEW #endif // CXBMPView IMPLEMENT_DYNCREATE(CXBMPView, CScrollView) BEGIN_MESSAGE_MAP(CXBMPView, CScrollView) // Standard printing commands ON_COMMAND(ID_FILE_PRINT, &CScrollView::OnFilePrint) ON_COMMAND(ID_FILE_PRINT_DIRECT, &CScrollView::OnFilePrint) ON_COMMAND(ID_FILE_PRINT_PREVIEW, &CXBMPView::OnFilePrintPreview) ON_WM_CONTEXTMENU() ON_WM_RBUTTONUP() END_MESSAGE_MAP() // CXBMPView construction/destruction CXBMPView::CXBMPView() { // TODO: add construction code here } CXBMPView::~CXBMPView() { } BOOL CXBMPView::PreCreateWindow(CREATESTRUCT& cs) { // TODO: Modify the Window class or styles here by modifying // the CREATESTRUCT cs return CScrollView::PreCreateWindow(cs); } // CXBMPView drawing void CXBMPView::OnDraw(CDC* /*pDC*/) { CXBMPDoc* pDoc = GetDocument(); ASSERT_VALID(pDoc); if (!pDoc) return; ShowBMP((void*)pDoc->m_Data, (int)pDoc->m_DataLen,0); ShowBMPGrayM((void*)pDoc->m_Data, (int)pDoc->m_DataLen); // TODO: add draw code for native data here } void CXBMPView::OnInitialUpdate() { CScrollView::OnInitialUpdate(); CSize sizeTotal; // TODO: calculate the total size of this view sizeTotal.cx = sizeTotal.cy = 100; SetScrollSizes(MM_TEXT, sizeTotal); } // CXBMPView printing void CXBMPView::OnFilePrintPreview() { #ifndef SHARED_HANDLERS AFXPrintPreview(this); #endif } BOOL CXBMPView::OnPreparePrinting(CPrintInfo* pInfo) { // default preparation return DoPreparePrinting(pInfo); } void CXBMPView::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/) { // TODO: add extra initialization before printing } void CXBMPView::OnEndPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/) { // TODO: add cleanup after printing } void CXBMPView::OnRButtonUp(UINT /* nFlags */, CPoint point) { ClientToScreen(&point); OnContextMenu(this, point); } void CXBMPView::OnContextMenu(CWnd* /* pWnd */, CPoint point) { #ifndef SHARED_HANDLERS theApp.GetContextMenuManager()->ShowPopupMenu(IDR_POPUP_EDIT, point.x, point.y, this, TRUE); #endif } // CXBMPView diagnostics #ifdef _DEBUG void CXBMPView::AssertValid() const { CScrollView::AssertValid(); } void CXBMPView::Dump(CDumpContext& dc) const { CScrollView::Dump(dc); } CXBMPDoc* CXBMPView::GetDocument() const // non-debug version is inline { ASSERT(m_pDocument->IsKindOf(RUNTIME_CLASS(CXBMPDoc))); return (CXBMPDoc*)m_pDocument; } #endif //_DEBUG // CXBMPView message handlers //BMP文件模式 void CXBMPView::ShowBMP(void* Data, int Len, int k) { unsigned char* pData; int DataLen; pData = (unsigned char*)Data; DataLen = Len; if (DataLen == 0) return; //if (DataLen > 10000) DataLen = 10000; BITMAPFILEHEADER hdr; LPBITMAPINFOHEADER lpbi; BITMAPINFOHEADER bmpInfo;//信息头 memcpy(&hdr, pData, 14); memcpy(&bmpInfo, pData + 14, 40); unsigned char* p = NULL; unsigned char* pBmpdata = NULL; p = (unsigned char*)pData + 14; pBmpdata = (unsigned char*)(pData + hdr.bfOffBits); lpbi = (LPBITMAPINFOHEADER)(p); /*CSize sizeTotal; sizeTotal.cx = lpbi->biWidth; sizeTotal.cy = lpbi->biHeight; SetScrollSizes(MM_TEXT, sizeTotal);*/ CRect rt; GetClientRect(&rt); CDC* pDC = GetDC(); pDC->SetStretchBltMode(COLORONCOLOR); int width = 0; int pos1 = 0; width = rt.Width() / 3; if (k == 0) { pos1 = 0; } else { pos1 = width; } StretchDIBits(pDC->GetSafeHdc(), pos1, 0, width, rt.Height(), 0, 0, lpbi->biWidth, lpbi->biHeight, pBmpdata, (BITMAPINFO*)&bmpInfo, DIB_RGB_COLORS, SRCCOPY); } //BMP文件模式 void CXBMPView::ShowBMPGrayM(void* Data, int Len) { unsigned char* pData; unsigned char* pData2; int DataLen; pData = (unsigned char*)Data; DataLen = Len; if (DataLen == 0) return; pData2 = new unsigned char[DataLen]; memcpy(pData2, (void*)pData, DataLen); //复制一份 BITMAPFILEHEADER hdr; LPBITMAPINFOHEADER lpbi; BITMAPINFOHEADER bmpInfo;//信息头 memcpy(&hdr, pData, 14); memcpy(&bmpInfo, pData + 14, 40); unsigned char* p = NULL; unsigned char* pBmpdata = NULL; p = (unsigned char*)pData + 14; pBmpdata = (unsigned char*)(pData + hdr.bfOffBits); lpbi = (LPBITMAPINFOHEADER)(p); /*CSize sizeTotal; sizeTotal.cx = lpbi->biWidth; sizeTotal.cy = lpbi->biHeight; SetScrollSizes(MM_TEXT, sizeTotal);*/ unsigned char* pBmpdata2 = NULL; pBmpdata2 = (unsigned char*)(pData2 + hdr.bfOffBits); int bfWidth, bfHeight; int biBitCount; int bytesPerLine; int skip; int imageSize; //DBI 数据信息 bfWidth = lpbi->biWidth; bfHeight = lpbi->biHeight; biBitCount = lpbi->biBitCount; //内存数据映射 bytesPerLine = ((bfWidth * biBitCount + 31) >> 5) << 2;//bfWidth *3+skip skip = 4 - ((bfWidth * biBitCount) >> 3) & 3; imageSize = bytesPerLine * bfHeight; int R, G, B, Gray; int door = 80; int tdoor[4][4] = //这是之前手工分区二值化用的，没有删掉 { 50, 60, 70, 80, 40, 50, 60, 70, 30, 40, 50, 60, 20, 30, 40, 50 }; int GK[256]; //各灰度值的像素点数 int GKSum = 0; //像素点总数 int GKMax1 = 0; //峰值 int GKMax2 = 0; for (int i = 0; i < 256; i++) GK[i] = 0; //处理内存数据 for (int i = 0; i < bfHeight; i++) //高度是反的 { for (int j = 0; j < bytesPerLine; j += 3) //宽度 转换为 行（列数） { B = pBmpdata[i * bytesPerLine + j]; G = pBmpdata[i * bytesPerLine + j + 1]; R = pBmpdata[i * bytesPerLine + j + 2]; Gray = (R * 19595 + G * 38469 + B * 7472) >> 16; GK[Gray]++; GKSum++; if (GK[Gray] > GKMax2) GKMax2 = GK[Gray]; if(GKMax2> GKMax1)GKMax1 = GKMax2; //door = tdoor[i * 4 / bfHeight][j * 4 / bytesPerLine]; //if (Gray > door) Gray = 255; //else Gray = 0; pBmpdata[i * bytesPerLine + j] = Gray; pBmpdata[i * bytesPerLine + j + 1] = Gray; pBmpdata[i * bytesPerLine + j + 2] = Gray; } } CString str; CString str2; double pk[256]; //灰度的概率 for (int i = 0; i < 256; i++) { pk[i] = (double)GK[i] / (double)GKSum; } double mg = 0; //图片的平均灰度值 for (int i = 0; i < 256; i++) { mg = mg + pk[i] * i; } PRINT("mg=%f", mg); int th; //阈值 double m1[256]; //对应每个阈值的 C1 均值 double m2[256]; //对应每个阈值的 C2 均值 double p1[256]; //p1[i]+p2[i] == 1 double p2[256]; for (th = 0; th < 256; th++) { m1[th] = 0; m2[th] = 0; p1[th] = 0; p2[th] = 0; for (int i = 0; i < th; i++) { m1[th] = m1[th] + pk[i] * i; p1[th] = p1[th] + pk[i]; } m1[th] = m1[th] / p1[th]; for (int i = th; i < 256; i++) { m2[th] = m2[th] + pk[i] * i; p2[th] = p2[th] + pk[i]; } m2[th] = m2[th] / p2[th]; } str = "P1+P2 == 1:"; //验证 P1+P2 == 1 for (int i = 0; i < 256; i++) { str2.Format(";%f+%f=%f", p1[i], p2[i],p1[i]+p2[i]); str += str2; } PRINT(str); str = "p1*m1+ p2*m2 - mg ==0:";//验证 p1*m1+ p2*m2 == mg for (int i = 0; i < 256; i++) { str2.Format(";%f", p1[i]*m1[i]+ p2[i]*m2[i]- mg); str += str2; } PRINT(str); double a[256];//类间方差 double amax; int thmax; amax = 0; thmax = 0; for (int i = 0; i < 256; i++) { a[i] = 0; //a[i]=p1[i] * p2[i] * (m1[i] - m2[i])* (m1[i] - m2[i]); a[i] = p1[i] * (m1[i] - mg) * (m1[i] - mg) + p2[i] * (m2[i] - mg) * (m2[i] - mg); if (a[i] > amax) { amax = a[i]; thmax = i; } } str = "a[i]:"; for (int i = 0; i < 256; i++) { str2.Format(";%d %f", i, a[i]); str += str2; } PRINT(str); PRINT("th=%d, amax= %f", thmax, amax); //pDC->Rectangle(t0 + thmax, 200, t0 + thmax + 1, 200 - GK[thmax] * 100 / GKMax1); //GKMax1！=0 //处理内存数据 for (int i = 0; i < bfHeight ; i++) //高度是反的 { for (int j = 0; j < bytesPerLine ; j += 3) //宽度 转换为 行（列数） { B = pBmpdata2[i * bytesPerLine + j]; G = pBmpdata2[i * bytesPerLine + j + 1]; R = pBmpdata2[i * bytesPerLine + j + 2]; Gray = (R * 19595 + G * 38469 + B * 7472) >> 16; if (Gray > thmax) Gray = 255; else Gray = 0; pBmpdata2[i * bytesPerLine + j] = Gray; pBmpdata2[i * bytesPerLine + j + 1] = Gray; pBmpdata2[i * bytesPerLine + j + 2] = Gray; } } //在右边显示灰度直方图 CRect rt; GetClientRect(&rt); CDC* pDC = GetDC(); int t0 = rt.Width()/3*2 + 10; int high = rt.Height() / 3 * 2; for (int i = 0; i < 256; i++) { pDC->Rectangle(t0 + i, high, t0 + i + 1, high - GK[i] * high / GKMax1); //GKMax1！=0 } PRINT("m1=%d,m2=%d, mg=%d", GKMax1, GKMax2, GKSum); pDC->Rectangle(t0 + thmax, high+20, t0 + thmax + 1, 0);
//在1号位置显示二值化后的图 ShowBMP(pData2, DataLen, 1); delete pData2; }
