円のHough変換をC言語で

つくってみたのですが、ループが多すぎて実行時間が30分以上はかかりそうです。
ループをどうやって減らせばよいか、アドバイスお願いします。
円は(x-a)^2+(y-b)^2=400です。
以下、現在のプログラムの手順

1.b[x][y][a]を初期化(x,y,aでループ)
2.画素値が0のとき、b[x][y][a]=(円の式をb=～～の形にしたときの右辺)を行う(x,y,aでループ)
3.b[x][y][a]の頻度をカウントする。(カウントするbの値,x,y,aでループ。ここが時間かかる。)
4.逆変換する。まだ試してないけどxでしかループしないので時間はかからないはず。

No.3 回答者： osamuy 回答日時： 2013/07/07 22:43

下記サイトのプログラムと見比べてみると、

> 2. b[x][y][a]=(円の式をb=～～の形にしたときの右辺)を行う
の表現法がムダっぽいような。
x,y,r(半径)で計算してみては。

「Visual C++ 2010 Express を用いた易しい画像認識（２）」
http://homepage3.nifty.com/ishidate/vcpp10_r2/vc …

170x60くらいなら、C++で書いたプログラムで、1秒かからなかった。
スクリプト言語の中でもとくに遅いといわれる)rubyでも、1秒かからないですね。
# 640x480x(半径)40のデータでも10秒ほど。

参考URL：http://ideone.com/13V2Y3

No.2 回答者： duke_kimura 回答日時： 2013/07/07 20:53

コンパイラは何を使用されていますか？

ダミー変数を憶測で用意してコンパイルしてみましたがVS2010では通りませんでした。

また、色々細工をして「こうだろう」で動かした結果、定義や内容が分からない変数(rやb)がネックで、どうもうまく動作しませんでした。
できれば全てのソースコードと実行結果の一例が欲しいのですが、無理そうでしたら今回はお力になれそうにありません。

一応性能改善のアドバイスとしては、
・bh配列の初期化　→　memset関数で初期化することで高速化可能

・Hough変換部分　→　変数iは加算なのでsqrt部を事前計算後、iを加算しながらループすることで高コストの計算回数を抑えられます(但しiのループを内側に持ってくる必要あり)。このループが変数rとの兼ね合いがある場合は難しいです。

・カウント部分　→　少なくとも「カウントするbの値,x,y,aでループ」から「カウントするbの値」を排除するようにプログラムを構成する必要があります。例えば「x,y,aでループ」の中でうまいこと数えてやるようなイメージです。

この回答への補足

わざわざ試行錯誤までありがとうございます・・・。
cygwinのgccです。
全コードは字数オーバーになるので・・・
rは画素値を入れるポインタです。
bhはbの値です。bがもともと前の部分で使われていたので・・・

カウント部分の改善を試みているのですが、うまいこと数えられません・・・。

補足日時：2013/07/07 21:23

No.1 回答者： duke_kimura 回答日時： 2013/07/07 14:58

ソースコードは公開できませんか？

アルゴリズムを知らないので的確なアドバイスはできませんが、論理の最適化であればできます。

ちなみに30分っていう時間はどうやって推測したのでしょう？
250 * 1000 * 1000の配列で、
No2(初期化)、No3(条件分岐と演算)を擬似的にやってみましたが、5秒程で終わりました。
(勿論高コストな計算が入るならそれだけかかっても不思議じゃないですが)

この回答への補足

該当部分(1,2,3,4)を抜き出しました。
時間推測は、カウントの部分で、printf("debug\n");(今はかいてない)をループさせ、その表示間隔からしました。間隔からして50分ぐらいかかりそうでしたが、出力しない場合は早まると思って、30分ぐらいかな、と。
以下がソースコードです。よろしくお願いします・・・。
/*初期化*/
int a,pic[170][60];
int (*bh)[170][60][11200]=malloc(sizeof(*bh));
for(i=0 ; i<pPic->y ; i++){
for(j=0 ; j<pPic->x ; j++){
for(a=0;a<size;a++){
(*bh)[j][i][a]=0;
}
}
}
// printf("debug\n");
for(i=0 ; i<pPic->y ; i++){
for(j=0 ; j<pPic->x ; j++){
if(*r==0){
for(a=0;a<size;a++){
(*bh)[j][i][a]=(int)(-sqrt(400-pow((j-a),2))+i);//Hough変換
}
r++;
}
}
}
//printf("debug\n");
int b_atai,B_MAX=0,A_MAX=0,COUNT_MAX=0;
/*カウント*/
int count[size],ah[size];
for(b_atai=0;b_atai<=size;b++){
for(i=0 ; i<pPic->y ; i++){
for(j=0 ; j<pPic->x ; j++){
for(a=0;a<size;a++){
if(b_atai==(*bh)[j][i][a]){
count[b_atai]++;
ah[b_atai]=a;
}
}
}
}
if(count[b_atai]>=COUNT_MAX){
COUNT_MAX=count[b_atai];
B_MAX=b_atai;
A_MAX=ah[b_atai];
}
}
// printf("debug\n");
/*Hough逆変換*/
int y=0;
for(j=0;j<pPic->x;j++){
y=(int)(sqrt(400-pow((j-A_MAX),2))+B_MAX);
printf("%d %d\n",j,y);
}

補足日時：2013/07/07 15:14