#include "mpi.h"
#include "stdlib.h"
#include "emdata.h"
#include "alignoptions.h"
Include dependency graph for utilcomm.h:
This graph shows which files directly or indirectly include this file:
Go to the source code of this file.
Functions | |
int | ReadVandBcast (MPI_Comm comm, EMData *volume, char *volfname) |
int | ReadStackandDist (MPI_Comm comm, EMData ***images2D, char *stackfname, int *nloc) |
int | ReadAngTrandDist (MPI_Comm comm, float *angleshift, char *paramfname, int nloc) |
int | CleanStack (MPI_Comm comm, EMData **image_stack, int nloc, int ri, Vec3i volsize, Vec3i origin) |
int | setpart (MPI_Comm comm, int nima, int *psize, int *nbase) |
int | ParseAlignOptions (MPI_Comm comm, AlignOptions &options, char *optionsfname, int nvoxels, EMData *&mask3D) |
int | asta2 (float *img, int nx, int ny, int ri, double *abaloc, int *klploc) |
|
Definition at line 12 of file utilnum.cpp. 00013 { 00014 // calculate the sum of the background pixels, and then set the intensity 00015 // of these pixels to zero. Also count the number of background pixels. 00016 // A background pixel is a pixel outside of the circle with radius ri 00017 // This is done for images assigned to a processor 00018 // 00019 int xcent = (nx / 2) + 1; 00020 int ycent = (ny / 2) + 1; 00021 int r_squared = ri*ri; 00022 00023 int x_summand, y_summand; 00024 00025 for ( int i = 0 ; i < nx ; ++i ) { 00026 x_summand = (i-xcent) * (i-xcent); 00027 for ( int j = 0 ; j < ny ; ++j ) { 00028 y_summand = (j-ycent) * (j-ycent); 00029 if ( x_summand + y_summand > r_squared ) { 00030 *abaloc += (double) img[j*nx + i]; 00031 //chao set the background to zero 00032 img[j*nx+i]=0.0; 00033 ++*klploc; 00034 } 00035 } 00036 } 00037 00038 return 0; 00039 }
|
|
Definition at line 233 of file utilcomm.cpp. References asta2(), EMDeleteArray(), EMAN::EMData::get_data(), imgdata, nx, ny, rhs, and EMAN::Vec3i. Referenced by main(). 00234 { 00235 int nx = volsize[0]; 00236 int ny = volsize[1]; 00237 00238 float * rhs = new float[nx*ny]; 00239 float * imgdata; 00240 00241 // Calculate average "background" from all pixels strictly outside of radius 00242 double aba, abaloc; // average background 00243 aba = 0.0; 00244 abaloc = 0.0; 00245 int klp, klploc; // number of pixels in background 00246 klp = 0; 00247 klploc = 0; 00248 00249 // calculate avg background in parallel 00250 for ( int i = 0 ; i < nloc ; ++i ) { 00251 imgdata = image_stack[i]->get_data(); 00252 asta2(imgdata, nx, ny, ri, &abaloc, &klploc); 00253 } 00254 00255 MPI_Allreduce(&abaloc, &aba, 1, MPI_DOUBLE, MPI_SUM, comm); 00256 MPI_Allreduce(&klploc, &klp, 1, MPI_INT, MPI_SUM, comm); 00257 00258 aba /= klp; 00259 00260 // subtract off the average background from pixels weakly inside of radius 00261 int x_summand, y_summand; 00262 int r_squared = ri * ri; 00263 for ( int i = 0 ; i < nloc ; ++i ) { 00264 imgdata = image_stack[i]->get_data(); 00265 for ( int j = 0 ; j < nx ; ++j) { 00266 x_summand = (j - origin[0]) * (j - origin[0]); 00267 for ( int k = 0 ; k < ny ; ++k ) { 00268 y_summand = (k - origin[1]) * (k - origin[1]); 00269 if ( x_summand + y_summand <= r_squared) { 00270 imgdata[j*ny + k] -= aba; 00271 } 00272 } 00273 } 00274 } 00275 EMDeleteArray(rhs); 00276 return 0; 00277 }
|
|
Definition at line 8 of file utilparse.cpp. 00010 { 00011 int ncpus, mypid, ierr; 00012 00013 MPI_Comm_size(comm,&ncpus); 00014 MPI_Comm_rank(comm,&mypid); 00015 00016 std::ifstream option_stream; 00017 std::string current_option; 00018 char option_buffer[100]; 00019 int int_option, parserr = 0; 00020 float float_option; 00021 00022 float * mask_data; 00023 if ( mypid == 0 ) { // read data from the options file 00024 option_stream.open(optionsfname); 00025 while ( option_stream >> current_option ) { 00026 if ( current_option == "maskfile" ) { 00027 option_stream >> current_option; 00028 mask3D = new EMData(); 00029 mask3D->read_image(current_option); 00030 options.set_mask3D(mask3D); 00031 } 00032 else if ( current_option == "inner_ring" ) { 00033 option_stream >> int_option; 00034 options.set_first_ring(int_option); 00035 std::cout << "first_ring = " << int_option << std::endl; 00036 } 00037 else if ( current_option == "outer_ring" ) { 00038 option_stream >> int_option; 00039 options.set_last_ring(int_option); 00040 std::cout << "last_ring = " << int_option << std::endl; 00041 } 00042 else if ( current_option == "rstep" ) { 00043 option_stream >> int_option; 00044 options.set_rstep(int_option); 00045 std::cout << "rstep = " << int_option << std::endl; 00046 } 00047 else if ( current_option == "radius" ) { // perhaps this is the same as outer_ring? 00048 option_stream >> int_option; 00049 options.set_ri(int_option); 00050 std::cout << "radius = " << int_option << std::endl; 00051 } 00052 else if ( current_option == "x_range" ) { 00053 option_stream >> float_option; 00054 options.set_xrng(float_option); 00055 std::cout << "x_range = " << float_option << std::endl; 00056 } 00057 else if ( current_option == "y_range" ) { 00058 option_stream >> float_option; 00059 options.set_yrng(float_option); 00060 std::cout << "y_range = " << float_option << std::endl; 00061 } 00062 else if ( current_option == "translation_step" ) { 00063 option_stream >> float_option; 00064 options.set_step(float_option); 00065 std::cout << "step = " << float_option << std::endl; 00066 } 00067 else if ( current_option == "theta_step" ) { 00068 option_stream >> float_option; 00069 options.set_dtheta(float_option); 00070 std::cout << "theta_step = " << float_option << std::endl; 00071 } 00072 else if ( current_option == "CTF" ) { 00073 option_stream >> current_option; 00074 if ( current_option == "true" ) { 00075 options.set_CTF(true); 00076 std::cout << "CTF = true" << std::endl; 00077 } 00078 else { // anything else sets it to false 00079 options.set_CTF(false); 00080 std::cout << "CTF = false" << std::endl; 00081 } 00082 } 00083 else if ( current_option == "snr" ) { 00084 option_stream >> float_option; 00085 options.set_snr(float_option); 00086 std::cout << "snr = " << float_option << std::endl; 00087 } 00088 else if ( current_option == "ref_a" ) { 00089 option_stream >> current_option; 00090 if ( current_option == "P" ) { 00091 options.set_ref_angle_type("P"); 00092 } 00093 else if ( current_option == "S" ){ // should add support for this 00094 std::cerr << "Currently only support Penczek-type reference angles..." << std::endl; 00095 options.set_ref_angle_type("P"); 00096 } 00097 else { // Currently default to "P", will eventually default to "S" 00098 } 00099 } 00100 else if ( current_option == "symmetry" ) { 00101 option_stream >> current_option; 00102 options.set_symmetry(current_option); 00103 std::cout << "symmetry = " << current_option << std::endl; 00104 } 00105 else if ( current_option == "use_sirt" ) { 00106 option_stream >> current_option; 00107 if ( current_option == "true" ) { 00108 options.set_use_sirt(true); 00109 std::cout << "use_sirt = true" << std::endl; 00110 } 00111 else { // anything else sets it to false 00112 options.set_use_sirt(false); 00113 std::cout << "use_sirt = false" << std::endl; 00114 } 00115 } 00116 else if ( current_option == "sirt_tol" ) { 00117 option_stream >> float_option; 00118 options.set_sirt_tol(float_option); 00119 std::cout << "sirt_tol = " << float_option << std::endl; 00120 } 00121 else if ( current_option == "sirt_lam" ) { 00122 option_stream >> float_option; 00123 options.set_sirt_lam(float_option); 00124 std::cout << "sirt_lam = " << float_option << std::endl; 00125 } 00126 else if ( current_option == "sirt_maxit" ) { 00127 option_stream >> int_option; 00128 options.set_sirt_maxit(int_option); 00129 std::cout << "sirt_maxit = " << int_option << std::endl; 00130 } 00131 else if ( current_option == "maxit" ) { 00132 option_stream >> int_option; 00133 options.set_maxit(int_option); 00134 std::cout << "maxit = " << int_option << std::endl; 00135 } 00136 else { // catch-all 00137 std::cerr << "Unsupported option " << current_option << "..." << std::endl; 00138 parserr = 1; 00139 } 00140 } 00141 option_stream.close(); 00142 if (parserr) { 00143 printf("The supported options are:\n"); 00144 printf(" maskfile \n"); 00145 printf(" inner_ring\n"); 00146 printf(" outer_ring\n"); 00147 printf(" rstep\n"); 00148 printf(" radius\n"); 00149 printf(" x_range\n"); 00150 printf(" y_range\n"); 00151 printf(" translation_step\n"); 00152 printf(" CTF\n"); 00153 printf(" snr\n"); 00154 printf(" ref_a\n"); 00155 printf(" S\n"); 00156 printf(" symmetry\n"); 00157 printf(" sirt_tol\n"); 00158 printf(" sirt_lam\n"); 00159 printf(" sirt_maxit\n"); 00160 printf(" maxit\n"); 00161 } 00162 } 00163 // Then broadcast all the data that was read 00164 ierr = MPI_Bcast(&mask3D, 1, MPI_INT, 0, comm); 00165 // if it's not NULL, need to allocate and bcast its data 00166 // NOTE: this is sending over the master's address for mask3D ONLY as a test to see if it's NULL or not. 00167 // DO NOT USE THIS POINTER ON ANY OTHER NODES EXCEPT FOR THIS TEST! 00168 if ( mask3D != NULL ) { 00169 if ( mypid != 0 ) { 00170 mask3D = new EMData(); 00171 } 00172 mask_data = mask3D->get_data(); 00173 ierr = MPI_Bcast(mask_data, nvoxels, MPI_FLOAT, 0, comm); 00174 options.set_mask3D(mask3D); 00175 } 00176 00177 int_option = options.get_first_ring(); 00178 ierr = MPI_Bcast(&int_option, 1, MPI_INT, 0, comm); 00179 options.set_first_ring(int_option); 00180 00181 int_option = options.get_last_ring(); 00182 ierr = MPI_Bcast(&int_option, 1, MPI_INT, 0, comm); 00183 options.set_last_ring(int_option); 00184 00185 int_option = options.get_rstep(); 00186 ierr = MPI_Bcast(&int_option, 1, MPI_INT, 0, comm); 00187 options.set_rstep(int_option); 00188 00189 int_option = options.get_ri(); 00190 ierr = MPI_Bcast(&int_option, 1, MPI_INT, 0, comm); 00191 options.set_ri(int_option); 00192 00193 float_option = options.get_xrng(); 00194 ierr = MPI_Bcast(&float_option, 1, MPI_FLOAT, 0, comm); 00195 options.set_xrng(float_option); 00196 00197 float_option = options.get_yrng(); 00198 ierr = MPI_Bcast(&float_option, 1, MPI_FLOAT, 0, comm); 00199 options.set_yrng(float_option); 00200 00201 float_option = options.get_step(); 00202 ierr = MPI_Bcast(&float_option, 1, MPI_FLOAT, 0, comm); 00203 options.set_step(float_option); 00204 00205 float_option = options.get_dtheta(); 00206 ierr = MPI_Bcast(&float_option, 1, MPI_FLOAT, 0, comm); 00207 options.set_dtheta(float_option); 00208 00209 int_option = options.get_CTF(); 00210 ierr = MPI_Bcast(&int_option, 1, MPI_INT, 0, comm); 00211 options.set_CTF(int_option); 00212 00213 float_option = options.get_snr(); 00214 ierr = MPI_Bcast(&float_option, 1, MPI_FLOAT, 0, comm); 00215 options.set_snr(float_option); 00216 00217 current_option = options.get_ref_angle_type(); 00218 strcpy(option_buffer, current_option.c_str()); 00219 ierr = MPI_Bcast(option_buffer, current_option.size(), MPI_CHAR, 0, comm); 00220 options.set_ref_angle_type(option_buffer); 00221 00222 current_option = options.get_symmetry(); 00223 strcpy(option_buffer, current_option.c_str()); 00224 ierr = MPI_Bcast(option_buffer, current_option.size(), MPI_CHAR, 0, comm); 00225 options.set_symmetry(option_buffer); 00226 00227 int_option = options.get_use_sirt(); 00228 ierr = MPI_Bcast(&int_option, 1, MPI_INT, 0, comm); 00229 options.set_use_sirt(int_option); 00230 00231 float_option = options.get_sirt_tol(); 00232 ierr = MPI_Bcast(&float_option, 1, MPI_FLOAT, 0, comm); 00233 options.set_sirt_tol(float_option); 00234 00235 float_option = options.get_sirt_lam(); 00236 ierr = MPI_Bcast(&float_option, 1, MPI_FLOAT, 0, comm); 00237 options.set_sirt_lam(float_option); 00238 00239 int_option = options.get_sirt_maxit(); 00240 ierr = MPI_Bcast(&int_option, 1, MPI_INT, 0, comm); 00241 options.set_sirt_maxit(int_option); 00242 00243 int_option = options.get_maxit(); 00244 ierr = MPI_Bcast(&int_option, 1, MPI_INT, 0, comm); 00245 options.set_maxit(int_option); 00246 00247 return 0; // ParseAlignOptions 00248 }
|
|
Definition at line 157 of file utilcomm.cpp. References EMDeleteArray(), and ierr. Referenced by main(). 00158 { 00159 // read a parameter file that contains a list of angles and shifts and distribute 00160 // them evenly among different processors 00161 00162 int mypid, ierr=0, mpierr=0, nimgs=0, ncpus; 00163 double t0; 00164 MPI_Status mpistatus; 00165 FILE *fp=NULL; 00166 00167 float *iobuffer = new float[5*nloc]; 00168 if (!iobuffer) { 00169 fprintf(stderr,"failed to allocate buffer to read angles shifts\n"); 00170 ierr = -1; 00171 goto EXIT; 00172 } 00173 00174 MPI_Comm_rank(comm,&mypid); 00175 MPI_Comm_size(comm,&ncpus); 00176 00177 t0 = MPI_Wtime(); 00178 00179 if (mypid ==0) { 00180 fp = fopen(paramfname,"r"); 00181 if (!fp) ierr = 1; 00182 } 00183 MPI_Bcast(&ierr, 1, MPI_INT, 0, comm); 00184 00185 if ( ierr ) { 00186 if (mypid ==0) fprintf(stderr,"failed to open %s\n", paramfname); 00187 } 00188 else { 00189 if (mypid == 0) { 00190 for (int iproc = 0; iproc < ncpus; iproc++) { 00191 // figure out the number of images assigned to processor iproc 00192 if (iproc > 0) { 00193 MPI_Recv(&nimgs, 1, MPI_INT, iproc, iproc, comm, &mpistatus); 00194 // Read the next nimgs set of angles and shifts 00195 for (int i = 0; i < nimgs; i++) { 00196 fscanf(fp,"%f %f %f %f %f", 00197 &iobuffer[5*i+0], 00198 &iobuffer[5*i+1], 00199 &iobuffer[5*i+2], 00200 &iobuffer[5*i+3], 00201 &iobuffer[5*i+4]); 00202 } 00203 MPI_Send(iobuffer,5*nimgs,MPI_FLOAT,iproc,iproc,comm); 00204 } 00205 else { 00206 for (int i = 0; i < nloc; i++) { 00207 fscanf(fp,"%f %f %f %f %f", 00208 &angleshift[5*i+0], 00209 &angleshift[5*i+1], 00210 &angleshift[5*i+2], 00211 &angleshift[5*i+3], 00212 &angleshift[5*i+4]); 00213 } 00214 } 00215 } 00216 fclose(fp); 00217 } 00218 else { 00219 // send image count to the master processor (mypid = 0) 00220 MPI_Send(&nloc, 1, MPI_INT, 0, mypid, comm); 00221 // Receive angleshifts 00222 MPI_Recv(angleshift, 5*nloc, MPI_FLOAT, 0, mypid, comm, &mpistatus); 00223 } 00224 } 00225 EMDeleteArray(iobuffer); 00226 if (mypid == 0) 00227 printf("I/O time for reading angles & shifts = %11.3e\n", 00228 MPI_Wtime() - t0); 00229 EXIT: 00230 return ierr; 00231 }
|
|
Definition at line 48 of file utilcomm.cpp. References EMDeleteArray(), EMDeletePtr(), EMAN::EMData::get_attr(), EMAN::EMData::get_data(), EMAN::EMUtil::get_image_count(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), ierr, img_ptr, imgdata, nx, ny, EMAN::EMData::read_image(), and setpart(). Referenced by main(). 00049 { 00050 int ncpus, mypid, ierr, mpierr=0; 00051 MPI_Status mpistatus; 00052 EMUtil *my_util; 00053 int nima; 00054 FILE *fp=NULL; 00055 00056 MPI_Comm_size(comm,&ncpus); 00057 MPI_Comm_rank(comm,&mypid); 00058 00059 ierr = 0; 00060 if (mypid == 0) { 00061 fp = fopen(stackfname,"r"); 00062 if (!fp) { 00063 ierr = 1; 00064 printf("failed to open %s\n", stackfname); 00065 } 00066 else { 00067 fclose(fp); 00068 nima = my_util->get_image_count(stackfname); 00069 } 00070 } 00071 mpierr = MPI_Bcast(&ierr, 1, MPI_INT, 0, comm); 00072 if (ierr == 0) { 00073 int *psize = new int[ncpus]; 00074 int *nbase = new int[ncpus]; 00075 00076 mpierr = MPI_Bcast(&nima, 1, MPI_INT, 0, comm); 00077 00078 *nloc = setpart(comm, nima, psize, nbase); 00079 *images2D = new EMData*[*nloc]; // NB!: whoever calls ReadStackandDist must delete this! 00080 00081 EMData *img_ptr; 00082 int img_index; 00083 float *imgdata; 00084 00085 // read the first image to get size 00086 img_ptr = new EMData(); 00087 img_ptr->read_image(stackfname, 0); 00088 int nx = img_ptr->get_xsize(); 00089 int ny = img_ptr->get_ysize(); 00090 00091 float s2x, s2y; 00092 00093 if (mypid == 0) { 00094 printf("Master node reading and distributing %d images...\n", nima); 00095 for ( int ip = 0 ; ip < ncpus ; ++ip ) { 00096 for ( int i = 0 ; i < psize[ip] ; ++i ) { 00097 img_index = nbase[ip] + i; 00098 if (ip != 0) { 00099 img_ptr->read_image(stackfname, img_index); 00100 // get a pointer to the image's data 00101 imgdata = img_ptr->get_data(); 00102 // find the x/y shift values if it has them, otherwise set them to 0.0 00103 try { 00104 s2x = (*images2D)[i]->get_attr("s2x"); 00105 } catch ( std::exception& e ) { 00106 s2x = 0.0; 00107 } 00108 try { 00109 s2y = (*images2D)[i]->get_attr("s2y"); 00110 } catch ( std::exception& e ) { 00111 s2y = 0.0; 00112 } 00113 // send these to processor ip 00114 MPI_Send(imgdata, nx*ny, MPI_FLOAT, ip, ip, comm); 00115 MPI_Send(&s2x, 1, MPI_FLOAT, ip, ip, comm); 00116 MPI_Send(&s2y, 1, MPI_FLOAT, ip, ip, comm); 00117 } else { // ip == 0 00118 (*images2D)[i] = new EMData(); 00119 (*images2D)[i]->read_image(stackfname, img_index); 00120 try { 00121 s2x = (*images2D)[i]->get_attr("s2x"); 00122 } catch ( std::exception& e ) { 00123 (*images2D)[i]->set_attr("s2x",0.0); 00124 } 00125 try { 00126 s2y = (*images2D)[i]->get_attr("s2y"); 00127 } catch ( std::exception& e ) { 00128 (*images2D)[i]->set_attr("s2y",0.0); 00129 } 00130 } 00131 } 00132 printf("finished reading data for processor %d\n", ip); 00133 } 00134 } else { // mypid != 0 : everyone else receives and reads in their data 00135 for ( int i = 0 ; i < psize[mypid] ; ++i ) { 00136 (*images2D)[i] = new EMData(); 00137 (*images2D)[i]->set_size(nx, ny, 1); 00138 imgdata = (*images2D)[i]->get_data(); 00139 MPI_Recv(imgdata, nx*ny, MPI_FLOAT, 0, mypid, comm, &mpistatus); 00140 MPI_Recv(&s2x, 1, MPI_FLOAT, 0, mypid, comm, &mpistatus); 00141 MPI_Recv(&s2y, 1, MPI_FLOAT, 0, mypid, comm, &mpistatus); 00142 (*images2D)[i]->set_attr("s2x",s2x); 00143 (*images2D)[i]->set_attr("s2y",s2y); 00144 } 00145 printf("received %d images for processor %d\n", *nloc, mypid); 00146 } 00147 if (mypid == 0) printf("finished reading and distributing data\n"); 00148 ierr = mpierr; 00149 00150 EMDeletePtr(img_ptr); 00151 EMDeleteArray(psize); 00152 EMDeleteArray(nbase); 00153 } 00154 return ierr; 00155 }
|
|
Definition at line 10 of file utilcomm.cpp. References EMAN::EMData::get_data(), EMAN::EMData::get_ndim(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), ierr, nx, ny, EMAN::EMData::read_image(), and EMAN::EMData::set_size(). Referenced by main(). 00011 { 00012 int mypid, ierr, ndim, nx, ny, nz, mpierr=0; 00013 MPI_Comm_rank(comm,&mypid); 00014 FILE *fp=NULL; 00015 00016 ierr = 0; 00017 if (mypid == 0) { 00018 // check to see if the file exists 00019 fp = fopen(volfname,"r"); 00020 if (!fp) { 00021 ierr = 1; 00022 printf("failed to open %s\n", volfname); 00023 } 00024 else { 00025 fclose(fp); 00026 volume->read_image(volfname); 00027 ndim = volume->get_ndim(); 00028 nx = volume->get_xsize(); 00029 ny = volume->get_ysize(); 00030 nz = volume->get_zsize(); 00031 printf("ndim = %d, nx = %d, ny = %d, nz = %d\n", ndim, nx, ny, nz); 00032 } 00033 } 00034 mpierr = MPI_Bcast(&ierr, 1, MPI_INT, 0, comm); 00035 if (ierr==0) { 00036 mpierr = MPI_Bcast(&nx, 1, MPI_INT, 0, comm); 00037 mpierr = MPI_Bcast(&ny, 1, MPI_INT, 0, comm); 00038 mpierr = MPI_Bcast(&nz, 1, MPI_INT, 0, comm); 00039 if (mypid !=0) volume->set_size(nx,ny,nz); 00040 00041 float * voldata = volume->get_data(); 00042 mpierr = MPI_Bcast(voldata, nx*ny*nz, MPI_FLOAT, 0, comm); 00043 ierr = mpierr; 00044 } 00045 return ierr; 00046 }
|
|
Definition at line 279 of file utilcomm.cpp. Referenced by ali3d_d(), fgcalc(), ReadStackandDist(), and unified(). 00280 { 00281 int ncpus, mypid, nangloc, nrem; 00282 00283 MPI_Comm_size(comm,&ncpus); 00284 MPI_Comm_rank(comm,&mypid); 00285 00286 nangloc = nang/ncpus; 00287 nrem = nang - ncpus*nangloc; 00288 if (mypid < nrem) nangloc++; 00289 00290 for (int i = 0; i < ncpus; i++) { 00291 psize[i] = nang/ncpus; 00292 if (i < nrem) psize[i]++; 00293 } 00294 00295 nbase[0] = 0; 00296 for (int i = 1; i < ncpus; i++) { 00297 nbase[i] = nbase[i-1] + psize[i-1]; 00298 } 00299 00300 return nangloc; 00301 }
|