全國最多中醫師線上諮詢網站-台灣中醫網
發文 回覆 瀏覽次數:9291
推到 Plurk!
推到 Facebook!

Otsu & Kapur 演算法求最佳閥值

 
m58610
初階會員


發表:22
回覆:83
積分:36
註冊:2003-09-07

發送簡訊給我
#1 引用回覆 回覆 發表時間:2003-11-28 23:04:45 IP:211.76.xxx.xxx 未訂閱
先前作的程式..
不過不符合我使用..
[1] 鍾國亮,影像處理與電腦視覺,台灣東華書局股份有限公司,民國91年,pp.84-88
我也是參考上面這本書寫的...
他在我學校資工系有開課...所以學長推薦我買他的書
程式有Otsu跟Kapur演算法...
我想應該沒有錯誤...裡面有一些變異數..或是期望值方面的統計學不好懂
當時我有問學商的老師...
所以才寫出來的...提供給大家 發表人 - m58610 於 2003/11/28 23:31:34
附加檔案:Otsu&Kapur.exe
編輯記錄
m58610 重新編輯於 2008-11-04 04:06:21, 註解 2008-11-04 修正BUG‧
hahalin
版主


發表:295
回覆:1698
積分:823
註冊:2002-04-14

發送簡訊給我
#2 引用回覆 回覆 發表時間:2003-11-30 12:52:05 IP:211.76.xxx.xxx 未訂閱
交流交流 剛好這學期修統計學 變異數是 ((樣本值-平均數)^2*次數)/樣本數 標準差則是 (變異數)^(1/2) 期望值則是 樣本發生機率 * 樣本值 有說錯的還請高手不吝指正
m58610
初階會員


發表:22
回覆:83
積分:36
註冊:2003-09-07

發送簡訊給我
#3 引用回覆 回覆 發表時間:2003-11-30 18:26:44 IP:211.76.xxx.xxx 未訂閱
引言: 交流交流 剛好這學期修統計學 變異數是 ((樣本值-平均數)^2*次數)/樣本數 標準差則是 (變異數)^(1/2) 期望值則是 樣本發生機率 * 樣本值 有說錯的還請高手不吝指正
謝謝你... 是怎樣我倒是忘了.. 那時就是照著算式寫程式而已...ㄏㄏ
pusong
一般會員


發表:2
回覆:7
積分:1
註冊:2003-02-27

發送簡訊給我
#4 引用回覆 回覆 發表時間:2003-12-08 22:24:25 IP:140.122.xxx.xxx 未訂閱
m58610大大.... 可以問你幾個問題嗎??    關於您的OTSU的一些程式問題 一.  
 
for(int i=0;i<=255;i  )
    if(tau_w[i]==0) tau_w[i]=9999999999999999999999999;
  tau_min=tau_w[0];
  p_min=0;
  for(int i=1;i<=255;i  )
  {
    if(tau_min>tau_w[i])
    {
      tau_min=tau_w[i];
      p_min=i;    
tau_w[i]=9999999999999999999999999 為何取那麼大呢及為何取? tau_min=tau_w[0]; 是什麼意思呢? 二.
 
for(int i=0;i<=1000;i  ) p[i]=0;
  for(int y=0;yHeight;y  )
  {
    ptr=(Byte*)TheBitmap->ScanLine[y];
    for(int x=0;xWidth;x  )
    {
       e=ptr[x*3];
       p[e]  ;
e=ptr[x*3]; 為何要乘三呢? 三.
 
 tau_w[i]=tau1 tau2;    
您好像忘了乘上w1 w2 應是 tau_w[i]=tau1*w1 tau2*w2; 麻煩了~~~~~~
chenliyan163
一般會員


發表:30
回覆:30
積分:12
註冊:2003-09-15

發送簡訊給我
#5 引用回覆 回覆 發表時間:2004-12-27 16:07:29 IP:221.12.xxx.xxx 未訂閱
能否把这个转成delphi的, 万分感谢
阿土
一般會員


發表:2
回覆:6
積分:1
註冊:2003-12-02

發送簡訊給我
#6 引用回覆 回覆 發表時間:2005-02-16 15:18:50 IP:211.22.xxx.xxx 未訂閱
請問此作品發表是否不能下載.試了好幾次都不行
------
陳土焱
GK
一般會員


發表:0
回覆:2
積分:0
註冊:2006-08-12

發送簡訊給我
#7 引用回覆 回覆 發表時間:2006-08-18 10:48:59 IP:163.13.xxx.xxx 未訂閱

已經不能下載了,請問可以重新分享一次嗎?謝謝...

m58610
初階會員


發表:22
回覆:83
積分:36
註冊:2003-09-07

發送簡訊給我
#8 引用回覆 回覆 發表時間:2008-11-04 04:07:24 IP:140.118.xxx.xxx 未訂閱
修改BUG...
已經重新上傳...^^
killghost
一般會員


發表:14
回覆:21
積分:7
註冊:2004-04-21

發送簡訊給我
#9 引用回覆 回覆 發表時間:2009-07-31 10:25:27 IP:222.211.xxx.xxx 訂閱
把OTSU算法转成了Delphi了的了。测试应用程序中的图片。得到的值是一致的。

但是,测试了几个图片都是。(把图片转化为灰度后)得到的阀值居然是0或者是255 (很明显,阀值不应该是这个值)

故,非常怀疑上面的代码存在严重问题。如果有人需要,我可以贴出转成Delphi部分的代码。

下面是原始的一份OTSU的经典算法。


[code cpp]
/*
the following code was send by Ryan Dibble <dibbler@umich.edu>

The algorithm is very simple but works good hopefully.

Compare the grayscale histogram with a mass density diagram:
I think the algorithm is a kind of
divide a body into two parts in a way that the mass
centers have the largest distance from each other,
the function is weighted in a way that same masses have a advantage

TODO:
RGB: do the same with all colors (CMYG?) seperately
test: hardest case = two colors
bbg: test done, using a two color gray file. Output:
# OTSU: thresholdValue = 43 gmin=43 gmax=188
my changes:
- float -> double
- debug option added (vvv & 1..2)
- **image => *image, &image[i][1] => &image[i*cols 1]
Joerg.Schulenburg@physik.uni-magdeburg.de
ToDo:
- measure contrast
- detect low-contrast regions
*/
#include
#include
/*======================================================================*/
/* OTSU global thresholding routine */
/* takes a 2D unsigned char array pointer, number of rows, and */
/* number of cols in the array. returns the value of the threshold */
/*======================================================================*/
int
otsu (unsigned char *image, int rows, int cols, int x0, int y0, int dx, int dy, int vvv) {
unsigned char *np; // pointer to position in the image we are working with
int thresholdValue=1; // value we will threshold at
int ihist[256]; // image histogram
int i, j, k; // various counters
int n, n1, n2, gmin, gmax;
double m1, m2, sum, csum, fmax, sb;
// zero out histogram ...
memset(ihist, 0, sizeof(ihist));
gmin=255; gmax=0; k=dy/512 1;
// generate the histogram
for (i = 0; i < dy ; i =k) {
np = &image[(y0 i)*cols x0];
for (j = 0; j < dx ; j ) {
ihist[*np] ;
if(*np > gmax) gmax=*np;
if(*np < gmin) gmin=*np;
np ; /* next pixel */
}
}
// set up everything
sum = csum = 0.0;
n = 0;
for (k = 0; k <= 255; k ) {
sum = (double) k * (double) ihist[k]; /* x*f(x) mass moment */
n = ihist[k]; /* f(x) mass */
}
if (!n) {
// if n has no value we have problems...
fprintf (stderr, "NOT NORMAL thresholdValue = 160\n");
return (160);
}
// do the otsu global thresholding method
fmax = -1.0;
n1 = 0;
for (k = 0; k < 255; k ) {
n1 = ihist[k];
if (!n1) { continue; }
n2 = n - n1;
if (n2 == 0) { break; }
csum = (double) k *ihist[k];
m1 = csum / n1;
m2 = (sum - csum) / n2;
sb = (double) n1 *(double) n2 *(m1 - m2) * (m1 - m2);
/* bbg: note: can be optimized. */
if (sb > fmax) {
fmax = sb;
thresholdValue = k;
}
}
// at this point we have our thresholding value
// debug code to display thresholding values
if ( vvv & 1 )
fprintf(stderr,"# OTSU: thresholdValue = %d gmin=%d gmax=%d\n",
thresholdValue, gmin, gmax);
return(thresholdValue);
}
[/code]
編輯記錄
killghost 重新編輯於 2009-07-31 10:28:37, 註解 無‧
killghost 重新編輯於 2009-07-31 10:29:21, 註解 無‧
killghost 重新編輯於 2009-07-31 10:29:48, 註解 無‧
killghost
一般會員


發表:14
回覆:21
積分:7
註冊:2004-04-21

發送簡訊給我
#10 引用回覆 回覆 發表時間:2009-07-31 11:25:52 IP:222.211.xxx.xxx 訂閱
OTSU改良算法,稍后提供一份Delphi的。


[code cpp]
OTSU方法计算图像二值化的自适应阈值
/*
OTSU 算法可以说是自适应计算单阈值(用来转换灰度图像为二值图像)的简单高效方法。下面的代码最早由 Ryan Dibble提供,此后经过多人Joerg.Schulenburg, R.Z.Liu 等修改,补正。

转自:http://forum.assuredigit.com/display_topic_threads.asp?ForumID=8&TopicID=3480

算法对输入的灰度图像的直方图进行分析,将直方图分成两个部分,使得两部分之间的距离最大。划分点就是求得的阈值。

parameter: *image --- buffer for image
rows, cols --- size of image
x0, y0, dx, dy --- region of vector used for computing threshold
vvv --- debug option, is 0, no debug information outputed
*/
/*======================================================================*/
/* OTSU global thresholding routine */
/* takes a 2D unsigned char array pointer, number of rows, and */
/* number of cols in the array. returns the value of the threshold */
/*======================================================================*/
int otsu (unsigned char *image, int rows, int cols, int x0, int y0, int dx, int dy, int vvv)
{

unsigned char *np; // 图像指针
int thresholdValue=1; // 阈值
int ihist[256]; // 图像直方图,256个点

int i, j, k; // various counters
int n, n1, n2, gmin, gmax;
double m1, m2, sum, csum, fmax, sb;

// 对直方图置零...
memset(ihist, 0, sizeof(ihist));

gmin=255; gmax=0;
// 生成直方图
for (i = y0 1; i < y0 dy - 1; i ) {
np = ℑ[i*cols x0 1];
for (j = x0 1; j < x0 dx - 1; j ) {
ihist[*np] ;
if(*np > gmax) gmax=*np;
if(*np < gmin) gmin=*np;
np ; /* next pixel */
}
}

// set up everything
sum = csum = 0.0;
n = 0;

for (k = 0; k <= 255; k ) {
sum = (double) k * (double) ihist[k]; /* x*f(x) 质量矩*/
n = ihist[k]; /* f(x) 质量 */
}

if (!n) {
// if n has no value, there is problems...
fprintf (stderr, "NOT NORMAL thresholdValue = 160\n");
return (160);
}

// do the otsu global thresholding method
fmax = -1.0;
n1 = 0;
for (k = 0; k < 255; k ) {
n1 = ihist[k];
if (!n1) { continue; }
n2 = n - n1;
if (n2 == 0) { break; }
csum = (double) k *ihist[k];
m1 = csum / n1;
m2 = (sum - csum) / n2;
sb = (double) n1 *(double) n2 *(m1 - m2) * (m1 - m2);
/* bbg: note: can be optimized. */
if (sb > fmax) {
fmax = sb;
thresholdValue = k;
}
}

// at this point we have our thresholding value

// debug code to display thresholding values
if ( vvv & 1 )
fprintf(stderr,"# OTSU: thresholdValue = %d gmin=%d gmax=%d\n",
thresholdValue, gmin, gmax);

return(thresholdValue);
}

[/code]
cheng8210
一般會員


發表:2
回覆:3
積分:1
註冊:2010-12-21

發送簡訊給我
#11 引用回覆 回覆 發表時間:2010-12-27 16:52:58 IP:140.116.xxx.xxx 訂閱
我參考m大提供的程式碼再對照otsu的理論  修改了一些地方
分別使用修改前和修改後的程式去對同一張圖片算otsu值
發現顯示出來的黑白圖其實都差不多



void __fastcall TForm3::Button2Click(TObject *Sender)
{
Graphics::TBitmap *Bmp = new Graphics::TBitmap();
Bmp->Assign(Form1->Image1->Picture->Bitmap);
int e,p_min;
float w1,w2,u1,u2,sigma1,sigma2,sigma_w[256],p[256],sigma_min;
Byte* ptr;
//統計影像中灰階值為i的像素量
for(int i=0;i<256;i )
p[i]=0;
for(int y=0;yHeight;y )
{
ptr=(Byte*)Bmp->ScanLine[y];
for(int x=0;xWidth;x )
{
p[i]=p[i]/(Bmp->Width*Bmp->Height);
//分別計算門檻值為i時, w1, u1, sigma1
for(int i=0;i<256;i )
{
w1=0;
w2=0;
u1=0;
u2=0;
sigma1=0;
sigma2=0;
for(int j=0;j<=i;j )
w1=w1 p[j];
if(w1!=0)
{
for(int j=0;j<=i;j )
u1=u1 p[j]*j/w1;
for(int j=0;j<=i;j )
sigma1=sigma1 (j-u1)*(j-u1)*p[j]/w1;
}
for(int j=i 1;j<256;j )
w2=w2 p[j];
if(w2!=0)
{
for(int j=i 1;j<256;j )
u2=u2 p[j]*j/w2;
for(int j=i 1;j<256;j )
sigma2=sigma2 (j-u2)*(j-u2)*p[j]/w2;
}
sigma_min=99999999999999999;
p_min=0;
for(int i=0;i<256;i )
{
if(sigma_min>sigma_w[i])
{
sigma_min=sigma_w[i];
p_min=i;
}
}
Form3->Edit1->Text=IntToStr(p_min);
}
系統時間:2024-11-23 3:44:40
聯絡我們 | Delphi K.Top討論版
本站聲明
1. 本論壇為無營利行為之開放平台,所有文章都是由網友自行張貼,如牽涉到法律糾紛一切與本站無關。
2. 假如網友發表之內容涉及侵權,而損及您的利益,請立即通知版主刪除。
3. 請勿批評中華民國元首及政府或批評各政黨,是藍是綠本站無權干涉,但這裡不是政治性論壇!