src/Tools/rec3d/mrc2Dto3D/src/mrc2Dto3D.c

   1 /*
   2 # %M% %Y% %I%
   3 # The latest update : %G% at %U%
   4 #
   5 #%Z% mrc2Dto3D ver %I%
   6 #%Z% Created by
   7 #%Z%
   8 #%Z% Usage : mrc2Dto3D
   9 #%Z% Attention
  10 #%Z%
  11 */
  12 static char __sccs_id[] = "%Z%mrc2Dto3D ver%I%; Date:%D% %Z%";
  13 #include <stdio.h>
  14 #include <stdlib.h>
  15 #include <string.h>
  16 #include <math.h>
  17 #define GLOBAL_DECLARATION
  18 #undef DEBUG
  19 #include "../inc/config.h"
  20 #include "genUtil.h"
  21 #include "Memory.h"
  22 #include "eosString.h"
  23 #include "mrcImage.h"
  24 #include "Matrix3D.h"
  25 #include "eosPThread.h"
  26 #include "lmrcImageRhoFiltering.h"
  27 #include "lmrc2Dto3D.h"
  28
  29 extern void lmrc2Dto2DMulti(mrcImage* out, char** filename, int number, FILE* fptInfo, int mode);
  30 extern void lmrc2Dto2DEach(lmrc2Dto3DInfo* linfo, char** filename, int number, FILE* fptInfo, int mode);
  31
  32 int
  33 main(int argc, char* argv[])
  34 {
  35         mrc2Dto3DInfo info;
  36         static lmrc2Dto3DInfo linfo;
  37         static lmrc2Dto3DSIRTInfo llinfo;
  38         static mrcImage  In;
  39         static mrcImage Out;
  40
  41         init0(&info);
  42         argCheck(&info, argc, argv);
  43     init1(&info);
  44
  45         linfo.Rmax = info.Rmax;
  46         linfo.mode = info.mode;
  47         linfo.singleTiltMode = info.single;
  48         linfo.singleFilterMode = info.singleFilter;
  49         linfo.flagDouble     = info.Double;
  50         linfo.flagPlusXrot   = info.PlusXrot;
  51         linfo.InterpolationMode = info.InterpolationMode;
  52         linfo.rhoInfo.counterThreshold = info.CounterThreshold;
  53         linfo.rhoInfo.counterThresholdMode = info.CounterThresholdMode;
  54         linfo.rhoInfo.weightMode     = info.WeightMode;
  55         linfo.rhoInfo.SubSampling    = info.SubSampling;
  56         linfo.rhoInfo.flagThicknessWeight = info.thicknessWeight;
  57         __eosPThread__ = info.flagpthreadMax;
  58         __eosPThreadNum__ = info.pthreadMax;
  59
  60         if(info.flagNx) {
  61                 linfo.flagOutputSize = 1;
  62                 linfo.Nx = info.Nx;
  63                 linfo.Ny = info.Ny;
  64                 linfo.Nz = info.Nz;
  65         }
  66         if(info.SIRT) {
  67                 llinfo.maxIter = info.maxIter;
  68                 llinfo.rms = info.rms;
  69                 llinfo.lambda = 0.1;
  70         } else {
  71         }
  72
  73         if(info.flagIn) {
  74                 mrcFileRead(&In, info.In, "in main", 0);
  75         } else if (info.flagIn2) {
  76                 if(info.each) {
  77                         lmrc2Dto2DEach(&linfo, info.In2, info.flagIn2, info.fptIn2List, 0);
  78                 } else {
  79                         lmrc2Dto2DMulti(&In, info.In2, info.flagIn2, info.fptIn2List, 0);
  80                 }
  81         } else {
  82                 fprintf(stderr, "-i or -I is required\n");
  83                 exit(EXIT_FAILURE);
  84         }
  85
  86         if(info.SIRT) {
  87                 lmrc2Dto3DSIRT(&Out, &In, &linfo, &llinfo, 0);
  88         } else {
  89                 if(info.each) {
  90                         lmrc2Dto3D(&Out, NULL, &linfo, 0);
  91                 } else {
  92                         lmrc2Dto3D(&Out, &In, &linfo, 0);
  93                 }
  94         }
  95
  96         mrcFileWrite(&Out, info.Out, "in main", 0);
  97         if(info.flagOut2) {
  98                 mrcFileWrite(&In, info.Out2, "in main", 0);
  99         }
 100         if(info.Double && info.flagDoubleCounter) {
 101                 if(NULL!=linfo.CounterForWeight) {
 102                         mrcFileWrite(linfo.CounterForWeight, info.DoubleCounter, "in main", 0);
 103                 }
 104         }
 105         exit(EXIT_SUCCESS);
 106 }
 107
 108 void
 109 additionalUsage()
 110 {
 111         fprintf(stderr, "----- Additional Usage -----\n");
 112         fprintf(stderr, "-m Option\n");
 113         fprintf(stderr, "    %d:SimpleBackProjection\n", mrc2Dto3DModeSimpleBackProjection);
 114         fprintf(stderr, "    %d:FilteredBackProjection(Fourier Space)\n", mrc2Dto3DModeFilteredBackProjection);
 115         fprintf(stderr, "    %d:WeightedBackProjection(Real Space)\n", mrc2Dto3DModeWeightedBackProjection);
 116         fprintf(stderr, "-single 0|1 \n");
 117         fprintf(stderr, "    0: tilt axis is parallel to x-axis\n");
 118         fprintf(stderr, "    1: tilt axis is parallel to y-axis\n");
 119         fprintf(stderr, "-singleFilter 0|1|2 \n");
 120         fprintf(stderr, "    0: simple rho filter\n");
 121         fprintf(stderr, "    1: Ram-Lak Filter\n");
 122         fprintf(stderr, "    2: Shepp-Logan Filter\n");
 123         fprintf(stderr, "-Double \n");
 124         fprintf(stderr, "    Double Tilt \n");
 125         fprintf(stderr, "-WeightMode \n");
 126         fprintf(stderr, "    1 : RealSpace: Circle(same density)\n");
 127         fprintf(stderr, "    2 : RealSpace: Circle(weighted density)\n");
 128         fprintf(stderr, "    3 : RealSpace: Square(weighted density) \n");
 129         fprintf(stderr, "    4 : Fourier Space : Plane(same density) -CounterThreshold 0.5\n");
 130         fprintf(stderr, "    5 : Fourier Space : Plane(Linear Gradient)-CounterThreshold 0.5\n");
 131         fprintf(stderr, "    6 : Fourier Space : Plane(Cosine Gradient) Current Recommende using -CounterThreshold 0.5\n");
 132         fprintf(stderr, "-------------------------------\n");
 133         fprintf(stderr, "-I Option file format\n");
 134         fprintf(stderr, "filename0 RotationOrder0 rot1 rot2 rot3\n");
 135         fprintf(stderr, "filename1 RotationOrder1 rot1 rot2 rot3\n");
 136         fprintf(stderr, ".......................................\n");
 137         fprintf(stderr, ">>> RotationOrder : Eular Angle Expression <<< \n");
 138         fprintf(stderr, "Example: YOYS : RotY(rot3)RotX(rot2)RotY(rot1)*v\n");
 139         fprintf(stderr, "First Rotation  : y-axis : Y: Y     : [X|Y|Z] Axis used initially\n");
 140         fprintf(stderr, "Second Rotation : x-axis : O: Odd   : [O|E]   Parity of axis permutation\n");
 141         fprintf(stderr, "Last Rotation   : z-axis : Y: Yes   : [Y|N]   Repetition of initial axis as last\n");
 142         fprintf(stderr, "v1 = A v0                : S: Staic : [S|R]   Frame from which axes are taken\n");
 143
 144 }
 145
 146 void
 147 lmrc2Dto2DMulti(mrcImage* out, char** filename, int number, FILE* fptInfo, int mode)
 148 {
 149         long i;
 150         mrcImage* in;
 151         mrcImageParaTypeInteger Nx, Ny;
 152         mrcImageParaTypeReal    Length;
 153         mrcImageParaTypeReal    srcx, srcy, srcz;
 154         mrcImageParaTypeReal    dstx, dsty, dstz;
 155         char s[1024];
 156         double data;
 157
 158         in = (mrcImage*)memoryAllocate(sizeof(mrcImage)*number, "in lmrc2Dto2DMulti");
 159         Nx = 0;
 160         Ny = 0;
 161         Length = 1e6;
 162         for(i=0; i<number; i++) {
 163                 DEBUGPRINT1("Opening:%s\n", filename[i]);
 164                 mrcFileRead(&(in[i]), filename[i], "in lmrc2Dto2DMulti", 0);
 165                 if(Nx < in[i].HeaderN.x) {
 166                         Nx = in[i].HeaderN.x;
 167                 }
 168                 if(Ny < in[i].HeaderN.y) {
 169                         Ny = in[i].HeaderN.y;
 170                 }
 171                 if(in[i].HeaderLength.x < Length) {
 172                         Length = in[i].HeaderLength.x;
 173                 }
 174                 if(in[i].HeaderLength.y < Length) {
 175                         Length = in[i].HeaderLength.y;
 176                 }
 177         }
 178         out->HeaderN.x = Nx;
 179         out->HeaderN.y = Ny;
 180         out->HeaderN.z = number;
 181         out->HeaderMode = mrcFloatImage;
 182         out->HeaderLength.x = Length;
 183         out->HeaderLength.y = Length;
 184         out->HeaderLength.z = Length;
 185         mrcInit(out, NULL);
 186         srcz = 0;
 187         for(dstz=0; dstz<number; dstz++) {
 188                 for(dstx=0; dstx<in->HeaderN.x; dstx++) {
 189                         for(dsty=0; dsty<in->HeaderN.y; dsty++) {
 190                                 srcx = dstx - (out->HeaderN.x - in->HeaderN.x)/2.0;
 191                                 srcy = dsty - (out->HeaderN.y - in->HeaderN.y)/2.0;
 192                                 mrcPixelDataGet(&(in[(int)dstz]), srcx, srcy, srcz, &data, mrcPixelRePart, mrcPixelHowLinear);
 193                                 mrcPixelDataSet(out,              dstx, dsty, dstz,  data, mrcPixelRePart);
 194                         }
 195                 }
 196         }
 197         out->numTailer = number;
 198         out->Tailer = (mrcImageTailer*)memoryAllocate(sizeof(mrcImageTailer)*number, "in lmrc2Dto2DMulti");
 199         fseek(fptInfo, 0L, SEEK_SET);
 200         for(i=0; i<number; i++) {
 201                 stringGetFromFile(s, "", fptInfo, stdout, 3);
 202                 out->Tailer[i].Cont.Mode = mrcImageTailerMode2DProjection;
 203                 stringCopy(out->Tailer[i].Cont.EulerAngleMode, stringGetNthWord(s, 2, " \t,"), 4);
 204                 out->Tailer[i].Cont.Rot1 = stringGetNthRealData(s, 3, " ,\t")*RADIAN;
 205                 out->Tailer[i].Cont.Rot2 = stringGetNthRealData(s, 4, " ,\t")*RADIAN;
 206                 out->Tailer[i].Cont.Rot3 = stringGetNthRealData(s, 5, " ,\t")*RADIAN;
 207         }
 208         for(i=0; i<number; i++) {
 209                 mrcImageFree(&(in[i]), "in lmrc2Dto2DMulti");
 210         }
 211         free(in);
 212 }
 213
 214 void
 215 lmrc2Dto2DEach(lmrc2Dto3DInfo* linfo, char** filename, int number, FILE* fptInfo, int mode)
 216 {
 217         char s[1024];
 218         long i;
 219
 220         linfo->inFileNum = number;
 221         linfo->inFileList= filename;
 222         linfo->Tailer = (mrcImageTailer*)memoryAllocate(sizeof(mrcImageTailer)*number, "in lmrc2Dto2DEach");
 223         fseek(fptInfo, 0L, SEEK_SET);
 224         for(i=0; i<number; i++) {
 225                 stringGetFromFile(s, "", fptInfo, stdout, 3);
 226                 linfo->Tailer[i].Cont.Mode = mrcImageTailerMode2DProjection;
 227                 stringCopy(linfo->Tailer[i].Cont.EulerAngleMode, stringGetNthWord(s, 2, " \t,"), 4);
 228                 linfo->Tailer[i].Cont.Rot1 = stringGetNthRealData(s, 3, " ,\t")*RADIAN;
 229                 linfo->Tailer[i].Cont.Rot2 = stringGetNthRealData(s, 4, " ,\t")*RADIAN;
 230                 linfo->Tailer[i].Cont.Rot3 = stringGetNthRealData(s, 5, " ,\t")*RADIAN;
 231         }
 232 }
 233