src/Tools/mrcImage/mrcImageHornSchunck/src/mrcImageHornSchunck.c

   1 /*
   2 # mrcImageHornSchunck : $Revision$
   3 # $Date$
   4 # Created by $Author$
   5 # Usage : mrcImageHornSchunck
   6 # Attention
   7 #   $Loccker$
   8 #       $State$
   9 #
  10 */
  11 #include <stdio.h>
  12 #include <stdlib.h>
  13 #include <string.h>
  14 #include <math.h>
  15 #define GLOBAL_DECLARATION
  16 #include "../inc/config.h"
  17
  18 #define DEBUG
  19 #include "genUtil.h"
  20 #include "mrcImage.h"
  21 #include "mrcImageHornSchunck.h"
  22
  23
  24 /////////////////////////コード描き始め////////////////////////
  25
  26
  27 //#define WIDTH  128  // ボリュームの幅
  28 //#define HEIGHT 128  // ボリュームの高さ
  29 //#define DEPTH  128  // ボリュームの奥行き
  30 //#define SCALE_FACTOR 1.0  // スケールファクタ
  31 //#define MAX_ITERATIONS 10  // 反復回数
  32 //#define ALPHA 1.0  // パラメータ α
  33
  34
  35
  36 int main(int argc, char* argv[])
  37 {
  38     mrcImageHornSchunckInfo info;
  39     mrcImage in1;
  40         mrcImage in2;
  41         mrcImage in3;
  42         mrcImage outIx;
  43         mrcImage outIy;
  44         mrcImage outIz;
  45         mrcImage outIt;
  46         mrcImage outU;
  47         mrcImage outV;
  48         mrcImage outW;
  49         mrcImage Size;
  50
  51         mrcImage Size1;
  52         mrcImage Size2;
  53         mrcImage Size3;
  54         mrcImage Size4;
  55         mrcImage Size5;
  56         mrcImage Size6;
  57         mrcImage Size7;
  58         mrcImage Size8;
  59         mrcImage Size9;
  60         mrcImage Size10;
  61
  62         //mrcImage Conv;
  63
  64
  65         init0(&info);
  66         argCheck(&info, argc, argv);
  67         init1(&info);
  68
  69         // ファイル読み込み
  70         DEBUGPRINT("Program Start\n");
  71         mrcFileRead(&in1, info.In1, "in main", 0);
  72         mrcFileRead(&in2, info.In2, "in main", 0);
  73         if(info.flagIn3) {
  74                 mrcFileRead(&in3, info.In3, "in main", 0);
  75         }
  76
  77         // ファイルの大きさを揃えてる
  78         outIx.Header = in1.Header;
  79         outIy.Header = in1.Header;
  80         outIz.Header = in1.Header;
  81         outIt.Header = in1.Header;
  82     outU.Header  = in1.Header;
  83     outV.Header  = in1.Header;
  84     outW.Header  = in1.Header;
  85     Size.Header  = in1.Header;
  86         Size1.Header  = in1.Header;
  87         Size2.Header  = in1.Header;
  88         Size3.Header  = in1.Header;
  89         Size4.Header  = in1.Header;
  90         Size5.Header  = in1.Header;
  91         Size6.Header  = in1.Header;
  92         Size7.Header  = in1.Header;
  93         Size8.Header  = in1.Header;
  94         Size9.Header  = in1.Header;
  95         Size10.Header  = in1.Header;
  96
  97         //Conv.Header  = in1.Header;
  98
  99         //初期化
 100     mrcInit(&outIx, NULL);
 101     mrcInit(&outIy, NULL);
 102     mrcInit(&outIz, NULL);
 103     mrcInit(&outIt, NULL);
 104     mrcInit(&outU, NULL);
 105     mrcInit(&outV, NULL);
 106     mrcInit(&outW, NULL);
 107     mrcInit(&Size, NULL);
 108         mrcInit(&Size1, NULL);
 109         mrcInit(&Size2, NULL);
 110         mrcInit(&Size3, NULL);
 111         mrcInit(&Size4, NULL);
 112         mrcInit(&Size5, NULL);
 113         mrcInit(&Size6, NULL);
 114         mrcInit(&Size7, NULL);
 115         mrcInit(&Size8, NULL);
 116         mrcInit(&Size9, NULL);
 117         mrcInit(&Size10, NULL);
 118
 119         //mrcInit(&Conv, NULL);
 120
 121         //各画像の勾配のための変数
 122         double Ix, Iy, Iz, It;
 123         Ix = Iy = Iz = It = 0;
 124
 125         // Ix,Iy,Izの計算の一次データ保管場所としての変数
 126         double data1, data2, data3, data4, data5, data6, data7, data8, data9, data10, data11, data12, data13, data14, data15, data16;
 127
 128         data1 = data2 = data3 = data4 = data5 = data6 = data7 = data8 = data9 = data10 = data11 = data12 = data13 = data14 = data15 = data16 = 0;
 129
 130         //u,v,wの平均値の格納場所
 131         double u_ave, v_ave, w_ave;
 132
 133         //オプティカルフローの一回前の値の一時的な格納場所
 134         //double size_before = 0;
 135
 136         //収束を見るためのもの
 137         //double conv = 0;
 138
 139         //最終的に求めたいもの
 140         double u, v, w;
 141         u = v = w = 0;
 142
 143         //Horn-Schunck法の3D拡張
 144         int l;
 145         float x, y, z;
 146         double size;
 147
 148         // inputfileがある場合勾配の計算をする(Ix,Iy,Iz,Itの)
 149                 for (x = 1; x < in1.HeaderN.x - 1; x++) {
 150                         for (y = 1; y < in1.HeaderN.y - 1; y++) {
 151                                 for(z =1; z < in1.HeaderN.z - 1; z++) {
 152
 153                                 // オプテカルフローIx,Iy,Iz,Itの計算（平滑性制約に倣って）
 154                                 // 前フレームのIxの値取得
 155                                 mrcPixelDataGet(&in1, x  , y  , z  , &data1, mrcPixelRePart, mrcPixelHowNearest);
 156                                 mrcPixelDataGet(&in1, x+1, y  , z  , &data2, mrcPixelRePart, mrcPixelHowNearest);
 157                                 mrcPixelDataGet(&in1, x  , y+1, z  , &data3, mrcPixelRePart, mrcPixelHowNearest);
 158                                 mrcPixelDataGet(&in1, x+1, y+1, z  , &data4, mrcPixelRePart, mrcPixelHowNearest);
 159                                 mrcPixelDataGet(&in1, x  , y  , z+1, &data5, mrcPixelRePart, mrcPixelHowNearest);
 160                                 mrcPixelDataGet(&in1, x+1, y  , z+1, &data6, mrcPixelRePart, mrcPixelHowNearest);
 161                                 mrcPixelDataGet(&in1, x  , y+1, z+1, &data7, mrcPixelRePart, mrcPixelHowNearest);
 162                                 mrcPixelDataGet(&in1, x+1, y+1, z+1, &data8, mrcPixelRePart, mrcPixelHowNearest);
 163                                 // 後フレームの値取得
 164                                 mrcPixelDataGet(&in2, x  , y  , z  , &data9, mrcPixelRePart, mrcPixelHowNearest);
 165                 mrcPixelDataGet(&in2, x+1, y  , z  , &data10, mrcPixelRePart, mrcPixelHowNearest);
 166                 mrcPixelDataGet(&in2, x  , y+1, z  , &data11, mrcPixelRePart, mrcPixelHowNearest);
 167                 mrcPixelDataGet(&in2, x+1, y+1, z  , &data12, mrcPixelRePart, mrcPixelHowNearest);
 168                                 mrcPixelDataGet(&in2, x  , y  , z+1, &data13, mrcPixelRePart, mrcPixelHowNearest);
 169                                 mrcPixelDataGet(&in2, x+1, y  , z+1, &data14, mrcPixelRePart, mrcPixelHowNearest);
 170                                 mrcPixelDataGet(&in2, x  , y+1, z+1, &data15, mrcPixelRePart, mrcPixelHowNearest);
 171                                 mrcPixelDataGet(&in2, x+1, y+1, z+1, &data16, mrcPixelRePart, mrcPixelHowNearest);
 172                                 // Ixの計算
 173                                 Ix = 0.125 * ((data2-data1)+(data4-data3)+(data6-data5)+(data8-data7)+(data10-data9)+(data12-data11)+(data14-data13)+(data16-data15));
 174
 175                 // Iyの計算
 176                                 Iy = 0.125 * ((data3-data1)+(data4-data2)+(data7-data5)+(data8-data6)+(data11-data9)+(data12-data10)+(data15-data13)+(data16-data14));
 177
 178                 // Izの計算
 179                                 Iz = 0.125 * ((data5-data1)+(data6-data2)+(data7-data3)+(data8-data4)+(data13-data9)+(data14-data10)+(data15-data11)+(data16-data12));
 180
 181
 182                                 // 前フレームのItの計算
 183                 mrcPixelDataGet(&in1, x  , y  , z  , &data1, mrcPixelRePart, mrcPixelHowNearest);
 184                 mrcPixelDataGet(&in2, x  , y  , z  , &data2, mrcPixelRePart, mrcPixelHowNearest);
 185
 186                                 mrcPixelDataGet(&in1, x+1, y  , z  , &data3, mrcPixelRePart, mrcPixelHowNearest);
 187                 mrcPixelDataGet(&in2, x+1, y  , z  , &data4, mrcPixelRePart, mrcPixelHowNearest);
 188
 189                                 mrcPixelDataGet(&in1, x  , y+1, z  , &data5, mrcPixelRePart, mrcPixelHowNearest);
 190                 mrcPixelDataGet(&in2, x  , y+1, z  , &data6, mrcPixelRePart, mrcPixelHowNearest);
 191
 192                                 mrcPixelDataGet(&in1, x+1, y+1, z  , &data7, mrcPixelRePart, mrcPixelHowNearest);
 193                 mrcPixelDataGet(&in2, x+1, y+1, z  , &data8, mrcPixelRePart, mrcPixelHowNearest);
 194
 195
 196                 mrcPixelDataGet(&in1, x  , y  , z+1, &data9, mrcPixelRePart, mrcPixelHowNearest);
 197                 mrcPixelDataGet(&in2, x  , y  , z+1, &data10, mrcPixelRePart, mrcPixelHowNearest);
 198
 199                                 mrcPixelDataGet(&in1, x+1, y  , z+1, &data11, mrcPixelRePart, mrcPixelHowNearest);
 200                 mrcPixelDataGet(&in2, x+1, y  , z+1, &data12, mrcPixelRePart, mrcPixelHowNearest);
 201
 202                                 mrcPixelDataGet(&in1, x  , y+1, z+1, &data13, mrcPixelRePart, mrcPixelHowNearest);
 203                 mrcPixelDataGet(&in2, x  , y+1, z+1, &data14, mrcPixelRePart, mrcPixelHowNearest);
 204
 205                                 mrcPixelDataGet(&in1, x+1, y+1, z+1, &data15, mrcPixelRePart, mrcPixelHowNearest);
 206                 mrcPixelDataGet(&in2, x+1, y+1, z+1, &data16, mrcPixelRePart, mrcPixelHowNearest);
 207
 208                                 // Itの計算
 209                                 It = 0.125 * ((data2-data1)+(data4-data3)+(data6-data5)+(data8-data7)+(data10-data9)+(data12-data11)+(data14-data13)+(data16-data15));
 210
 211
 212                                 //DEBUGPRINT("Program Start %f", &x);
 213                                 mrcPixelDataSet(&outIx, x, y, z, Ix, mrcPixelRePart);
 214                                 mrcPixelDataSet(&outIy, x, y, z, Iy, mrcPixelRePart);
 215                                 mrcPixelDataSet(&outIz, x, y, z, Iz, mrcPixelRePart);
 216                                 mrcPixelDataSet(&outIt, x, y, z, It, mrcPixelRePart);
 217                                 }
 218                         }
 219                 }
 220
 221
 222         //outputfileがオプション指定されている場合、情報を書き込む
 223         if(info.flagOutIx) {
 224                 mrcFileWrite(&outIx, info.OutIx, "in main", 0);
 225                 mrcFileWrite(&outIy, info.OutIy, "in main", 0);
 226                 mrcFileWrite(&outIz, info.OutIz, "in main", 0);
 227                 mrcFileWrite(&outIt, info.OutIt, "in main", 0);
 228         }
 229
 230
 231
 232         // MAX_ITERATIONS従って反復して計算を行う
 233         for(l = 1; l <= info.MAX_ITERATIONS; l++) {
 234                 for(x = 1; x < in1.HeaderN.x - 1; x++){
 235                         for(y = 1; y < in1.HeaderN.y - 1; y++){
 236                                 for(z = 1; z < in1.HeaderN.z - 1; z++){
 237
 238                                         // オプティカルフローの更新を行う
 239                                         mrcPixelDataGet(&outIx, x, y, z, &Ix, mrcPixelRePart, mrcPixelHowNearest);
 240                                         mrcPixelDataGet(&outIx, x, y, z, &Iy, mrcPixelRePart, mrcPixelHowNearest);
 241                                         mrcPixelDataGet(&outIz, x, y, z, &Iz, mrcPixelRePart, mrcPixelHowNearest);
 242                                         mrcPixelDataGet(&outIt, x, y, z, &It, mrcPixelRePart, mrcPixelHowNearest);
 243
 244                                         //一回目の反復におけるu,v,w,の平均値の計算
 245                                         mrcPixelDataGet(&outU, x+1, y  , z  , &data1, mrcPixelRePart, mrcPixelHowNearest);
 246                         mrcPixelDataGet(&outU, x-1, y  , z  , &data2, mrcPixelRePart, mrcPixelHowNearest);
 247                         mrcPixelDataGet(&outU, x  , y+1, z  , &data3, mrcPixelRePart, mrcPixelHowNearest);
 248                         mrcPixelDataGet(&outU, x  , y-1, z  , &data4, mrcPixelRePart, mrcPixelHowNearest);
 249                                         mrcPixelDataGet(&outU, x  , y  , z+1, &data5, mrcPixelRePart, mrcPixelHowNearest);
 250                     mrcPixelDataGet(&outU, x  , y  , z-1, &data6, mrcPixelRePart, mrcPixelHowNearest);
 251
 252                                         u_ave = (data1 + data2 + data3 + data4 + data5 + data6) / 6;
 253
 254                                         mrcPixelDataGet(&outV, x+1, y  , z  , &data1, mrcPixelRePart, mrcPixelHowNearest);
 255                     mrcPixelDataGet(&outV, x-1, y  , z  , &data2, mrcPixelRePart, mrcPixelHowNearest);
 256                     mrcPixelDataGet(&outV, x  , y+1, z  , &data3, mrcPixelRePart, mrcPixelHowNearest);
 257                     mrcPixelDataGet(&outV, x  , y-1, z  , &data4, mrcPixelRePart, mrcPixelHowNearest);
 258                     mrcPixelDataGet(&outV, x  , y  , z+1, &data5, mrcPixelRePart, mrcPixelHowNearest);
 259                     mrcPixelDataGet(&outV, x  , y  , z-1, &data6, mrcPixelRePart, mrcPixelHowNearest);
 260
 261                     v_ave = (data1 + data2 + data3 + data4 + data5 + data6) / 6;
 262
 263                     mrcPixelDataGet(&outW, x+1, y  , z  , &data1, mrcPixelRePart, mrcPixelHowNearest);
 264                     mrcPixelDataGet(&outW, x-1, y  , z  , &data2, mrcPixelRePart, mrcPixelHowNearest);
 265                     mrcPixelDataGet(&outW, x  , y+1, z  , &data3, mrcPixelRePart, mrcPixelHowNearest);
 266                     mrcPixelDataGet(&outW, x  , y-1, z  , &data4, mrcPixelRePart, mrcPixelHowNearest);
 267                     mrcPixelDataGet(&outW, x  , y  , z+1, &data5, mrcPixelRePart, mrcPixelHowNearest);
 268                     mrcPixelDataGet(&outW, x  , y  , z-1, &data6, mrcPixelRePart, mrcPixelHowNearest);
 269
 270                     w_ave = (data1 + data2 + data3 + data4 + data5 + data6) / 6;
 271
 272                                         //更新されたu,v,w の計算
 273                                         //u = u_ave - (Ix * (Ix * u_ave + Iy * v_ave + Iz * w_ave + It)) / ((1 / info.ALPHA) + Ix * Ix + Iy * Iy + Iz * Iz);
 274                                         //v = v_ave - (Iy * (Ix * u_ave + Iy * v_ave + Iz * w_ave + It)) / ((1 / info.ALPHA) + Ix * Ix + Iy * Iy + Iz * Iz);
 275                                         //w = w_ave - (Iz * (Ix * u_ave + Iy * v_ave + Iz * w_ave + It)) / ((1 / info.ALPHA) + Ix * Ix + Iy * Iy + Iz * Iz);
 276
 277                                         u = u_ave - (Ix * (Ix * u_ave + Iy * v_ave + Iz * w_ave + It)) / (info.ALPHA * info.ALPHA + Ix * Ix + Iy * Iy + Iz * Iz);
 278                                         v = v_ave - (Iy * (Ix * u_ave + Iy * v_ave + Iz * w_ave + It)) / (info.ALPHA * info.ALPHA + Ix * Ix + Iy * Iy + Iz * Iz);
 279                                         w = w_ave - (Iz * (Ix * u_ave + Iy * v_ave + Iz * w_ave + It)) / (info.ALPHA * info.ALPHA + Ix * Ix + Iy * Iy + Iz * Iz);
 280
 281                                         size = sqrt(u*u + v*v + w*w);
 282
 283                                         //オプティカルフロー の差を出してどのぐらい収束するのか知るために...
 284                                         if(info.flagSize1) {
 285                                                 if(l == 0){
 286                                                         mrcPixelDataSet(&Size1, x, y, z, size, mrcPixelRePart);
 287                                                 }
 288                                                 if(l == 1){
 289                                 mrcPixelDataSet(&Size2, x, y, z, size, mrcPixelRePart);
 290                         }
 291                         if(l == 2){
 292                                 mrcPixelDataSet(&Size3, x, y, z, size, mrcPixelRePart);
 293                         }
 294                         if(l == 3){
 295                                 mrcPixelDataSet(&Size4, x, y, z, size, mrcPixelRePart);
 296                         }
 297                         if(l == 9){
 298                                 mrcPixelDataSet(&Size5, x, y, z, size, mrcPixelRePart);
 299                         }
 300                         if(l == 19){
 301                                 mrcPixelDataSet(&Size6, x, y, z, size, mrcPixelRePart);
 302                         }
 303                         if(l == 29){
 304                                 mrcPixelDataSet(&Size7, x, y, z, size, mrcPixelRePart);
 305                         }
 306                         if(l == 39){
 307                                 mrcPixelDataSet(&Size8, x, y, z, size, mrcPixelRePart);
 308                         }
 309                         if(l == 49){
 310                                 mrcPixelDataSet(&Size9, x, y, z, size, mrcPixelRePart);
 311                         }
 312                         if(l == 99){
 313                                 mrcPixelDataSet(&Size10, x, y, z, size, mrcPixelRePart);
 314                         }
 315                                         }
 316                                         mrcPixelDataSet(&outU, x, y, z, u,    mrcPixelRePart);
 317                                         mrcPixelDataSet(&outV, x, y, z, v,    mrcPixelRePart);
 318                                         mrcPixelDataSet(&outW, x, y, z, w,    mrcPixelRePart);
 319                                         mrcPixelDataSet(&Size, x, y, z, size, mrcPixelRePart);
 320                                         //printf("反復数%d回目\n", l);
 321                                 }
 322                         }
 323                 }
 324                 printf("反復数%d回目\n", l);
 325         }
 326
 327
 328         if(info.flagoutU) {
 329                 mrcFileWrite(&outU, info.outU, "in main", 0);
 330                 mrcFileWrite(&outV, info.outV, "in main", 0);
 331                 mrcFileWrite(&outW, info.outW, "in main", 0);
 332         }
 333         if(info.flagSize) {
 334                 mrcFileWrite(&Size, info.Size, "in main", 0);
 335         }
 336     if(info.flagSize1) {
 337         mrcFileWrite(&Size1, info.Size1, "in main", 0);
 338                 mrcFileWrite(&Size2, info.Size2, "in main", 0);
 339                 mrcFileWrite(&Size3, info.Size3, "in main", 0);
 340                 mrcFileWrite(&Size4, info.Size4, "in main", 0);
 341                 mrcFileWrite(&Size5, info.Size5, "in main", 0);
 342                 mrcFileWrite(&Size6, info.Size6, "in main", 0);
 343                 mrcFileWrite(&Size7, info.Size7, "in main", 0);
 344                 mrcFileWrite(&Size8, info.Size8, "in main", 0);
 345                 mrcFileWrite(&Size9, info.Size9, "in main", 0);
 346                 mrcFileWrite(&Size10, info.Size10, "in main", 0);
 347         }
 348
 349         /*
 350         if(info.flagConv) {
 351         mrcFileWrite(&Conv, info.Conv, "in main", 0);
 352     }
 353         */
 354
 355         return 0;
 356 }
 357
 358
 359 void
 360 additionalUsage()
 361 {
 362         fprintf(stderr, "----- Additional Usage -----\n");
 363 }