00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <mpi.h>
00004 #include <math.h>
00005 #include <string.h>
00006
00007
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
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
00030
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
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
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
00054 dims[ROW] = atoi(argv[1]);
00055 dims[COL] = atoi(argv[2]);
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;
00063
00064 MPI_Cart_create(comm, 2, dims, periods, 1, &comm_2d);
00065 MPI_Comm_rank(comm_2d, &my2dpid);
00066 MPI_Cart_coords(comm_2d, my2dpid, 2, mycoords);
00067
00068
00069 keep_dims[ROW] = 0;
00070 keep_dims[COL] = 1;
00071 MPI_Cart_sub(comm_2d, keep_dims, &comm_row);
00072
00073
00074 keep_dims[ROW] = 1;
00075 keep_dims[COL] = 0;
00076 MPI_Cart_sub(comm_2d, keep_dims, &comm_col);
00077
00078
00079
00080 if (mycoords[ROW] == 0 && mycoords[COL] == 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
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
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
00120
00121 if (mycoords[COL] == 0 && mycoords[ROW] == 0) {
00122 for(int ip = 0; ip < dims[ROW]; ++ip){
00123 int begidx = nbase[ip]*nx*ny;
00124 if (ip !=0){
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{
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
00139 }
00140 }
00141 }
00142
00143 if (mycoords[COL] == 0 && mycoords[ROW] != 0) {
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
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
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
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
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
00210 MPI_Allreduce (bvol_loc, vol_sphloc, nnzloc, MPI_FLOAT, MPI_SUM, comm_col);
00211
00212
00213 if (mycoords[COL] != 0 && mycoords[ROW] == 0) {
00214
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
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
00234
00235
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
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
00278 MPI_Allreduce (newimagesloc, reprojloc, nangsloc*nx*ny, MPI_FLOAT, MPI_SUM, comm_row);
00279
00280 delete [] newimagesloc;
00281
00282 if (mycoords[ROW] != 0 && mycoords[COL] == 0) {
00283
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
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
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
00322 int idx;
00323 float err=0, max =0;
00324 for (int i=0; i< nangs*nx*ny; i++){
00325
00326
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 }