//u,v,wの平均値の格納場所
double u_ave, v_ave, w_ave;
-
+ u_ave = v_ave = w_ave = 0;
+
//オプティカルフローの一回前の値の一時的な格納場所
- //double size_before = 0;
+ double size_before = 0;
//収束を見るためのもの
- //double conv = 0;
-
+ double conv = 0;
+ double density = 0;
+ double I = 0;
//最終的に求めたいもの
double u, v, w;
u = v = w = 0;
//Horn-Schunck法の3D拡張
int l;
float x, y, z;
- double size;
+ double size = 0;
// inputfileがある場合勾配の計算をする(Ix,Iy,Iz,Itの)
for (x = 1; x < in1.HeaderN.x - 1; x++) {
//u = u_ave - (Ix * (Ix * u_ave + Iy * v_ave + Iz * w_ave + It)) / ((1 / info.ALPHA) + Ix * Ix + Iy * Iy + Iz * Iz);
//v = v_ave - (Iy * (Ix * u_ave + Iy * v_ave + Iz * w_ave + It)) / ((1 / info.ALPHA) + Ix * Ix + Iy * Iy + Iz * Iz);
//w = w_ave - (Iz * (Ix * u_ave + Iy * v_ave + Iz * w_ave + It)) / ((1 / info.ALPHA) + Ix * Ix + Iy * Iy + Iz * Iz);
-
+
+ //1つ前の速度
+ size_before = size;
+
u = u_ave - (Ix * (Ix * u_ave + Iy * v_ave + Iz * w_ave + It)) / (info.ALPHA * info.ALPHA + Ix * Ix + Iy * Iy + Iz * Iz);
v = v_ave - (Iy * (Ix * u_ave + Iy * v_ave + Iz * w_ave + It)) / (info.ALPHA * info.ALPHA + Ix * Ix + Iy * Iy + Iz * Iz);
w = w_ave - (Iz * (Ix * u_ave + Iy * v_ave + Iz * w_ave + It)) / (info.ALPHA * info.ALPHA + Ix * Ix + Iy * Iy + Iz * Iz);
size = sqrt(u*u + v*v + w*w);
+ //breakのための差異計算
+ mrcPixelDataGet(&in1, x , y , z , &I, mrcPixelRePart, mrcPixelHowNearest);
+ conv = conv + fabs( (size - size_before) * I );
+ density = density + I;
+
//オプティカルフロー の差を出してどのぐらい収束するのか知るために...
if(info.flagSize1) {
if(l == 0){
}
}
printf("反復数%d回目\n", l);
+ if(conv/density < 0.0001)
+ break;
}