src/Tools/rec3d/mrcImageOrientationSearch/src/mrcImageOrientationSearch.c

   1 /*
   2 # mrcImageOrientationSearch : $Revision$
   3 # $Date$
   4 # Created by $Author$
   5 # Usage : mrcImageOrientationSearch
   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 "Memory.h"
  22 #include "eosString.h"
  23 #include "lCommonLineCalculation.h"
  24
  25 #define WORDLEN 1024
  26
  27 typedef struct SinogramPositionInfo {
  28         int     PositionNum;
  29         char*   filename1;
  30         char*   filename2;
  31         /* DEGREE */
  32         double* CommonLinePosition1;
  33         double* CommonLinePosition2;
  34         int*    mode;
  35 } SinogramPositionInfo;
  36
  37 typedef struct EulerAngleInfo {
  38         char*  filename;
  39         int    flag;
  40         /* DEGREE */
  41         double theta;
  42         double phi;
  43         double psi;
  44 } EulerAngleInfo;
  45
  46 typedef struct lmrcImageOrientationSearchInfo {
  47         char** CorrelationFilename;
  48         int    Correlationfilenum;
  49         int    filenum;
  50
  51         SinogramPositionInfo* position;
  52         EulerAngleInfo**      angle;
  53         char*                 rotationalMode;
  54 } lmrcImageOrientationSearchInfo;
  55
  56 /*
  57 typedef enum lmrcImageOrientationSearchMode {
  58         a=0,
  59         b=1
  60 } lmrcImageOrientationSearchMode;
  61 */
  62
  63 void
  64 __FileListReadModeSelect(lmrcImageOrientationSearchInfo* linfo, FILE* fpt, int mode);
  65
  66 void
  67 __FileListAndCommonLinePositionRead(lmrcImageOrientationSearchInfo* linfo, FILE* fpt, int mode);
  68
  69 void
  70 __FileListRead(lmrcImageOrientationSearchInfo* linfo, FILE* fpt, int mode);
  71
  72 void
  73 lmrcImageOrientationSearch(lmrcImageOrientationSearchInfo* linfo, int mode);
  74
  75 int
  76 lmrcImageOrientationSearchOut(FILE* fpt, lmrcImageOrientationSearchInfo* linfo, int mode);
  77
  78 void
  79 lEulerAngleCalculationFromThreeCommonLinePosition(lmrcImageOrientationSearchInfo* linfo, int number, int* flag, int mode);
  80
  81 void
  82 lEulerAngleCalculationFromAdditionlCommonLinePosition(lmrcImageOrientationSearchInfo* linfo, int num, int mode);
  83
  84
  85 int
  86 main(int argc, char* argv[])
  87 {
  88         mrcImageOrientationSearchInfo  info;
  89         lmrcImageOrientationSearchInfo linfo;
  90         int i;
  91
  92         init0(&info);
  93     argCheck(&info, argc, argv);
  94     init1(&info);
  95
  96         DEBUGPRINT("Program Start\n");
  97
  98         /* file and correlation file number copy in linfo*/
  99         linfo.Correlationfilenum = info.flagIn1;
 100         linfo.filenum            = info.flagIn2;
 101         //linfo.Correlationfilenum = linfo.filenum*(linfo.filenum-1)/2;
 102
 103         /* input file2 */
 104         linfo.angle    = MemoryAllocate(EulerAngleInfo*,            2, "in main");
 105         linfo.angle[0] = MemoryAllocate(EulerAngleInfo, linfo.filenum, "in main");
 106         linfo.angle[1] = MemoryAllocate(EulerAngleInfo, linfo.filenum, "in main");
 107         for(i=0; i<linfo.filenum; i++){
 108                 linfo.angle[0][i].filename = MemoryAllocate(char, WORDLEN, "in main");
 109                 linfo.angle[1][i].filename = MemoryAllocate(char, WORDLEN, "in main");
 110                 stringCopy(linfo.angle[0][i].filename, info.In2[i], WORDLEN-1);
 111                 stringCopy(linfo.angle[1][i].filename, info.In2[i], WORDLEN-1);
 112                 linfo.angle[0][i].flag = 1;
 113                 linfo.angle[1][i].flag = 1;
 114         }
 115
 116         switch(info.mode){
 117                 case 0 :{
 118                         break;
 119                 }
 120                 case 1 :
 121                 case 2 :{
 122                         linfo.CorrelationFilename = MemoryAllocate(char*, linfo.Correlationfilenum, "in main");
 123                         for(i=0; i<linfo.Correlationfilenum; i++){
 124                                 linfo.CorrelationFilename[i] = MemoryAllocate(char, WORDLEN, "in main");
 125                                 strcpy(linfo.CorrelationFilename[i], info.In1[i]);
 126                         }
 127                         break;
 128                 }
 129         }
 130
 131         /* input file1 */
 132         __FileListReadModeSelect(&linfo, info.fptIn1List, info.mode);
 133
 134         /* determine Euler angle */
 135         lmrcImageOrientationSearch(&linfo, info.mode);
 136
 137         /* Output */
 138         lmrcImageOrientationSearchOut(info.fptOut, &linfo, info.mode);
 139         DEBUGPRINT("Program End\n");
 140
 141         return EXIT_SUCCESS;
 142 }
 143
 144 void
 145 additionalUsage()
 146 {
 147         fprintf(stderr, "\n");
 148         fprintf(stderr, "----- mode -----\n");
 149         fprintf(stderr, "-m : 0 : There is one common line position in input file (using Test)\n");
 150         fprintf(stderr, "     1 : There is    correlation file name in input file (common line position is determined Max correlation)\n");
 151         fprintf(stderr, "     2 : There is    correlation file name in input file (common line position is determined Threshold), -T\n");
 152         fprintf(stderr, "\n");
 153         fprintf(stderr, "----- input file 1 format -----\n");
 154         fprintf(stderr, "-m : 0 : filename1 filename2 commonlineposition1 commonlineposition2 commonlinepositionMode\n");
 155         fprintf(stderr, "             .         .             .                   .                     .           \n");
 156         fprintf(stderr, "             .         .             .                   .                     .           \n");
 157         fprintf(stderr, "             .         .             .                   .                     .           \n");
 158         fprintf(stderr, "             .         .             .                   .                     .           \n");
 159         fprintf(stderr, "\n");
 160         fprintf(stderr, "     1 : correlationfilename filename1 filename2\n");
 161         fprintf(stderr, "     2 :         .               .         .    \n");
 162         fprintf(stderr, "                 .               .         .    \n");
 163         fprintf(stderr, "                 .               .         .    \n");
 164         fprintf(stderr, "                 .               .         .    \n");
 165         fprintf(stderr, "\n");
 166         fprintf(stderr, "----- input file 2 format -----\n");
 167         fprintf(stderr, "-m : 0 : filename \n");
 168         fprintf(stderr, "     1 :    .     \n");
 169         fprintf(stderr, "     2 :    .     \n");
 170         fprintf(stderr, "            .     \n");
 171         fprintf(stderr, "            .     \n");
 172         fprintf(stderr, "\n");
 173 }
 174
 175 void
 176 __FileListReadModeSelect(lmrcImageOrientationSearchInfo* linfo, FILE* fpt, int mode)
 177 {
 178         switch(mode){
 179                 case 0 :{
 180                         __FileListAndCommonLinePositionRead(linfo, fpt, 0);
 181                         break;
 182                 }
 183                 case 1 :
 184                 case 2 :{
 185                         __FileListRead(linfo, fpt, 0);
 186                         break;
 187                 }
 188         }
 189 }
 190
 191 void
 192 __FileListAndCommonLinePositionRead(lmrcImageOrientationSearchInfo* linfo, FILE* fpt, int mode)
 193 {
 194         int   i;
 195         char* s;
 196         char* m;
 197
 198         s = MemoryAllocate(char, WORDLEN, "in __FileListAndCommonLinePositionRead");
 199         m = MemoryAllocate(char, WORDLEN, "in __FileListAndCommonLinePositionRead");
 200
 201         linfo->position = MemoryAllocate(SinogramPositionInfo, linfo->Correlationfilenum, "in __FileListAndCommonLinePositionRead");
 202
 203         fseek(fpt, 0L, SEEK_SET);
 204         for(i=0; i<linfo->Correlationfilenum; i++){
 205                 linfo->position[i].PositionNum         = 1;
 206                 linfo->position[i].filename1           = MemoryAllocate(char  , WORDLEN, "in __FileListAndCommonLinePositionRead");
 207                 linfo->position[i].filename2           = MemoryAllocate(char  , WORDLEN, "in __FileListAndCommonLinePositionRead");
 208                 linfo->position[i].CommonLinePosition1 = MemoryAllocate(double, 1      , "in __FileListAndCommonLinePositionRead");
 209                 linfo->position[i].CommonLinePosition2 = MemoryAllocate(double, 1      , "in __FileListAndCommonLinePositionRead");
 210                 linfo->position[i].mode                = MemoryAllocate(int   , 1      , "in __FileListAndCommonLinePositionRead");
 211                 stringGetFromFile(s, "", fpt, stdout, 3);
 212                 stringCopy(linfo->position[i].filename1, stringGetNthWord(s, 2, " \t,"), WORDLEN-1);
 213                 stringCopy(linfo->position[i].filename2, stringGetNthWord(s, 3, " \t,"), WORDLEN-1);
 214                 linfo->position[i].CommonLinePosition1[0] = stringGetNthRealData(s, 4, " \t,");
 215                 linfo->position[i].CommonLinePosition2[0] = stringGetNthRealData(s, 5, " \t,");
 216                 stringCopy(m, stringGetNthWord(s, 6, " \t,"), WORDLEN-1);
 217                 if(strcoll(m, "1") == 0){
 218                         linfo->position[i].mode = 0;
 219                 }
 220                 else if(strcoll(m, "z-rotation") == 0){
 221                         linfo->position[i].mode[0] = 1;
 222                 }
 223                 else{
 224                         linfo->position[i].mode[0] = 2;
 225                         DEBUGPRINT1("Is not supported Common line position mode : in input file line >> %d\n", i);
 226                 }
 227         }
 228         free(m);
 229         free(s);
 230 }
 231
 232 void
 233 __FileListRead(lmrcImageOrientationSearchInfo* linfo, FILE* fpt, int mode)
 234 {
 235         int   i;
 236         char* s;
 237
 238         s = MemoryAllocate(char, WORDLEN, "in __FileListRead");
 239
 240         linfo->position = MemoryAllocate(SinogramPositionInfo, linfo->Correlationfilenum, "in __FileListRead");
 241
 242         fseek(fpt, 0L, SEEK_SET);
 243         for(i=0; i<linfo->Correlationfilenum; i++){
 244                 linfo->position[i].filename1 = MemoryAllocate(char, WORDLEN, "in __FileListRead");
 245                 linfo->position[i].filename2 = MemoryAllocate(char, WORDLEN, "in __FileListRead");
 246                 stringGetFromFile(s, "", fpt, stdout, 3);
 247                 stringCopy(linfo->position[i].filename1, stringGetNthWord(s, 2, " \t,"), WORDLEN-1);
 248                 stringCopy(linfo->position[i].filename2, stringGetNthWord(s, 3, " \t,"), WORDLEN-1);
 249         }
 250         free(s);
 251 }
 252
 253
 254 int
 255 lmrcImageOrientationSearchOut(FILE* fpt, lmrcImageOrientationSearchInfo* linfo, int mode)
 256 {
 257         int i, j;
 258         int flag;
 259
 260         flag = 1;
 261         i    = 2;
 262
 263         switch(mode){
 264                 case 0 :{
 265                         for(i=0; i<2; i++){
 266                         for(j=0; j<linfo->filenum; j++){
 267                                 if(linfo->angle[i][j].flag == 0){
 268                                         fprintf(fpt, "%s %s %f %f %f\n", linfo->angle[i][j].filename, linfo->rotationalMode,
 269                                                 linfo->angle[i][j].psi, linfo->angle[i][j].theta, linfo->angle[i][j].phi);
 270                                 }
 271                         }
 272                         }
 273                         break;
 274                 }
 275                 case 1 :{
 276                         break;
 277                 }
 278                 case 2 :{
 279                         break;
 280                 }
 281                 default: {
 282                         break;
 283                 }
 284         }
 285 }
 286
 287 void
 288 lmrcImageOrientationSearch(lmrcImageOrientationSearchInfo* linfo, int mode)
 289 {
 290         int i, j;
 291         int flag;
 292
 293         flag = 1;
 294         i    = 2;
 295         switch(mode){
 296                 case 0 :{
 297                         while(flag == 1 && i < linfo->filenum){
 298                                 lEulerAngleCalculationFromThreeCommonLinePosition(linfo, i, &flag, 0);
 299                                 i++;
 300                         }
 301                         DEBUGPRINT("\n");
 302                         for(i=0; i<2; i++){
 303                         for(j=0; j<linfo->filenum; j++){
 304                                 if(linfo->angle[i][j].flag == 0){
 305                                         DEBUGPRINT5("filename : %s >>> %s (%f, %f, %f)\n", linfo->angle[i][j].filename, linfo->rotationalMode,
 306                                                 linfo->angle[i][j].psi, linfo->angle[i][j].theta, linfo->angle[i][j].phi);
 307                                 }
 308                         }
 309                         DEBUGPRINT("\n");
 310                         }
 311                         //DEBUGPRINT("program end\n");
 312                         //exit(1);
 313                         for(j=i; j<linfo->Correlationfilenum; j++){
 314                                 lEulerAngleCalculationFromAdditionlCommonLinePosition(linfo, j, 0);
 315                         }
 316                         break;
 317                 }
 318                 case 1 :{
 319                         break;
 320                 }
 321                 case 2 :{
 322                         break;
 323                 }
 324         }
 325 }
 326
 327 void
 328 lEulerAngleCalculationFromThreeCommonLinePosition(lmrcImageOrientationSearchInfo* linfo, int number, int* flag, int mode)
 329 {
 330         int*   wnum;
 331         int    i;
 332         double* angle;
 333         double* cosine;
 334         double* sine;
 335         double* Xrotation1;
 336         double  cosineXr1;
 337         double* Xrotation2;
 338         double  cosineXr2;
 339         double  ratio;
 340
 341         wnum = MemoryAllocate(int , 3, "lEulerAngleCalculationFromThreeCommonLinePosition");
 342         for(i=0; i<linfo->Correlationfilenum; i++){
 343                 if(strcoll(linfo->angle[0][0].filename, linfo->position[i].filename1) == 0 &&
 344                    strcoll(linfo->angle[0][1].filename, linfo->position[i].filename2) == 0)
 345                 {
 346                         wnum[0] = i;
 347                 }
 348                 if(strcoll(linfo->angle[0][0].filename, linfo->position[i].filename1) == 0 &&
 349                    strcoll(linfo->angle[0][number].filename, linfo->position[i].filename2) == 0)
 350                 {
 351                         wnum[1] = i;
 352                 }
 353                 if(strcoll(linfo->angle[0][1].filename, linfo->position[i].filename1) == 0 &&
 354                    strcoll(linfo->angle[0][number].filename, linfo->position[i].filename2) == 0)
 355                 {
 356                         wnum[2] = i;
 357                 }
 358         }
 359         angle = MemoryAllocate(double, 6, "lEulerAngleCalculationFromThreeCommonLinePosition");
 360         for(i=0 ;i<3; i++){
 361                 angle[i*2]   = linfo->position[wnum[i]].CommonLinePosition1[0]*RADIAN;
 362                 angle[i*2+1] = linfo->position[wnum[i]].CommonLinePosition2[0]*RADIAN;
 363         }
 364         cosine = MemoryAllocate(double, 6, "lEulerAngleCalculationFromThreeCommonLinePosition");
 365         sine   = MemoryAllocate(double, 6, "lEulerAngleCalculationFromThreeCommonLinePosition");
 366         for(i=0; i<6; i++){
 367                 cosine[i] = cos(angle[i]);
 368                 sine[i]   = sin(angle[i]);
 369                 DEBUGPRINT4("angle[%d] = %f >> cosine = %f sine = %f\n", i, angle[i], cosine[i], sine[i]);
 370         }
 371
 372         cosineXr1 = 0;
 373         cosineXr2 = 0;
 374
 375         cosineXr1 = ((cosine[4]*cosine[5]+sine[4]*sine[5])-((cosine[2]*cosine[3]+sine[2]*sine[3])*(cosine[0]*cosine[1]+sine[0]*sine[1])))
 376                    /((cosine[2]*sine[3]-sine[2]*cosine[3])*(cosine[0]*sine[1]-sine[0]*cosine[1]));
 377
 378         cosineXr2 = (((cosine[2]*sine[3]-sine[2]*cosine[3])*(cosine[0]*cosine[1]+sine[0]*sine[1])*cosineXr1)
 379                     +((cosine[2]*cosine[3]+sine[2]*sine[3])*(sine[0]*cosine[1]-cosine[0]*sine[1])))
 380                    /(cosine[4]*sine[5]-sine[4]*cosine[5]);
 381
 382         DEBUGPRINT2("cosineXr1 = %f cosineXr2 = %f \n", cosineXr1, cosineXr2);
 383
 384         ratio = (cosine[4]*sine[5]-sine[4]*cosine[5])/(cosine[2]*sine[3]-sine[2]*cosine[3]);
 385
 386         DEBUGPRINT1("ratio = %f \n", ratio);
 387
 388         Xrotation1 = MemoryAllocate(double, 2, "in lEulerAngleCalculationFromThreeCommonLinePosition");
 389         Xrotation2 = MemoryAllocate(double, 2, "in lEulerAngleCalculationFromThreeCommonLinePosition");
 390         if(ratio < -0.0000001 || ratio > 0.0000001){
 391                 Xrotation1[0] = acos(cosineXr1);
 392                 Xrotation1[1] = 2*M_PI-Xrotation1[0];
 393                 Xrotation2[0] = acos(cosineXr2);
 394                 Xrotation2[1] = 2*M_PI-Xrotation2[0];
 395                 if(ratio*sin(Xrotation1[0] < 0)){
 396                         Xrotation2[0] = 2*M_PI-Xrotation2[0];
 397                         Xrotation2[1] = Xrotation2[0];
 398                 }
 399                 DEBUGPRINT2("Xrotation1 = %f Xrotation2 = %f\n",Xrotation1[0], Xrotation2[0]);
 400                 linfo->rotationalMode = MemoryAllocate(char, 5, "in lEulerAngleCalculationFromThreeCommonLinePosition");
 401                 stringCopy(linfo->rotationalMode, "ZENR", 4); /* Z->X->Z rotation */
 402                 linfo->angle[0][0].theta   = 0;
 403                 linfo->angle[0][0].phi     = 0;
 404                 linfo->angle[0][0].psi     = 0;
 405                 linfo->angle[0][1].theta   = linfo->position[wnum[0]].CommonLinePosition1[0];
 406                 linfo->angle[0][1].phi     = Xrotation1[0]*DEGREE;
 407                 linfo->angle[0][1].psi     = linfo->position[wnum[0]].CommonLinePosition2[0];
 408                 linfo->angle[0][number].theta = linfo->position[wnum[1]].CommonLinePosition1[0];
 409                 linfo->angle[0][number].phi   = Xrotation2[0]*DEGREE;
 410                 linfo->angle[0][number].psi   = linfo->position[wnum[1]].CommonLinePosition2[0];
 411
 412                 linfo->angle[1][0].theta   = 0;
 413                 linfo->angle[1][0].phi     = 0;
 414                 linfo->angle[1][0].psi     = 0;
 415                 linfo->angle[1][1].theta   = linfo->position[wnum[0]].CommonLinePosition1[0];
 416                 linfo->angle[1][1].phi     = Xrotation1[1]*DEGREE;
 417                 linfo->angle[1][1].psi     = linfo->position[wnum[0]].CommonLinePosition2[0];
 418                 linfo->angle[1][number].theta = linfo->position[wnum[1]].CommonLinePosition1[0];
 419                 linfo->angle[1][number].phi   = Xrotation2[1]*DEGREE;
 420                 linfo->angle[1][number].psi   = linfo->position[wnum[1]].CommonLinePosition2[0];
 421
 422                 /* angle is determined */
 423                 *flag = 0;
 424                 linfo->angle[0][0].flag = 0;
 425                 linfo->angle[0][1].flag = 0;
 426                 linfo->angle[0][number].flag = 0;
 427                 linfo->angle[1][0].flag = 0;
 428                 linfo->angle[1][1].flag = 0;
 429                 linfo->angle[1][number].flag = 0;
 430         } else {
 431                 DEBUGPRINT3("Can't determine Euler Angle : %s %s %s\n", linfo->angle[0][0].filename, linfo->angle[0][1].filename,
 432                                                                        linfo->angle[0][number].filename);
 433                 /* angle is not determined */
 434                 *flag = 1;
 435                 linfo->angle[0][0].flag = 1;
 436                 linfo->angle[0][1].flag = 1;
 437                 linfo->angle[0][number].flag = 1;
 438                 linfo->angle[1][0].flag = 1;
 439                 linfo->angle[1][1].flag = 1;
 440                 linfo->angle[1][number].flag = 1;
 441         }
 442
 443         free(wnum);
 444         free(angle);
 445         free(cosine);
 446         free(sine);
 447         free(Xrotation1);
 448         free(Xrotation2);
 449 }
 450
 451 void
 452 lEulerAngleCalculationFromAdditionlCommonLinePosition(lmrcImageOrientationSearchInfo* linfo, int num, int mode)
 453 {
 454
 455 }