00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "log.h"
00033 #include "emdata.h"
00034 #include "xydata.h"
00035 #include "assert.h"
00036 #include "projector.h"
00037
00038 #include <sys/types.h>
00039 #include <sys/times.h>
00040 #include <time.h>
00041 #include <sys/time.h>
00042
00043 #ifndef CLK_TCK
00044 #define CLK_TCK 60
00045 #endif
00046
00047 double mytimer()
00048 {
00049 struct tms use;
00050 double tmp;
00051 times(&use);
00052 tmp = use.tms_utime;
00053 tmp += use.tms_stime;
00054 return (double)(tmp) / CLK_TCK;
00055 }
00056
00057 using namespace EMAN;
00058
00059 using std::cout;
00060 using std::cin;
00061 using std::endl;
00062 using std::string;
00063
00064 #include "spidutil.h"
00065
00066 int main()
00067 {
00068 EMData *volume = new EMData();
00069 EMData *expimg = new EMData();
00070
00071 Dict refparams;
00072
00073 APMQopt options;
00074
00075 vector <float> angles;
00076 float delta, tlb, tub, plb, pub;
00077 int ndim, nx, ny, nz, nang, status;
00078 int nref, nimg, i;
00079 double t0, t1;
00080 float *newangles, *shifts;
00081
00082 printf("reading a 3-D volume...\n");
00083 volume->read_image("vol001.tfc");
00084
00085 ndim = volume->get_ndim();
00086 nx = volume->get_xsize();
00087 ny = volume->get_ysize();
00088 nz = volume->get_zsize();
00089
00090
00091 printf("ndim = %d, nx = %d, ny = %d, nz = %d\n", ndim, nx, ny, nz);
00092
00093
00094
00095 delta = 15.0;
00096 tlb = 0.0;
00097 tub = 90.0;
00098 plb = 0.0;
00099 pub = 359.9;
00100
00101
00102 angles = Util::voea(delta, tlb, tub, plb, pub);
00103 nang = angles.size()/3;
00104 printf("size(angles) = %d\n", nang);
00105
00106 refparams["anglelist"] = angles;
00107 refparams["angletype"] = "SPIDER";
00108 refparams["radius"] = 30.0;
00109 t0 = mytimer();
00110 printf("generating reference projections...\n");
00111 EMData* refprj = volume->project("pawel", refparams);
00112 t1 = mytimer() - t0;
00113 nref = refprj->get_zsize();
00114 printf("nref = %d\n", nref);
00115
00116 printf("time used in ref projection calculation = %11.3e\n", t1);
00117
00118
00119 expimg->read_image("tf2d0001.tfc");
00120 nimg = expimg->get_zsize();
00121
00122
00123 printf("running APMQ...\n");
00124
00125 options.mr = 5;
00126 options.nr = 25;
00127 options.shrange = 4;
00128 options.istep = 1;
00129 options.mode = 'F';
00130 options.range = 0.0;
00131 options.angexps = NULL;
00132
00133 newangles = (float*) calloc(3*nimg,sizeof(float));
00134 shifts = (float*) calloc(2*nimg,sizeof(float));
00135
00136 status = apmq(refprj, refparams, expimg, options, newangles, shifts);
00137 if (status != 0) {
00138 fprintf(stderr, "APMD failed! status = %d\n", status);
00139 goto EXIT;
00140 }
00141
00142 for (i = 1; i<=nimg; i++)
00143 printf("psi = %g, theta = %g, phi = %g, sx = %g, sy = %g\n",
00144 newangles(1,i), newangles(2,i), newangles(3,i),
00145 shifts(1,i), shifts(2,i));
00146
00147 EXIT:
00148 if (volume) delete volume;
00149 if (expimg) delete expimg;
00150 if (refprj) delete refprj;
00151 if (newangles) delete newangles;
00152 if (shifts) delete shifts;
00153 return status;
00154 }