EMAN2
runcartrec.cpp
Go to the documentation of this file.
00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <mpi.h>
00004 #include <math.h>
00005 #include <string.h>
00006 
00007 // GLOBAL Variables
00008 int nangs, nangsloc, nx, ny, myptrstart;
00009 int nnz, nrays, nnzloc, nraysloc;
00010 float *imagesloc, *images;
00011 float *newimagesloc, *reprojloc;
00012 float *angles, *anglesloc;
00013 float *vol_sph;
00014 float *onevol_sph, *reproj;
00015 
00016 // Additional functions that are found in backproj_Cart.cpp
00017 int setpart_gc1(MPI_Comm comm_2d, int nangs, int *psize, int *nbase);
00018 int getnnz(int *volsize, int ri, int *origin, int *nrays, int *nnz);
00019 int setpart_gr1(MPI_Comm comm_2d, int nnz, int *nnzpart, int *nnzbase);
00020 int getcb2sph(int *volsize, int ri, int *origin, int nnz0, int *ptrs, int *cord);
00021 int sphpart(MPI_Comm comm_2d, int nrays, int *ptrs, int *nnzbase, int *ptrstart);
00022 int make_proj_mat(float phi, float theta, float psi, float * dm);
00023 int ifix(float a);
00024 int fwdpj3_Cart(int *volsize, int nraysloc, int nnzloc, float *dm, int *origin, int ri, int *ptrs, int *cord, int myptrstart, float *x, float *y);
00025 int bckpj3_Cart(int *volsize, int nraysloc, int nnzloc, float *dm, int *origin, int ri, int *ptrs, int *cord, int myptrstart, float *x, float *y);
00026 int fwdpj3(int *volsize, int nrays, int   nnz, float *dm, int *origin, int ri, int *ptrs, int *cord, float *x, float  *y);
00027 int bckpj3(int *volsize, int nrays, int   nnz, float *dm, int *origin, int ri, int *ptrs, int *cord, float *x, float *y);
00028 
00029 // ---------------- MAIN FUNCTION --------------------------------
00030 // This code will read in the 2D images and perform the backprojection operation in parallel, and then a parallel forward projection operation
00031 
00032 int main(int argc, char *argv[])
00033 {
00034   int ncpus, mypid, nrem, ierr;
00035   MPI_Status mpistatus;
00036   MPI_Comm comm = MPI_COMM_WORLD;
00037   
00038   // Variables needed for the Cartesian topology.
00039   int ROW = 0, COL = 1;
00040   int dims[2], periods[2], keep_dims[2];
00041   int my2dpid, mycoords[2], srcoords[2], otherpid;
00042   MPI_Comm comm_2d, comm_row, comm_col; 
00043                   
00044 // Initialize MPI.
00045   MPI_Init(&argc, &argv);
00046   MPI_Comm_size(comm, &ncpus);
00047   MPI_Comm_rank(comm, &mypid);
00048   
00049   if ( argc < 3  ) {
00050           printf ("ERROR: %s requires Cartesian dimensions input\n", argv[0]);
00051           return -1;
00052   }
00053 // Set up a Cartesian virtual topology and get the rank and coordinates of the processes in the topology. 
00054   dims[ROW] = atoi(argv[1]); // Row dimension of the topology
00055   dims[COL] = atoi(argv[2]); // Column dimension of the topology
00056   
00057   if (dims[ROW]*dims[COL] != ncpus){
00058         printf("ERROR: Row dim and col dim not equal to ncpus\n");
00059         return -1;
00060   }
00061   
00062   periods[ROW] = periods[COL] = 1; // Set the periods for wrap-around
00063   
00064   MPI_Cart_create(comm, 2, dims, periods, 1, &comm_2d);
00065   MPI_Comm_rank(comm_2d, &my2dpid); //Get my pid in the new 2D topology
00066   MPI_Cart_coords(comm_2d, my2dpid, 2, mycoords); // Get my coordinates
00067   
00068   /* Create the row-based sub-topology */ 
00069   keep_dims[ROW] = 0; 
00070   keep_dims[COL] = 1; 
00071   MPI_Cart_sub(comm_2d, keep_dims, &comm_row); 
00072  
00073   /* Create the column-based sub-topology */ 
00074   keep_dims[ROW] = 1; 
00075   keep_dims[COL] = 0; 
00076   MPI_Cart_sub(comm_2d, keep_dims, &comm_col); 
00077 
00078 // STEP 1: Have processor (0,0) read in the entire set of 2D images, divide up the images, and send corresponding images to processors in processor group: g_c_0  Do the same for the angles.
00079   
00080   if (mycoords[ROW] == 0 && mycoords[COL] == 0){ //I'm processor (0,0)
00081     FILE *fp, *fpa;
00082     char imagefname[80]="tf2d84.raw", anglesfname[80]="angles.dat";
00083           
00084     fp = fopen(imagefname,"r");
00085     fread(&nangs, sizeof(int), 1, fp);
00086     fread(&nx, sizeof(int), 1, fp);
00087     fread(&ny, sizeof(int), 1, fp);
00088     
00089     images = new float[nx*ny*nangs];
00090     fread(images, sizeof(float), nx*ny*nangs, fp);
00091     fclose(fp);
00092     
00093     fpa = fopen(anglesfname,"r");
00094     angles = new float[3*nangs];
00095     for (int i = 0; i< 3*nangs; i++)
00096       fscanf(fpa, "%f",&angles[i]);
00097        
00098     fclose(fpa);
00099     printf("There are %d 2D images of size %d x %d\n", nangs, nx, ny);
00100   }
00101   
00102   // Broadcast variables nangs, nx, ny to all processors
00103   srcoords[ROW] = srcoords[COL] = 0;
00104   MPI_Cart_rank(comm_2d, srcoords, &otherpid); 
00105   
00106   MPI_Bcast (&nangs, 1, MPI_INT, otherpid, comm_2d);
00107   MPI_Bcast (&nx, 1, MPI_INT, otherpid, comm_2d);
00108   MPI_Bcast (&ny, 1, MPI_INT, otherpid, comm_2d);
00109   
00110   // Send images and angles from Processor (0,0) to processors in group g_c_0
00111   int *psize = new int[dims[ROW]];
00112   int *nbase = new int[dims[ROW]];
00113   
00114   nangsloc = setpart_gc1(comm_2d, nangs, psize, nbase);
00115   imagesloc = new float[psize[mycoords[ROW]]*nx*ny];
00116   reprojloc = new float[psize[mycoords[ROW]]*nx*ny];
00117   anglesloc = new float[psize[mycoords[ROW]]*3];
00118   
00119 // printf("My coords are (%d,%d) and nangsloc = %d\n", mycoords[ROW], mycoords[COL], nangsloc);
00120   
00121   if (mycoords[COL] == 0 && mycoords[ROW] == 0) { //I'm Proc. (0,0)
00122     for(int ip = 0; ip < dims[ROW]; ++ip){
00123       int begidx = nbase[ip]*nx*ny;
00124       if (ip !=0){ // Proc (0,0) sends images and angle data to other processors
00125          srcoords[COL] = 0;
00126          srcoords[ROW] = ip;
00127          MPI_Cart_rank(comm_2d, srcoords, &otherpid);
00128          MPI_Send(&images[begidx],psize[ip]*nx*ny, MPI_FLOAT, otherpid, otherpid, comm_2d);
00129          MPI_Send(&angles[nbase[ip]*3],psize[ip]*3, MPI_FLOAT, otherpid, otherpid, comm_2d);
00130       }
00131       else{ // ip = 0: Proc (0,0) needs to copy images and angles into its imagesloc and anglesloc
00132         for (int i = 0; i < psize[ip]*nx*ny; i++){
00133           imagesloc[i] = images[begidx+i];
00134         }
00135         for (int i = 0; i < psize[ip]*3; i++){
00136                 anglesloc[i] = angles[nbase[ip]*3 + i];
00137         }
00138         //printf("Finished copying to Proc (0,0) local");
00139       }
00140     } //End for loop
00141   } //End if
00142   
00143   if (mycoords[COL] == 0 && mycoords[ROW] != 0) { //I'm in g_c_0 and I'm not Processor (0,0) so I should receive data.
00144     MPI_Recv(imagesloc, psize[mycoords[ROW]]*nx*ny, MPI_FLOAT, 0, mypid, comm_2d, &mpistatus);
00145     MPI_Recv(anglesloc, psize[mycoords[ROW]]*3, MPI_FLOAT, 0, mypid, comm_2d, &mpistatus);
00146   }
00147   // Now have all the processors in group g_c_0 broadcast the images and angles along the row communicator
00148   srcoords[ROW] = 0;
00149   MPI_Cart_rank(comm_row, srcoords, &otherpid);
00150   MPI_Bcast(imagesloc, nangsloc*nx*ny, MPI_FLOAT, otherpid , comm_row);
00151   MPI_Bcast(anglesloc, nangsloc*3, MPI_FLOAT, otherpid , comm_row);
00152  
00153 // Now distribute the volume (in spherical format) among columns of processors and use nnz to determine the splitting.  Note: ptrs and coord are on all processors
00154   int radius;
00155   int volsize[3], origin[3];
00156   volsize[0] = nx;
00157   volsize[1] = nx;
00158   volsize[2] = nx;
00159   origin[0] = nx/2+1;
00160   origin[1] = nx/2+1;
00161   origin[2] = nx/2+1;
00162   radius = nx/2-1;
00163    
00164   ierr = getnnz( volsize, radius, origin, &nrays, &nnz);
00165    
00166   int * ptrs = new int[nrays+1];
00167   int * cord = new int[3*nrays];
00168   ierr = getcb2sph(volsize, radius, origin, nnz, ptrs, cord);
00169                   
00170   int *nnzpart = new int[dims[COL]];
00171   int *nnzbase = new int[dims[COL]+1]; 
00172   nnzloc = setpart_gr1(comm_2d, nnz, nnzpart, nnzbase);
00173   
00174   int *ptrstart = new int[dims[COL]+1];
00175   nraysloc = sphpart(comm_2d, nrays, ptrs, nnzbase, ptrstart);
00176   
00177   myptrstart = ptrstart[mycoords[COL]];
00178   int nnzall[dims[COL]];
00179   for (int i = 0; i<dims[COL]; i++)
00180     nnzall[i] = ptrs[ptrstart[i+1]] - ptrs[ptrstart[i]];
00181   
00182   nnzloc = nnzall[mycoords[COL]];
00183   
00184   // Print some stuff.
00185  printf("My coords are (%d,%d) and nangsloc = %d, nraysloc = %d, myptrstart = %d, nnzloc = %d\n", mycoords[ROW], mycoords[COL], nangsloc, nraysloc, myptrstart, nnzloc);
00186   
00187   float *bvol_loc = new float[nnzloc];
00188   float *vol_sphloc = new float[nnzloc];
00189   for (int i=0; i< nnzloc; i++)
00190     bvol_loc[i] = 0.0;
00191   
00192   // STEP 2: Have everyone perform the backprojection operation for their assigned images and portions of the volume.  Then perform an Allreduce along the columns.
00193   
00194   float phi, theta, psi;
00195   float dm[8];
00196   
00197   for (int i=0; i<nangsloc; i++){
00198     phi = anglesloc[3*i+0];
00199     theta = anglesloc[3*i+1];
00200     psi = anglesloc[3*i+2];
00201     dm[6] = 0;
00202     dm[7] = 0;
00203  
00204     make_proj_mat(phi, theta, psi, dm);
00205 
00206     ierr = bckpj3_Cart(volsize, nraysloc, nnzloc, dm, origin, radius, ptrs, cord, myptrstart, &imagesloc[nx*ny*i], bvol_loc);
00207   }
00208   
00209   // Now an all reduce along the columns
00210   MPI_Allreduce (bvol_loc, vol_sphloc, nnzloc, MPI_FLOAT, MPI_SUM, comm_col);
00211   
00212   // For testing purposes, we bring all the portions of the volume back together onto Proc (0,0). Note: we only need to deal with the first row of processors.
00213   if (mycoords[COL] != 0 && mycoords[ROW] == 0) {
00214         //Send data to Processor (0,0)
00215         srcoords[COL] = srcoords[ROW] = 0;
00216         MPI_Cart_rank(comm_2d, srcoords, &otherpid);
00217         MPI_Send(vol_sphloc, nnzloc, MPI_FLOAT, otherpid, otherpid, comm_2d);
00218   }
00219   float *onevol_sph = new float[nnz];
00220   if (mycoords[COL] == 0 && mycoords[ROW] ==0){
00221           //Copy data and recieve data
00222         float *vol_sph = new float[nnz];
00223         for (int i=0; i<nnzloc; i++)
00224           vol_sph[i] = vol_sphloc[i];
00225         
00226         for (int i=1; i<dims[COL]; i++){
00227           srcoords[ROW] = 0;
00228           srcoords[COL] = i;
00229           MPI_Cart_rank(comm_2d, srcoords, &otherpid);
00230                 
00231           MPI_Recv(&vol_sph[ptrs[ptrstart[i]]-1], nnzall[i], MPI_FLOAT, otherpid, mypid, comm_2d, &mpistatus);
00232         }
00233         //printf("Finished combining all volume parts\n");
00234   
00235   //Now compute the back projection serially on one processor (0,0)
00236     for (int i=0; i< nnz; i++)
00237       onevol_sph[i] = 0.0;
00238     
00239     for (int i=0; i<nangs; i++){
00240             phi = angles[3*i+0];
00241             theta = angles[3*i+1];
00242             psi = angles[3*i+2];
00243             dm[6] = 0;
00244             dm[7] = 0;
00245  
00246             make_proj_mat(phi, theta, psi, dm);
00247   
00248             ierr = bckpj3(volsize, nrays, nnz, dm, origin, radius, ptrs, cord, &images[nx*ny*i], onevol_sph);
00249     }
00250     
00251     float err=0;
00252     for (int i=0; i< nnz; i++){
00253             err = err+(onevol_sph[i]-vol_sph[i])*(onevol_sph[i]-vol_sph[i]);
00254     }
00255     err = sqrt(err);
00256     printf("Cumulative error for backprojection is %f\n", err);
00257     delete [] vol_sph;
00258   }
00259   
00260     // STEP 3: Now perform a forward projection operation for the assigned images and portions of the volume.  Then perform an all_reduce along the rows.
00261   float * newimagesloc = new float[nangsloc*nx*ny];
00262   for (int i=0; i<nangsloc*nx*ny; i++)
00263           newimagesloc[i] = 0.0;
00264         
00265   for (int i=0; i<nangsloc; i++){
00266   phi = anglesloc[3*i+0];
00267           theta = anglesloc[3*i+1];
00268           psi = anglesloc[3*i+2];
00269           dm[6] = 0;
00270           dm[7] = 0;
00271     
00272           make_proj_mat(phi, theta, psi, dm);
00273   
00274          ierr = fwdpj3_Cart(volsize, nraysloc, nnzloc, dm, origin, radius, ptrs, cord, myptrstart, vol_sphloc, &newimagesloc[nx*ny*i]);
00275   }
00276 
00277   // Now an all reduce along the rows
00278   MPI_Allreduce (newimagesloc, reprojloc, nangsloc*nx*ny, MPI_FLOAT, MPI_SUM, comm_row);
00279  
00280   delete [] newimagesloc;
00281   // For testing purposes, we bring all the 2D images together onto Proc (0,0). Note: we only need to deal with the first column of processors.
00282   if (mycoords[ROW] != 0 && mycoords[COL] == 0) {
00283         //Send data to Processor (0,0)
00284           srcoords[COL] = srcoords[ROW] = 0;
00285           MPI_Cart_rank(comm_2d, srcoords, &otherpid);
00286           MPI_Send(reprojloc, nangsloc*nx*ny, MPI_FLOAT, otherpid, otherpid, comm_2d);
00287   }
00288   if (mycoords[COL] == 0 && mycoords[ROW] ==0){
00289           //Copy data and recieve data
00290           float *reproj = new float[nangs*nx*ny];
00291           for (int i=0; i<nangsloc*nx*ny; i++)
00292                   reproj[i] = reprojloc[i];
00293         
00294           for (int i=1; i<dims[ROW]; i++){
00295                   srcoords[COL] = 0;
00296                   srcoords[ROW] = i;
00297                   MPI_Cart_rank(comm_2d, srcoords, &otherpid);
00298                 
00299                   MPI_Recv(&reproj[nbase[i]*nx*ny], psize[i]*nx*ny, MPI_FLOAT, otherpid, mypid, comm_2d, &mpistatus);
00300           }
00301   delete [] reprojloc;
00302   // Now compute the forward projection serially on one processor (0,0)
00303 
00304           float *allimages = new float[nangs*nx*ny];
00305           for (int i=0; i< nangs*nx*ny; i++)
00306                   allimages[i] = 0.0;
00307     
00308           for (int i=0; i<nangs; i++){
00309             phi = angles[3*i+0];
00310             theta = angles[3*i+1];
00311             psi = angles[3*i+2];
00312             dm[6] = 0;
00313             dm[7] = 0;
00314     
00315             make_proj_mat(phi, theta, psi, dm);
00316                   
00317             ierr = fwdpj3(volsize, nrays, nnz, dm, origin, radius, ptrs, cord, onevol_sph, &allimages[nx*ny*i]);
00318 
00319           }
00320 
00321           //  Now compute the overall error.
00322 int idx;        
00323   float err=0, max =0;
00324           for (int i=0; i< nangs*nx*ny; i++){
00325                 //if (allimages[i]!=reproj[i] && i < 256)
00326                 //      printf("i= %d\n",i);
00327                 
00328             err = err+(allimages[i]-reproj[i])*(allimages[i]-reproj[i]);
00329 
00330                 if (fabs(allimages[i]-reproj[i]) > max){
00331                         max = fabs(allimages[i]-reproj[i]);
00332 idx = i;
00333 }
00334           }
00335 err = sqrt(err);  
00336           printf("Cumulative error for forward projection is %f with max error of %f occuring at %d\n", err, max, idx);
00337 printf("Max error: compare %f and %f\n",allimages[idx], reproj[idx]);
00338 
00339           delete [] reproj;
00340           delete [] allimages;
00341           delete [] angles;
00342           delete [] images;
00343   }
00344   delete [] onevol_sph;
00345   delete [] vol_sphloc;
00346   delete [] bvol_loc;  
00347   delete [] ptrs;
00348   delete [] cord;
00349   delete [] nnzpart;
00350   delete [] nnzbase;
00351   delete [] ptrstart;
00352   delete [] anglesloc;
00353   delete [] imagesloc;
00354   delete [] nbase;
00355   delete [] psize;        
00356   
00357   MPI_Comm_free(&comm_2d);
00358   MPI_Comm_free(&comm_row);
00359   MPI_Comm_free(&comm_col);
00360   
00361   MPI_Finalize();
00362 }