001: //-----------------------------------------------------------------
002: // hough.cpp:
003: //     Hough変換用関数 (直線検出*2, 円検出)
004: //                 Last Update: <2004/12/10 09:55:14 A.Murakami>
005: //-----------------------------------------------------------------
006: #include <stdio.h>
007: #include <stdlib.h>
008: #include <math.h>
009: #include <windows.h>
010: #include "wingui.h"
011: #include "ipcommon.h"
012: #include "hough.h"
013: //-----------------------------------------------------------------
014: extern UINT iHeight,iWidth,iLength,iSize;// 高さ,幅,１ラインの長さ,サイズ
015: extern LPBYTE lpOrgBMP,lpBMP;            // オリジナル画像, 表示画像
016: //-----------------------------------------------------------------
017: // for qsort
018: int cir_compare(const void *c1,const void *c2)
019: {
020:     return (((HC_POINT *)c2)->fv - ((HC_POINT *)c1)->fv);
021: }
022: // calc circle distance
023: int calc_cirdis(HC_POINT p1,HC_POINT p2)
024: {
025:     return rint(sqrt((p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y)));
026: }
027: //-----------------------------------------------------------------
028: // Hough 変換による検出結果の表示
029: //-----------------------------------------------------------------
030: void draw_hough_line(LPBYTE inCBuf,int x1,int y1,int x2,int y2)
031: {
032:     int xp,yp;
033:     double distance,step,t;
034:     int dx1,dx2,dy1,dy2;
035:     dx1 = x1; dx2 = x2; dy1 = y1; dy2 = y2;
036:     // 描画位置のクリッピング処理
037:     if(scrossln(&dx1,&y1,&x2,&y2)==0) return;
038:     // 描画開始
039:     t = 0.0;
040:     distance = calc_distance(dx1,dy1,dx2,dy2);
041:     step = 1.0/(1.5*distance);
042:     while(t<1.0){
043:         xp = rint(t*dx2 + (1.0-t)*dx1);
044:         yp = rint(t*dy2 + (1.0-t)*dy1);
045:         if(xp>=0 && xp<iWidth && yp>=0 && yp<iHeight){
046:             inCBuf[xp*3+yp*iLength + 0] = iWHITE;
047:             inCBuf[xp*3+yp*iLength + 1] = iWHITE;
048:             inCBuf[xp*3+yp*iLength + 2] = iWHITE;
049:         }
050:         t = t + step;
051:     }
052: }
053: // 円
054: void draw_hough_cir(LPBYTE inCBuf,int cir_n,LPPOINT cir_p,int radius)
055: {
056:     int i,ci,x,y;
057:     int th_max,*cir_x,*cir_y;
058:     double th_step;
059:     //-------------------------------------------------------
060:     // 円形フィルタの作成:[sin,cos 計算の高速化]
061:     //-------------------------------------------------------
062:     // １ピクセル毎の表示となるように
063:     th_max  = rint(2*M_PI*radius);
064:     th_step = 360./th_max;
065:     cir_x   = (int*)GlobalAlloc(GPTR,(th_max+1)*sizeof(int));
066:     cir_y   = (int*)GlobalAlloc(GPTR,(th_max+1)*sizeof(int));
067:     for(i=0;i<th_max;i++){
068:         cir_x[i] = rint(radius * cos(deg2rad(th_step*i)));
069:         cir_y[i] = rint(radius * sin(deg2rad(th_step*i)));
070:     }
071:     // 複製
072:     CopyMemory(inCBuf,lpOrgBMP,iLength*iHeight);
073:     // 円の表示
074:     for(ci=0;ci<cir_n;ci++) for(i=0;i<th_max;i++){
075:         x = cir_p[ci].x + cir_x[i];
076:         y = cir_p[ci].y + cir_y[i];
077:         if(x>=0 && x<iWidth && y>=0 && y<iHeight){
078:             inCBuf[x*3+y*iLength + 0] = iWHITE;
079:             inCBuf[x*3+y*iLength + 1] = iWHITE;
080:             inCBuf[x*3+y*iLength + 2] = iWHITE;
081:         }
082:     }
083:     // 後片付け
084:     GlobalFree(cir_x); GlobalFree(cir_y);
085: }
086: //-----------------------------------------------------------------
087: // Hough 変換実行: [画像中から直線の抽出]
088: //-----------------------------------------------------------------
089: void LHough_proc(LPBYTE inBuf,double lnMag,int lhough_opt)
090: {
091:     if(lhough_opt == HOUGH_LINE_MC){
092:         //--------------------------------------------------
093:         // m,c 空間への写像
094:         // x-y, y-x 座標での写像によりm->∞の問題に対応
095:         //--------------------------------------------------
096:         int m,c;                             // m:傾き*0.01, c:切片
097:         int skip_n = 1;                      // 写像スキップ数
098:         int hpx_limit = rint(iWidth*lnMag);  // 検出頻度
099:         int hpy_limit = rint(iHeight*lnMag); // 検出頻度
100:         int max_m=200,max_c=iHeight*2;       // 最大傾き,切片
101:         int hps_size=max_m * max_c;          // 写像空間の大きさ
102:         int upmax =rint(iHeight*1.5);        // 切片最大値
103:         int lowmax=rint(iHeight*0.5);        // 切片最小値
104:         int *hps_mcx,*hps_mcy;               // パラメータ空間
105:         int lx1,lx2,ly1,ly2;
106:         hps_mcx = (int*)GlobalAlloc(GPTR,hps_size*sizeof(int));
107:         hps_mcy = (int*)GlobalAlloc(GPTR,hps_size*sizeof(int));
108:         // パラメータ空間への写像
109:         Hough_line_mc(inBuf,hps_mcx,hps_mcy,max_m,max_c,skip_n);
110:         // 直線の描画
111:         for(m=0;m<max_m;m++) for(c=0;c<max_c;c++){
112:             if(hps_mcy[m*max_c+c] > hpy_limit){// 横方向
113:                 ly1 = c-lowmax;
114:                 ly2 = rint(0.01*(m-max_m/2)*iWidth+(c-lowmax));
115:                 draw_hough_line(lpBMP,0,ly1,iWidth,ly2);
116:             }
117:             if(hps_mcx[m*max_c+c] > hpx_limit){// 縦方向
118:                 lx1 = c-lowmax;
119:                 lx2 = rint(0.01*(m-max_m/2)*iHeight+(c-lowmax));
120:                 draw_hough_line(lpBMP,lx1,0,lx2,iHeight);
121:             }
122:         }
123:         GlobalFree(hps_mcx); GlobalFree(hps_mcy);
124:     } else {
125:         //--------------------------------------------------
126:         // ρ,θ 空間への写像
127:         //--------------------------------------------------
128:         int r,th;                            // r:距離, th:角度
129:         int skip_n = 1;                      // 写像スキップ数
130:         int hp_limit = rint(iWidth*lnMag);   // 検出頻度
131:         int max_r=rint(calc_distance(0,0,iWidth,iHeight));
132:         int max_th=180;                      // 最大距離,角度
133:         int hps_size=max_r * max_th;         // 写像空間の大きさ
134:         int *hps_rth;                        // パラメータ空間
135:         int ly1,ly2;
136:         hps_rth = (int*)GlobalAlloc(GPTR,hps_size*sizeof(int));
137:         // パラメータ空間への写像
138:         Hough_line_rth(inBuf,hps_rth,max_r,max_th,skip_n);
139:         // 検出直線の表示
140:         for(r=0;r<max_r;r++) for(th=0;th<max_th;th++){
141:             if(hps_rth[r*max_th+th] > hp_limit){
142:                 if(th==0){
143:                     ly1 = rint(r/sin(deg2rad(th)));
144:                     ly2 = rint(-1e32*iWidth+r/sin(deg2rad(th)));
145:                     draw_hough_line(lpBMP,0,ly1,iWidth,ly2);
146:                 } else if(th == 90){
147:                     ly1 = rint(r/sin(deg2rad(th)));
148:                     ly2 = rint(r/sin(deg2rad(th)));
149:                     draw_hough_line(lpBMP,0,ly1,iWidth,ly2);
150:                 } else {
151:                     ly1 = rint(r/sin(deg2rad(th)));
152:                     ly2 = rint(-1/tan(deg2rad(th))*iWidth
153:                                + r/sin(deg2rad(th)));
154:                     draw_hough_line(lpBMP,0,ly1,iWidth,ly2);
155:                 }
156:             }
157:         }
158:         GlobalFree(hps_rth);
159:     }
160: }
161: //-----------------------------------------------------------------
162: // Hough: m-cパラメータ空間への写像(直線検出)
163: //     m: 傾き, c: 切片
164: //-----------------------------------------------------------------
165: void Hough_line_mc(LPBYTE inBuf,int *hps_mcx,int *hps_mcy,
166:                    int max_m,int max_c,int step)
167: //-----------------------------------------------------------------
168: // inBuf: 白黒画像
169: // hps_mcx: HPS空間[x-y axis], hps_mcy: HPS空間[y-x axis]
170: // step: 写像スキップ数
171: //-----------------------------------------------------------------
172: {
173:     int x,y,m,c1,c2;
174:     int upmax=rint(iHeight*1.5),lowmax=rint(iHeight*0.5);
175:     int hps_size = max_m*max_c;
176:     // パラメータ空間の初期化
177:     FillMemory(hps_mcx,sizeof(int)*hps_size,0);
178:     FillMemory(hps_mcy,sizeof(int)*hps_size,0);
179:     //-------------------------------------------------------
180:     // Houghパラメータ空間(HPS:Hough Parameter Space)への写像
181:     //-------------------------------------------------------
182:     for(y=0;y<iHeight;y++) for(x=0;x<iWidth;x++){
183:         if(inBuf[x+y*iWidth] == iWHITE){
184:             for(m=0;m<max_m;m+=step){
185:                 // x-y 座標系から m-c 座標系への変換
186:                 c1 = rint(-0.01*(m-max_m/2)*x+y);
187:                 if(c1<upmax && c1>-lowmax)
188:                     hps_mcy[m*max_c+(c1+lowmax)] += 1;
189:                 // y-x 座標系から m-c 座標系への変換
190:                 c2 = rint(-0.01*(m-max_m/2)*y+x);
191:                 if(c2<upmax && c2>-lowmax)
192:                     hps_mcx[m*max_c+(c2+lowmax)] += 1;
193:             }
194:         }
195:     }
196: }
197: //-----------------------------------------------------------------
198: // Hough: ρ-θパラメータ空間への写像(直線検出)
199: //    ρ: 原点からの距離, θ: 角度
200: //-----------------------------------------------------------------
201: void Hough_line_rth(LPBYTE inBuf,int *rth,int max_r,int max_th,int thstep)
202: {
203:     int x,y,r,th,th_max;
204:     int hps_size = max_r*max_th;
205:     double *cos_x,*sin_y;
206:     // パラメータ空間の初期化
207:     FillMemory(rth,sizeof(int)*hps_size,0);
208:     //-------------------------------------------------------
209:     // 円形フィルタの作成:[sin,cos 計算の高速化]
210:     //-------------------------------------------------------
211:     // １ピクセル毎の計算となるように
212:     th_max = max_th;
213:     cos_x  = (double*)GlobalAlloc(GPTR,(th_max+1)*sizeof(double));
214:     sin_y  = (double*)GlobalAlloc(GPTR,(th_max+1)*sizeof(double));
215:     for(th=0;th<th_max;th++){
216:         cos_x[th] = cos(deg2rad(th));
217:         sin_y[th] = sin(deg2rad(th));
218:     }
219:     //-------------------------------------------------------
220:     // HPSへの写像
221:     //-------------------------------------------------------
222:     for(y=0;y<iHeight;y++) for(x=0;x<iWidth;x++){
223:         if(inBuf[x+y*iWidth] == iWHITE){
224:             // x-y 座標系から r-th 座標系への変換
225:             for(th=0;th<max_th;th+=thstep){
226:                 r = rint(x*cos_x[th]+y*sin_y[th]);
227:                 if(r<max_r && r>0) rth[r*max_th+th]++;
228:             }
229:         }
230:     }
231: }
232: //-----------------------------------------------------------------
233: // Hough 変換実行: [画像中から円の抽出]
234: //-----------------------------------------------------------------
235: int CHough_proc(LPBYTE inBuf,LPPOINT cpoints,
236:                 int radius,int FindMax,double lnMag)
237: {
238:     int i,ii,x,y,c_count,cf_count,hpsv;
239:     int hp_limit = rint(2*M_PI*radius*lnMag); // 検出頻度
240:     int near_dis = radius;                    // 最低近接距離
241:     HC_POINT p[1024];                         // 検出位置
242:     // Hough-パラメータ空間(HPS:Hough-Parameter Space)
243:     int* HPSwgt=(int*)GlobalAlloc(GPTR,iSize*sizeof(int));
244:     //-------------------------------------------------------
245:     // 初期化
246:     //-------------------------------------------------------
247:     for(i=0;i<1024;i++) p[i].flg = 1;
248:     FillMemory(HPSwgt,sizeof(int)*iSize,0);
249:     //-------------------------------------------------------
250:     // Hough: パラメータ空間への写像
251:     //-------------------------------------------------------
252:     Hough_cir(inBuf,HPSwgt,radius);
253:     //-------------------------------------------------------
254:     // 円の検出
255:     //-------------------------------------------------------
256:     c_count = 0;
257:     for(y=0;y<iHeight;y++) for(x=0;x<iWidth;x++) {
258:         if(c_count<1024){
259:             hpsv = HPSwgt[x+y*iWidth];
260:             if(hpsv > hp_limit){
261:                 p[c_count].x  = x; p[c_count].y  = y;
262:                 p[c_count].fv = hpsv;
263:                 c_count++;
264:             }
265:         }
266:     }
267:     //-------------------------------------------------------
268:     // 頻度の高い順にソート
269:     //-------------------------------------------------------
270:     qsort(p,c_count,sizeof(HC_POINT),cir_compare);
271:     //-------------------------------------------------------
272:     // 近接する円は除去
273:     //-------------------------------------------------------
274:     for(i=0;i<c_count;i++){
275:         if(p[i].flg==0) continue;
276:         for(ii=i+1;ii<c_count;ii++){
277:             if(p[ii].flg==0) continue;
278:             if(calc_cirdis(p[i],p[ii]) < near_dis){
279:                 p[ii].flg = 0;
280:             }
281:         }
282:     }
283:     cf_count = 0;
284:     for(i=0;i<c_count;i++){
285:         if(p[i].flg && i<FindMax){
286:             cpoints[cf_count].x = p[i].x;
287:             cpoints[cf_count].y = p[i].y;
288:             cf_count++;
289:         }
290:     }
291:     // 後片付け
292:     GlobalFree(HPSwgt);
293:     return cf_count;
294: }
295: //-----------------------------------------------------------------
296: // Hough: パラメータ空間への写像
297: //-----------------------------------------------------------------
298: void Hough_cir(LPBYTE inBuf,int* HPSwgt,int radius)
299: {
300:     int x,y,th,tx,ty;
301:     int yStep=1,xStep=1; // 元画像のスキップ数[n Skip毎に写像を行う]
302:     int cStep=1;         // 円形写像のスキップ
303:     int th_max,*cflt_x,*cflt_y;
304:     double cirstep;
305:     //-------------------------------------------------------
306:     // 円形フィルタの作成:[sin,cos 計算の高速化]
307:     //-------------------------------------------------------
308:     // １ピクセル毎の計算となるように
309:     th_max  = rint(2*M_PI*radius);
310:     cirstep = 360./th_max;
311:     cflt_x  = (int*)GlobalAlloc(GPTR,(th_max+1)*sizeof(int));
312:     cflt_y  = (int*)GlobalAlloc(GPTR,(th_max+1)*sizeof(int));
313:     for(th=0;th<th_max;th++){
314:         cflt_x[th] = rint(radius * cos(deg2rad(cirstep*th)));
315:         cflt_y[th] = rint(radius * sin(deg2rad(cirstep*th)));
316:     }
317:     //-------------------------------------------------------
318:     // パラメータ空間への写像
319:     //-------------------------------------------------------
320:     for(y=0;y<iHeight;y+=yStep) for(x=0;x<iWidth;x+=xStep){
321:         if(inBuf[x+y*iWidth]==iWHITE){
322:             for(th=0;th<th_max;th+=cStep){ // 円上
323:                 tx = x+cflt_x[th];
324:                 ty = y+cflt_y[th];
325:                 if(tx>=0 && tx<iWidth && ty>=0 && ty<iHeight){
326:                     HPSwgt[tx+ty*iWidth]+=1;
327:                 }
328:             }
329:         }
330:     }
331:     GlobalFree(cflt_x); GlobalFree(cflt_y);
332: }