001: //-----------------------------------------------------------------
002: // optmdthr.c:
003: //    最大尤度しきい値選定法
004: //    (Dynamic Programming を用いた大津のしきい値選定法)
005: //             Last Update: <2004/12/02 21:51:52 A.Murakami>
006: //-----------------------------------------------------------------
007: #include <stdio.h>
008: #include <stdlib.h>
009: #include <string.h>
010: #include "opthmdthr.h"
011: //-----------------------------------------------------------------
012: // しきい値の選定
013: //-----------------------------------------------------------------
014: //  multi: 多値化の数
015: //  limit: 画素の最大値
016: //  nn   : ピクセル数
017: //  hist : ヒストグラム
018: // throrg: しきい値を保存するための配列
019: void calc_opthmdthr(int multi,int limit,int nn,int* hist,int* throrg)
020: {
021:     int *thr=(int*)calloc((multi+1)*(multi+1),sizeof(int));
022:     dp_ootu(thr,multi,limit,hist,nn);
023:     dp_thr2thr(throrg,thr,limit,multi);
024:     free(thr);
025: }
026: //-----------------------------------------------------------------
027: // しきい値の変換
028: //-----------------------------------------------------------------
029: void dp_thr2thr(int *throrg,int *thr,int limit,int multi)
030: {
031:     int i,pos,height=(limit-1)/(multi-1);
032:     *(throrg)=0; *(throrg+multi)=(limit-1);
033:     pos = *(thr+multi*(multi+1)+(multi-1));
034:     *(throrg+multi-1)=pos;
035:     for(i=multi-2;i>0;i--){
036:         pos = *(thr+multi*(multi+1)+i);
037:         *(throrg+i)=pos;
038:     }
039: }
040: //-----------------------------------------------------------------
041: // DP + 大津による判別しきい値法
042: //-----------------------------------------------------------------
043: double dp_ootu(int* thr,int multi,int limit,int* hist,int nn)
044: {
045:     //--------------------------------------------------
046:     // prob: limitの値が低く多値化数が大きいときに注意.
047:     //--------------------------------------------------
048:     int i,pos,level,th0,th1;
049:     int maxp;
050:     int *cthr  =(int*)calloc(multi+1,sizeof(int));
051:     int *cphist=(int*)calloc(limit+1,sizeof(int));
052:     int *item  =(int*)calloc((multi+1)*(limit+1),sizeof(int));
053:     double maxv=0.,maxcv,vtmp;
054:     double *w =(double*)calloc((limit+1),sizeof(double));
055:     double *mu=(double*)calloc((limit+1),sizeof(double));
056:     double *value=(double*)calloc((multi+1)*(limit+1),sizeof(double));
057:     // 初期化
058:     cphist[0]=0;
059:     memcpy(cphist+1,hist,sizeof(int)*(limit));
060:     // ヒストグラム情報の計算
061:     hist_sumprob(cphist,w,limit,nn);
062:     hist_summean(cphist,mu,limit,nn);
063:     //--------------------------------------------------
064:     // DP による多値化計算
065:     //--------------------------------------------------
066:     value=(double*)calloc((multi+1)*(limit+1),sizeof(double));
067:     for(level=1;level<=multi;level++){  // しきい値数の増加(多値化レベルk)
068:         cthr[0]=0; cthr[level]=limit;
069:         for(th0=level-1;th0<=limit;th0++){ // しきい値レベルの増加(濃度値)
070:             if(level==1){
071:                 cthr[level] = th0;
072:                 *(value + (level-1)*(limit+1) + th0) =
073:                      hist_calccdv(level,cthr,w,mu);
074:             } else {
075:                 maxcv = 0;
076:                 for(th1=0;th1<th0;th1++){  // 最大分散を示すしきい値の選択
077:                     cthr[level-1] = th1; cthr[level] = th0;
078:                     vtmp = *(value+(level-2)*(limit+1)+th1) +
079:                         hist_calccdv(level,cthr,w,mu);
080:                     if(vtmp > maxcv){
081:                         maxcv = vtmp; maxp = th1;
082:                     }
083:                 }
084:                 *(value+(level-1)*(limit+1)+th0) = maxcv;
085:                 *(item+(level-1)*(limit+1)+th0)  = maxp;
086:             }
087:         }
088:     }
089:     free(cphist); free(cthr);
090:     free(w); free(mu);
091:     //--------------------------------------------------
092:     // しきい値の更新
093:     //--------------------------------------------------
094:     pos = 0;
095:     *(thr+multi*(multi+1)+0)=0;
096:     *(thr+multi*(multi+1)+multi)=(limit-1);
097:     for(level=multi;level>=2;level--){
098:         pos=*(item+(level-1)*(limit+1)+limit);
099:         *(thr+level*(multi+1)+(level-1))=pos-1; // 位置修正
100:         for(i=level-2;i>0;i--){
101:             pos = *(item + i*(limit+1) + pos);
102:             *(thr+level*(multi+1)+i)=pos-1;     // 位置修正
103:         }
104:     }
105:     maxv = *(value+(multi-1)*(limit+1)+(limit));
106:     free(item); free(value);
107:     return maxv;
108: }
109: //-----------------------------------------------------------------
110: // 頻度確率累積量
111: //-----------------------------------------------------------------
112: void hist_sumprob(int *hist,double *w,int limit,int nn)
113: {
114:     int i,j;
115:     memset(w,0,sizeof(double)*(limit+1));
116:     for(j=0;j<=limit;j++) for(i=0;i<=j;i++)
117:         w[j]+=(double)hist[i]/(double)(nn);
118: }
119: //-----------------------------------------------------------------
120: // 平均累積量の計算
121: //-----------------------------------------------------------------
122: void hist_summean(int *hist,double *mu,int limit,int nn)
123: {
124:     int i,j;
125:     memset(mu,0,sizeof(double)*(limit+1));
126:     for(j=0;j<=limit;j++) for(i=0;i<=j;i++)
127:         mu[j]+=i*(double)hist[i]/(double)(nn);
128: }
129: //-----------------------------------------------------------------
130: // クラス間分散値の評価関数
131: // [The evaluation function of a classes distribution value: CDV]
132: //-----------------------------------------------------------------
133: double hist_calccdv(int cls,int *thr,double *w,double *mu)
134: {
135:     double tmpw,tmpmu;
136:     tmpw  = w[thr[cls]]-w[thr[cls-1]];
137:     tmpmu = (tmpw)?(mu[thr[cls]]-mu[thr[cls-1]])/tmpw : 0;
138:     return tmpw*tmpmu*tmpmu;
139: }