reconstructor.cpp

Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  *
00034  * */
00035 
00036 #include "reconstructor.h"
00037 #include "plugins/reconstructor_template.h"
00038 #include "ctf.h"
00039 #include "emassert.h"
00040 #include "symmetry.h"
00041 #include <cstring>
00042 #include <fstream>
00043 #include <iomanip>
00044 #include <boost/bind.hpp>
00045 #include <boost/format.hpp>
00046 
00047 #include <gsl/gsl_statistics_double.h>
00048 #include <gsl/gsl_fit.h>
00049 
00050 #ifdef EMAN2_USING_CUDA
00051 #include "cuda/cuda_reconstructor.h"
00052 #endif
00053 
00054 using namespace EMAN;
00055 using std::complex;
00056 
00057 
00058 #include <iostream>
00059 using std::cerr;
00060 using std::endl;
00061 using std::cout; // for debug
00062 
00063 #include <algorithm>
00064 // find, for_each
00065 
00066 #include <iomanip>
00067 using std::setprecision;
00068 
00069 template < typename T > void checked_delete( T*& x )
00070 {
00071     typedef char type_must_be_complete[ sizeof(T)? 1: -1 ];
00072     (void) sizeof(type_must_be_complete);
00073     delete x;
00074     x = NULL;
00075 }
00076 
00077 const string FourierReconstructor::NAME = "fourier";
00078 const string FourierReconstructorSimple2D::NAME = "fouriersimple2D";
00079 const string WienerFourierReconstructor::NAME = "wiener_fourier";
00080 const string BackProjectionReconstructor::NAME = "back_projection";
00081 const string nn4Reconstructor::NAME = "nn4";
00082 const string nn4_rectReconstructor::NAME = "nn4_rect";
00083 const string nnSSNR_Reconstructor::NAME = "nnSSNR";
00084 const string nn4_ctfReconstructor::NAME = "nn4_ctf";
00085 const string nn4_ctf_rectReconstructor::NAME = "nn4_ctf_rect";
00086 const string nnSSNR_ctfReconstructor::NAME = "nnSSNR_ctf";
00087 
00088 template <> Factory < Reconstructor >::Factory()
00089 {
00090         force_add<FourierReconstructor>();
00091         force_add<FourierReconstructorSimple2D>();
00092 //      force_add(&BaldwinWoolfordReconstructor::NEW);
00093         force_add<WienerFourierReconstructor>();
00094         force_add<BackProjectionReconstructor>();
00095         force_add<nn4Reconstructor>();
00096         force_add<nn4_rectReconstructor>();
00097         force_add<nnSSNR_Reconstructor>();
00098         force_add<nn4_ctfReconstructor>();
00099         force_add<nn4_ctf_rectReconstructor>();
00100         force_add<nnSSNR_ctfReconstructor>();
00101 //      force_add<XYZReconstructor>();
00102 }
00103 
00104 void FourierReconstructorSimple2D::setup()
00105 {
00106         nx = params.set_default("nx",0);
00107 
00108         if ( nx < 0 ) throw InvalidValueException(nx, "nx must be positive");
00109 
00110         bool is_fftodd = (nx % 2 == 1);
00111 
00112         ny = nx;
00113         nx += 2-is_fftodd;
00114 
00115         image = new EMData();
00116         image->set_size(nx, ny);
00117         image->set_complex(true);
00118         image->set_fftodd(is_fftodd);
00119         image->set_ri(true);
00120 
00121         tmp_data = new EMData();
00122         tmp_data->set_size(nx/2, nx);
00123 }
00124 
00125 int FourierReconstructorSimple2D::insert_slice(const EMData* const slice, const Transform & euler, const float)
00126 {
00127 
00128         // Are these exceptions really necessary? (d.woolford)
00129         if (!slice) throw NullPointerException("EMData pointer (input image) is NULL");
00130 
00131         if ( slice->get_ndim() != 1 ) throw ImageDimensionException("Image dimension must be 1");
00132 
00133         // I could also test to make sure the image is the right dimensions...
00134         if (slice->is_complex()) throw ImageFormatException("The image is complex, expecting real");
00135 
00136         EMData* working_slice = slice->process("xform.phaseorigin.tocorner");
00137 
00138         // Fourier transform the slice
00139         working_slice->do_fft_inplace();
00140 
00141         float* rdata = image->get_data();
00142         float* norm = tmp_data->get_data();
00143         float* dat = working_slice->get_data();
00144 
00145         float g[4];
00146         int offset[4];
00147         float dt[2];
00148         offset[0] = 0; offset[1] = 2; offset[2] = nx; offset[3] = nx+2;
00149 
00150         float alt = -((float)(euler.get_rotation("2d"))["alpha"])*M_PI/180.0f;
00151         for (int x = 0; x < working_slice->get_xsize() / 2; x++) {
00152 
00153                 float rx = (float) x;
00154 
00155                 float xx = rx*cos(alt);
00156                 float yy = rx*sin(alt);
00157                 float cc = 1.0;
00158 
00159                 if (xx < 0) {
00160                         xx = -xx;
00161                         yy = -yy;
00162                         cc = -1.0;
00163                 }
00164 
00165                 yy += ny / 2;
00166 
00167 
00168                 dt[0] = dat[2*x];
00169                 dt[1] = cc * dat[2*x+1];
00170 
00171                 // PHIL IS INTERESTED FROM HERE DOWN
00172                 int x0 = (int) floor(xx);
00173                 int y0 = (int) floor(yy);
00174 
00175                 int i = 2*x0 + y0*nx;
00176 
00177                 float dx = xx - x0;
00178                 float dy = yy - y0;
00179 
00180                 g[0] = Util::agauss(1, dx, dy, 0, EMConsts::I2G);
00181                 g[1] = Util::agauss(1, 1 - dx, dy, 0, EMConsts::I2G);
00182                 g[2] = Util::agauss(1, dx, 1 - dy, 0, EMConsts::I2G);
00183                 g[3] = Util::agauss(1, 1 - dx, 1 - dy, 0, EMConsts::I2G);
00184 
00185                 // At the extreme we can only do some much...
00186                 if ( x0 == nx-2 ) {
00187                         int k = i + offset[0];
00188                         rdata[k] += g[0] * dt[0];
00189                         rdata[k + 1] += g[0] * dt[1];
00190                         norm[k/2] += g[0];
00191 
00192                         k = i + offset[2];
00193                         rdata[k] += g[2] * dt[0];
00194                         rdata[k + 1] += g[2] * dt[1];
00195                         norm[k/2] += g[2];
00196                         continue;
00197 
00198                 }
00199                 // capture and accommodate for periodic boundary conditions in the x direction
00200                 if ( x0 > nx-2 ) {
00201                         int dif = x0 - (nx-2);
00202                         x0 -= dif;
00203                 }
00204                 // At the extreme we can only do some much...
00205                 if ( y0 == ny -1 ) {
00206                         int k = i + offset[0];
00207                         rdata[k] += g[0] * dt[0];
00208                         rdata[k + 1] += g[0] * dt[1];
00209                         norm[k/2] += g[0];
00210 
00211                         k = i + offset[1];
00212                         rdata[k] += g[1] * dt[0];
00213                         rdata[k + 1] += g[1] * dt[1];
00214                         norm[k/2] += g[1];
00215                         continue;
00216                 }
00217                 // capture and accommodate for periodic boundary conditions in the y direction
00218                 if ( y0 > ny-1) {
00219                         int dif = y0 - (ny-1);
00220                         y0 -= dif;
00221                 }
00222 
00223                 if (x0 >= nx - 2 || y0 >= ny - 1) continue;
00224 
00225 
00226 
00227 
00228                 for (int j = 0; j < 4; j++)
00229                 {
00230                         int k = i + offset[j];
00231                         rdata[k] += g[j] * dt[0];
00232                         rdata[k + 1] += g[j] * dt[1];
00233                         norm[k/2] += g[j];
00234 
00235                 }
00236         }
00237 
00238         return 0;
00239 
00240 }
00241 
00242 EMData *FourierReconstructorSimple2D::finish(bool)
00243 {
00244         normalize_threed();
00245 
00246         image->process_inplace("xform.fourierorigin.tocorner");
00247         image->do_ift_inplace();
00248         image->depad();
00249         image->process_inplace("xform.phaseorigin.tocenter");
00250 
00251         EMData *ret = image;
00252         image = 0;
00253         return  ret;
00254 }
00255 
00256 void ReconstructorVolumeData::normalize_threed(const bool sqrt_damp,const bool wiener)
00257 // normalizes the 3-D Fourier volume. Also imposes appropriate complex conjugate relationships
00258 {
00259         float* norm = tmp_data->get_data();
00260         float* rdata = image->get_data();
00261 
00262 //      printf("%d %d %d %d %d %d\n",subnx,subny,subnz,image->get_xsize(),image->get_ysize(),image->get_zsize());
00263 
00264         // FIXME should throw a sensible error
00265         if ( 0 == norm ) throw NullPointerException("The normalization volume was null!");
00266         if ( 0 == rdata ) throw NullPointerException("The complex reconstruction volume was null!");
00267 
00268         // The math is a bit tricky to explain. Wiener filter is basically SNR/(1+SNR)
00269         // In this case, data have already been effectively multiplied by SNR (one image at a time),
00270         // so without Wiener filter, we just divide by total SNR. With Wiener filter, we just add
00271         // 1.0 to total SNR, and we're done             --steve
00272         float wfilt=0.0;
00273         if (wiener) wfilt=1.0;          
00274         
00275         for (size_t i = 0; i < (size_t)subnx * subny * subnz; i += 2) {
00276                 float d = norm[i/2];
00277                 if (sqrt_damp) d*=sqrt(d);
00278                 if (d == 0) {
00279                         rdata[i] = 0;
00280                         rdata[i + 1] = 0;
00281                 }
00282                 else {
00283 //                      rdata[i]=1.0/d;
00284 //                      rdata[i+1]=0.0;         // for debugging only
00285                         rdata[i] /= d+wfilt;
00286                         rdata[i + 1] /= d+wfilt;
00287                 }
00288         }
00289 
00290         // enforce complex conj, only works on subvolumes if the complete conjugate plane is in the volume
00291         if (subx0==0 && subnx>1 && subny==ny && subnz==nz) {
00292                 for (int z=0; z<=nz/2; z++) {
00293                         for (int y=1; y<=ny; y++) {
00294                                 if (y==0 && z==0) continue;
00295                                 // x=0
00296                                 size_t i=(size_t)(y%ny)*subnx+(size_t)(z%nz)*subnx*subny;
00297                                 size_t i2=(size_t)(ny-y)*subnx+(size_t)((nz-z)%nz)*subnx*subny;
00298                                 float ar=(rdata[i]+rdata[i2])/2.0f;
00299                                 float ai=(rdata[i+1]-rdata[i2+1])/2.0f;
00300                                 rdata[i]=ar;
00301                                 rdata[i2]=ar;
00302                                 rdata[i+1]=ai;
00303                                 rdata[i2+1]=-ai;
00304                         }
00305                 }
00306         }
00307 
00308         if (subx0+subnx==nx && subnx>1 && subny==ny && subnz==nz) {
00309                 for (int z=0; z<=nz/2; z++) {
00310                         for (int y=1; y<=ny; y++) {
00311                                 if (y==0 && z==0) continue;
00312                                 // x=0
00313                                 size_t i=(size_t)(y%ny)*subnx+(size_t)(z%nz)*subnx*subny+subnx-2;
00314                                 size_t i2=(size_t)(ny-y)*subnx+(size_t)((nz-z)%nz)*subnx*subny+subnx-2;
00315                                 float ar=(rdata[i]+rdata[i2])/2.0f;
00316                                 float ai=(rdata[i+1]-rdata[i2+1])/2.0f;
00317                                 rdata[i]=ar;
00318                                 rdata[i2]=ar;
00319                                 rdata[i+1]=ai;
00320                                 rdata[i2+1]=-ai;
00321                         }
00322                 }
00323         }
00324 }
00325 
00326 void FourierReconstructor::load_default_settings()
00327 {
00328         inserter=0;
00329         image=0;
00330         tmp_data=0;
00331 }
00332 
00333 void FourierReconstructor::free_memory()
00334 {
00335         if (image) { delete image; image=0; }
00336         if (tmp_data) { delete tmp_data; tmp_data=0; }
00337         if ( inserter != 0 )
00338         {
00339                 delete inserter;
00340                 inserter = 0;
00341         }
00342 }
00343 
00344 #include <sstream>
00345 
00346 void FourierReconstructor::load_inserter()
00347 {
00348 //      ss
00349 //      string mode = (string)params["mode"];
00350         Dict parms;
00351         parms["data"] = image;
00352         parms["norm"] = tmp_data->get_data();
00353         // These aren't necessary because we deal with them before calling the inserter
00354 //      parms["subnx"] = nx;
00355 //      parms["subny"] = ny;
00356 //      parms["subnx"] = nz;
00357 //      parms["subx0"] = x0;
00358 //      parms["suby0"] = y0;
00359 //      parms["subz0"] = z0;
00360 
00361         if ( inserter != 0 )
00362         {
00363                 delete inserter;
00364         }
00365 
00366         inserter = Factory<FourierPixelInserter3D>::get((string)params["mode"], parms);
00367         inserter->init();
00368 }
00369 
00370 void FourierReconstructor::setup()
00371 {
00372         // default setting behavior - does not override if the parameter is already set
00373         params.set_default("mode","gauss_2");
00374 
00375         vector<int> size=params["size"];
00376 
00377         nx = size[0];
00378         ny = size[1];
00379         nz = size[2];
00380         nx2=nx/2-1;
00381         ny2=ny/2;
00382         nz2=nz/2;
00383 
00384 
00385         // Adjust nx if for Fourier transform even odd issues
00386         bool is_fftodd = (nx % 2 == 1);
00387         // The Fourier transform requires one extra pixel in the x direction,
00388         // which is two spaces in memory, one each for the complex and the
00389         // real components
00390         nx += 2-is_fftodd;
00391 
00392         if (params.has_key("subvolume")) {
00393                 vector<int> sub=params["subvolume"];
00394                 subx0=sub[0];
00395                 suby0=sub[1];
00396                 subz0=sub[2];
00397                 subnx=sub[3];
00398                 subny=sub[4];
00399                 subnz=sub[5];
00400 
00401                 if (subx0<0 || suby0<0 || subz0<0 || subx0+subnx>nx || suby0+subny>ny || subz0+subnz>nz)
00402                         throw ImageDimensionException("The subvolume cannot extend outside the reconstructed volume");
00403 
00404         }
00405         else {
00406                 subx0=suby0=subz0=0;
00407                 subnx=nx;
00408                 subny=ny;
00409                 subnz=nz;
00410         }
00411 
00412 
00413         // Odd dimension support is here atm, but not above.
00414         if (image) delete image;
00415         image = new EMData();
00416         image->set_size(subnx, subny, subnz);
00417         image->set_complex(true);
00418         image->set_fftodd(is_fftodd);
00419         image->set_ri(true);
00420 //      printf("%d %d %d\n\n",subnx,subny,subnz);
00421         image->to_zero();
00422 
00423         if (params.has_key("subvolume")) {
00424                 image->set_attr("subvolume_x0",subx0);
00425                 image->set_attr("subvolume_y0",suby0);
00426                 image->set_attr("subvolume_z0",subz0);
00427                 image->set_attr("subvolume_full_nx",nx);
00428                 image->set_attr("subvolume_full_ny",ny);
00429                 image->set_attr("subvolume_full_nz",nz);
00430         }
00431         
00432         if (tmp_data) delete tmp_data;
00433         tmp_data = new EMData();
00434         tmp_data->set_size(subnx/2, subny, subnz);
00435         tmp_data->to_zero();
00436         tmp_data->update();
00437 
00438         load_inserter();
00439 
00440         if ( (bool) params["quiet"] == false )
00441         {
00442                 cout << "3D Fourier dimensions are " << nx << " " << ny << " " << nz << endl;
00443                 cout << "3D Fourier subvolume is " << subnx << " " << subny << " " << subnz << endl;
00444                 printf ("You will require approximately %1.3g GB of memory to reconstruct this volume\n",((float)subnx*subny*subnz*sizeof(float)*1.5)/1000000000.0);
00445         }
00446 }
00447 
00448 void FourierReconstructor::setup_seed(EMData* seed,float seed_weight) {
00449         // default setting behavior - does not override if the parameter is already set
00450         params.set_default("mode","gauss_2");
00451 
00452         vector<int> size=params["size"];
00453 
00454         nx = size[0];
00455         ny = size[1];
00456         nz = size[2];
00457         nx2=nx/2-1;
00458         ny2=ny/2;
00459         nz2=nz/2;
00460 
00461 
00462         // Adjust nx if for Fourier transform even odd issues
00463         bool is_fftodd = (nx % 2 == 1);
00464         // The Fourier transform requires one extra pixel in the x direction,
00465         // which is two spaces in memory, one each for the complex and the
00466         // real components
00467         nx += 2-is_fftodd;
00468 
00469         if (params.has_key("subvolume")) {
00470                 vector<int> sub=params["subvolume"];
00471                 subx0=sub[0];
00472                 suby0=sub[1];
00473                 subz0=sub[2];
00474                 subnx=sub[3];
00475                 subny=sub[4];
00476                 subnz=sub[5];
00477 
00478                 if (subx0<0 || suby0<0 || subz0<0 || subx0+subnx>nx || suby0+subny>ny || subz0+subnz>nz)
00479                         throw ImageDimensionException("The subvolume cannot extend outside the reconstructed volume");
00480 
00481         }
00482         else {
00483                 subx0=suby0=subz0=0;
00484                 subnx=nx;
00485                 subny=ny;
00486                 subnz=nz;
00487         }
00488 
00489         if (seed->get_xsize()!=subnx || seed->get_ysize()!=subny || seed->get_zsize()!=subnz || !seed->is_complex())
00490                 throw ImageDimensionException("The dimensions of the seed volume do not match the reconstruction size");
00491 
00492         // Odd dimension support is here atm, but not above.
00493         image = seed;
00494         if (params.has_key("subvolume")) {
00495                 image->set_attr("subvolume_x0",subx0);
00496                 image->set_attr("subvolume_y0",suby0);
00497                 image->set_attr("subvolume_z0",subz0);
00498                 image->set_attr("subvolume_full_nx",nx);
00499                 image->set_attr("subvolume_full_ny",ny);
00500                 image->set_attr("subvolume_full_nz",nz);
00501         }
00502 
00503         if (tmp_data) delete tmp_data;
00504         tmp_data = new EMData();
00505         tmp_data->set_size(subnx/2, subny, subnz);
00506         tmp_data->to_value(seed_weight);
00507 
00508         load_inserter();
00509 
00510         if ( (bool) params["quiet"] == false )
00511         {
00512                 cout << "Seeded direct Fourier inversion";
00513                 cout << "3D Fourier dimensions are " << nx << " " << ny << " " << nz << endl;
00514                 cout << "3D Fourier subvolume is " << subnx << " " << subny << " " << subnz << endl;
00515                 cout << "You will require approximately " << setprecision(3) << (subnx*subny*subnz*sizeof(float)*1.5)/1000000000.0 << "GB of memory to reconstruct this volume" << endl;
00516         }
00517 }
00518 
00519 void FourierReconstructor::clear()
00520 {
00521         bool zeroimage = true;
00522         bool zerotmpimg = true;
00523         
00524 #ifdef EMAN2_USING_CUDA
00525         if(EMData::usecuda == 1) {
00526                 if(image->getcudarwdata()) {
00527                         to_zero_cuda(image->getcudarwdata(),image->get_xsize(),image->get_ysize(),image->get_zsize());
00528                         zeroimage = false;
00529                 }
00530                 if(tmp_data->getcudarwdata()) {
00531                         //to_zero_cuda(tmp_data->getcudarwdata(),image->get_xsize(),image->get_ysize(),image->get_zsize());
00532                         zerotmpimg = false;
00533                 }
00534         }
00535 #endif
00536 
00537         if(zeroimage) image->to_zero();
00538         if(zerotmpimg) tmp_data->to_zero();
00539         
00540 }
00541 
00542 EMData* FourierReconstructor::preprocess_slice( const EMData* const slice,  const Transform& t )
00543 {
00544 #ifdef EMAN2_USING_CUDA
00545         if(EMData::usecuda == 1) {
00546                 if(!slice->getcudarwdata()) slice->copy_to_cuda_keepcpu(); //copy slice to cuda using the const version
00547         }
00548 #endif
00549         // Shift the image pixels so the real space origin is now located at the phase origin (at the bottom left of the image)
00550         EMData* return_slice = 0;
00551         Transform tmp(t);
00552         tmp.set_rotation(Dict("type","eman")); // resets the rotation to 0 implicitly, this way we only do 2d translation,scaling and mirroring
00553 
00554         if (tmp.is_identity()) return_slice=slice->copy();
00555         else return_slice = slice->process("xform",Dict("transform",&tmp));
00556 
00557         return_slice->process_inplace("xform.phaseorigin.tocorner");
00558 
00559 //      printf("@@@ %d\n",(int)return_slice->get_attr("nx"));
00560         // Fourier transform the slice
00561         
00562 #ifdef EMAN2_USING_CUDA
00563         if(return_slice->getcudarwdata()) {
00564                 return_slice->do_fft_inplace_cuda(); //a CUDA FFT inplace is quite slow as there is a lot of mem copying.
00565         }else{
00566                 return_slice->do_fft_inplace();
00567         }
00568 #else
00569         return_slice->do_fft_inplace();
00570 #endif
00571 
00572 //      printf("%d\n",(int)return_slice->get_attr("nx"));
00573 
00574         return_slice->mult((float)sqrt(1.0f/(return_slice->get_ysize())*return_slice->get_xsize()));
00575 
00576         // Shift the Fourier transform so that it's origin is in the center (bottom) of the image.
00577 //      return_slice->process_inplace("xform.fourierorigin.tocenter");
00578 
00579         return_slice->set_attr("reconstruct_preproc",(int)1);
00580         return return_slice;
00581 }
00582 
00583 
00584 int FourierReconstructor::insert_slice(const EMData* const input_slice, const Transform & arg, const float weight)
00585 {
00586         // Are these exceptions really necessary? (d.woolford)
00587         if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
00588 
00589 #ifdef EMAN2_USING_CUDA
00590         if(EMData::usecuda == 1) {
00591                 if(!input_slice->getcudarwdata()) input_slice->copy_to_cuda_keepcpu(); //copy slice to cuda using the const version
00592         }
00593 #endif
00594 
00595         Transform * rotation;
00596 /*      if ( input_slice->has_attr("xform.projection") ) {
00597                 rotation = (Transform*) (input_slice->get_attr("xform.projection")); // assignment operator
00598         } else {*/
00599         rotation = new Transform(arg); // assignment operator
00600 //      }
00601 
00602         EMData *slice;
00603         if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
00604         else slice = preprocess_slice( input_slice, *rotation);
00605 
00606 
00607         // We must use only the rotational component of the transform, scaling, translation and mirroring
00608         // are not implemented in Fourier space, but are in preprocess_slice
00609         rotation->set_scale(1.0);
00610         rotation->set_mirror(false);
00611         rotation->set_trans(0,0,0);
00612 
00613         // Finally to the pixel wise slice insertion
00614         //slice->copy_to_cuda();
00615         do_insert_slice_work(slice, *rotation, weight);
00616         
00617         delete rotation; rotation=0;
00618         delete slice;
00619 
00620 //      image->update();
00621         return 0;
00622 }
00623 
00624 void FourierReconstructor::do_insert_slice_work(const EMData* const input_slice, const Transform & arg,const float weight)
00625 {
00626         // Reload the inserter if the mode has changed
00627 //      string mode = (string) params["mode"];
00628 //      if ( mode != inserter->get_name() )     load_inserter();
00629 
00630 //      int y_in = input_slice->get_ysize();
00631 //      int x_in = input_slice->get_xsize();
00632 //      // Adjust the dimensions to account for odd and even ffts
00633 //      if (input_slice->is_fftodd()) x_in -= 1;
00634 //      else x_in -= 2;
00635 
00636         vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
00637 
00638         float inx=(float)(input_slice->get_xsize());            // x/y dimensions of the input image
00639         float iny=(float)(input_slice->get_ysize());
00640 
00641 #ifdef EMAN2_USING_CUDA
00642         if(EMData::usecuda == 1) {
00643                 if(!image->getcudarwdata()){
00644                         image->copy_to_cuda();
00645                         tmp_data->copy_to_cuda();
00646                 }
00647                 float * m = new float[12];
00648                 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
00649                         Transform t3d = arg*(*it);
00650                         t3d.copy_matrix_into_array(m);
00651                         //cout << "using CUDA " << image->getcudarwdata() << endl;
00652                         insert_slice_cuda(m,input_slice->getcudarwdata(),image->getcudarwdata(),tmp_data->getcudarwdata(),inx,iny,image->get_xsize(),image->get_ysize(),image->get_zsize(), weight);
00653                 }
00654                 delete m;
00655                 return;
00656         }
00657 #endif
00658         for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
00659                 Transform t3d = arg*(*it);
00660                 for (int y = -iny/2; y < iny/2; y++) {
00661                         for (int x = 0; x <=  inx/2; x++) {
00662 
00663                                 float rx = (float) x/(inx-2.0f);        // coords relative to Nyquist=.5
00664                                 float ry = (float) y/iny;
00665 
00666                                 Vec3f coord(rx,ry,0);
00667                                 coord = coord*t3d; // transpose multiplication
00668                                 float xx = coord[0]; // transformed coordinates in terms of Nyquist
00669                                 float yy = coord[1];
00670                                 float zz = coord[2];
00671 
00672                                 // Map back to real pixel coordinates in output volume
00673                                 xx=xx*(nx-2);
00674                                 yy=yy*ny;
00675                                 zz=zz*nz;
00676 
00677 //                              if (x==10 && y==0) printf("10,0 -> %1.2f,%1.2f,%1.2f\t(%5.2f %5.2f %5.2f   %5.2f %5.2f %5.2f   %5.2f %5.2f %5.2f) %1.0f %d\n",
00678 //                                      xx,yy,zz,t3d.at(0,0),t3d.at(0,1),t3d.at(0,2),t3d.at(1,0),t3d.at(1,1),t3d.at(1,2),t3d.at(2,0),t3d.at(2,1),t3d.at(2,2),inx,nx);
00679 //                              if (x==0 && y==10 FourierReconstructor:) printf("0,10 -> %1.2f,%1.2f,%1.2f\t(%5.2f %5.2f %5.2f   %5.2f %5.2f %5.2f   %5.2f %5.2f %5.2f)\n",
00680 //                                      xx,yy,zz,t3d.at(0,0),t3d.at(0,1),t3d.at(0,2),t3d.at(1,0),t3d.at(1,1),t3d.at(1,2),t3d.at(2,0),t3d.at(2,1),t3d.at(2,2));
00681 
00682                                 //printf("%3.1f %3.1f %3.1f\t %1.4f %1.4f\t%1.4f\n",xx,yy,zz,input_slice->get_complex_at(x,y).real(),input_slice->get_complex_at(x,y).imag(),weight);
00683 //                              if (floor(xx)==45 && floor(yy)==45 &&floor(zz)==0) printf("%d. 45 45 0\t %d %d\t %1.4f %1.4f\t%1.4f\n",(int)input_slice->get_attr("n"),x,y,input_slice->get_complex_at(x,y).real(),input_slice->get_complex_at(x,y).imag(),weight);
00684 //                              if (floor(xx)==21 && floor(yy)==21 &&floor(zz)==0) printf("%d. 21 21 0\t %d %d\t %1.4f %1.4f\t%1.4f\n",(int)input_slice->get_attr("n"),x,y,input_slice->get_complex_at(x,y).real(),input_slice->get_complex_at(x,y).imag(),weight);
00685                                 inserter->insert_pixel(xx,yy,zz,input_slice->get_complex_at(x,y),weight);
00686                         }
00687                 }
00688         }
00689 }
00690 
00691 int FourierReconstructor::determine_slice_agreement(EMData*  input_slice, const Transform & arg, const float weight,bool sub)
00692 {
00693         // Are these exceptions really necessary? (d.woolford)
00694         if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
00695 
00696 #ifdef EMAN2_USING_CUDA
00697         if(EMData::usecuda == 1) {
00698                 if(!input_slice->getcudarwdata()) input_slice->copy_to_cuda(); //copy slice to cuda using the const version
00699         }
00700 #endif
00701 
00702         Transform * rotation;
00703         rotation = new Transform(arg); // assignment operator
00704 
00705         EMData *slice;
00706         if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
00707         else slice = preprocess_slice( input_slice, *rotation);
00708 
00709 
00710         // We must use only the rotational component of the transform, scaling, translation and mirroring
00711         // are not implemented in Fourier space, but are in preprocess_slice
00712         rotation->set_scale(1.0);
00713         rotation->set_mirror(false);
00714         rotation->set_trans(0,0,0);
00715         if (sub) do_insert_slice_work(slice, *rotation, -weight);
00716         // Remove the current slice first (not threadsafe, but otherwise performance would be awful)
00717         
00718         // Compare
00719         do_compare_slice_work(slice, *rotation,weight);
00720 
00721         input_slice->set_attr("reconstruct_norm",slice->get_attr("reconstruct_norm"));
00722         input_slice->set_attr("reconstruct_absqual",slice->get_attr("reconstruct_absqual"));
00723 //      input_slice->set_attr("reconstruct_qual",slice->get_attr("reconstruct_qual"));
00724         input_slice->set_attr("reconstruct_weight",slice->get_attr("reconstruct_weight"));
00725 
00726         // Now put the slice back
00727         if (sub) do_insert_slice_work(slice, *rotation, weight);
00728 
00729         delete rotation;
00730         delete slice;
00731 
00732 //      image->update();
00733         return 0;
00734 
00735 }
00736 
00737 void FourierReconstructor::do_compare_slice_work(EMData* input_slice, const Transform & arg,float weight)
00738 {
00739 
00740         float dt[3];    // This stores the complex and weight from the volume
00741         float dt2[2];   // This stores the local image complex
00742         float *dat = input_slice->get_data();
00743         vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
00744 
00745         float inx=(float)(input_slice->get_xsize());            // x/y dimensions of the input image
00746         float iny=(float)(input_slice->get_ysize());
00747 
00748         double dot=0;           // summed pixel*weight dot product
00749         double vweight=0;               // sum of weights
00750         double power=0;         // sum of inten*weight from volume
00751         double power2=0;                // sum of inten*weight from image
00752         bool use_cpu = true;
00753         
00754 #ifdef EMAN2_USING_CUDA
00755         if(EMData::usecuda == 1) {
00756                 if(!input_slice->getcudarwdata()) input_slice->copy_to_cuda();
00757                 if(!image->getcudarwdata()){
00758                         image->copy_to_cuda();
00759                         tmp_data->copy_to_cuda();
00760                 }
00761                 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
00762                         Transform t3d = arg*(*it);
00763                         float * m = new float[12];
00764                         t3d.copy_matrix_into_array(m);
00765                         float4 stats = determine_slice_agreement_cuda(m,input_slice->getcudarwdata(),image->getcudarwdata(),tmp_data->getcudarwdata(),inx,iny,image->get_xsize(),image->get_ysize(),image->get_zsize(), weight);
00766                         dot = stats.x;
00767                         vweight = stats.y;
00768                         power = stats.z;
00769                         power2 = stats.w;
00770                         //cout << "CUDA stats " << stats.x << " " << stats.y << " " << stats.z << " " << stats.w << endl;
00771                         use_cpu = false;
00772                 }
00773         }
00774 #endif
00775         if(use_cpu) {
00776                 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
00777                         Transform t3d = arg*(*it);
00778                         for (int y = -iny/2; y < iny/2; y++) {
00779                                 for (int x = 0; x <=  inx/2; x++) {
00780                                         if (x==0 && y==0) continue;             // We don't want to use the Fourier origin
00781 
00782                                         float rx = (float) x/(inx-2);   // coords relative to Nyquist=.5
00783                                         float ry = (float) y/iny;
00784 
00785 //                                      if ((rx * rx + Util::square(ry - max_input_dim "xform.projection"/ 2)) > rl)
00786 //                                      continue;
00787 
00788                                         Vec3f coord(rx,ry,0);
00789                                         coord = coord*t3d; // transpose multiplication
00790                                         float xx = coord[0]; // transformed coordinates in terms of Nyquist
00791                                         float yy = coord[1];
00792                                         float zz = coord[2];
00793 
00794 
00795                                         if (fabs(xx)>0.5 || fabs(yy)>=0.5 || fabs(zz)>=0.5) continue;
00796 
00797                                         // Map back to actual pixel coordinates in output volume
00798                                         xx=xx*(nx-2);
00799                                         yy=yy*ny;
00800                                         zz=zz*nz;
00801 
00802 
00803                                         int idx = (int)(x * 2 + inx*(y<0?iny+y:y));
00804                                         dt2[0] = dat[idx];
00805                                         dt2[1] = dat[idx+1];
00806 
00807                                         // value returned indirectly in dt
00808                                         if (!pixel_at(xx,yy,zz,dt) || dt[2]==0) continue;
00809 
00810 //                                      printf("%f\t%f\t%f\t%f\t%f\n",dt[0],dt[1],dt[2],dt2[0],dt2[1]);
00811                                         dot+=(dt[0]*dt2[0]+dt[1]*dt2[1])*dt[2];
00812                                         vweight+=dt[2];
00813                                         power+=(dt[0]*dt[0]+dt[1]*dt[1])*dt[2];
00814                                         power2+=(dt2[0]*dt2[0]+dt2[1]*dt2[1])*dt[2];
00815                                 }
00816                         }
00817                         //cout << dot << " " << vweight << " " << power << " " << power2 << endl;
00818                 }
00819         }
00820         
00821         dot/=sqrt(power*power2);                // normalize the dot product
00822 //      input_slice->set_attr("reconstruct_norm",(float)(power2<=0?1.0:sqrt(power/power2)/(inx*iny)));
00823         input_slice->set_attr("reconstruct_norm",(float)(power2<=0?1.0:sqrt(power/power2)));
00824         input_slice->set_attr("reconstruct_absqual",(float)dot);
00825         float rw=weight<=0?1.0f:1.0f/weight;
00826         input_slice->set_attr("reconstruct_qual",(float)(dot*rw/((rw-1.0)*dot+1.0)));   // here weight is a proxy for SNR
00827         input_slice->set_attr("reconstruct_weight",(float)vweight/(float)(subnx*subny*subnz));
00828 //      printf("** %g\t%g\t%g\t%g ##\n",dot,vweight,power,power2);
00829         //printf("** %f %f %f ##\n",(float)(power2<=0?1.0:sqrt(power/power2)/(inx*iny)),(float)dot,(float)(dot*weight/((weight-1.0)*dot+1.0)));
00830 }
00831 
00832 bool FourierReconstructor::pixel_at(const float& xx, const float& yy, const float& zz, float *dt)
00833 {
00834         int x0 = (int) floor(xx);
00835         int y0 = (int) floor(yy);
00836         int z0 = (int) floor(zz);
00837         
00838         float *rdata=image->get_data();
00839         float *norm=tmp_data->get_data();
00840         float normsum=0,normsum2=0;
00841 
00842         dt[0]=dt[1]=dt[2]=0.0;
00843 
00844         if (nx==subnx) {                        // normal full reconstruction
00845                 if (x0<-nx2-1 || y0<-ny2-1 || z0<-nz2-1 || x0>nx2 || y0>ny2 || z0>nz2 ) return false;
00846 
00847                 // no error checking on add_complex_fast, so we need to be careful here
00848                 int x1=x0+1;
00849                 int y1=y0+1;
00850                 int z1=z0+1;
00851                 if (x0<-nx2) x0=-nx2;
00852                 if (x1>nx2) x1=nx2;
00853                 if (y0<-ny2) y0=-ny2;
00854                 if (y1>ny2) y1=ny2;
00855                 if (z0<-nz2) z0=-nz2;
00856                 if (z1>nz2) z1=nz2;
00857                 
00858                 size_t idx=0;
00859                 float r, gg;
00860                 for (int k = z0 ; k <= z1; k++) {
00861                         for (int j = y0 ; j <= y1; j++) {
00862                                 for (int i = x0; i <= x1; i ++) {
00863                                         r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
00864                                         idx=image->get_complex_index_fast(i,j,k);
00865                                         gg = Util::fast_exp(-r / EMConsts::I2G);
00866                                         
00867                                         dt[0]+=gg*rdata[idx];
00868                                         dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1];
00869                                         dt[2]+=norm[idx/2]*gg;
00870                                         normsum2+=gg;
00871                                         normsum+=gg*norm[idx/2];                                
00872                                 }
00873                         }
00874                 }
00875                 if (normsum==0) return false;
00876                 dt[0]/=normsum;
00877                 dt[1]/=normsum;
00878                 dt[2]/=normsum2;
00879 //              printf("%1.2f,%1.2f,%1.2f\t%1.3f\t%1.3f\t%1.3f\t%1.3f\t%1.3f\n",xx,yy,zz,dt[0],dt[1],dt[2],rdata[idx],rdata[idx+1]);
00880                 return true;
00881         } 
00882         else {                                  // for subvolumes, not optimized yet
00883                 size_t idx;
00884                 float r, gg;
00885                 for (int k = z0 ; k <= z0 + 1; k++) {
00886                         for (int j = y0 ; j <= y0 + 1; j++) {
00887                                 for (int i = x0; i <= x0 + 1; i ++) {
00888                                         r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
00889                                         idx=image->get_complex_index(i,j,k,subx0,suby0,subz0,nx,ny,nz);
00890                                         gg = Util::fast_exp(-r / EMConsts::I2G)*norm[idx/2];
00891                                         
00892                                         dt[0]+=gg*rdata[idx];
00893                                         dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1];
00894                                         dt[2]+=norm[idx/2];
00895                                         normsum+=gg;                            
00896                                 }
00897                         }
00898                 }
00899                 
00900                 if (normsum==0)  return false;
00901                 return true;
00902         }
00903 }
00904 
00905 
00906 EMData *FourierReconstructor::finish(bool doift)
00907 {
00908 //      float *norm = tmp_data->get_data();
00909 //      float *rdata = image->get_data();
00910 #ifdef EMAN2_USING_CUDA
00911         if(image->getcudarwdata()){
00912                 cout << "copy back from CUDA" << endl;
00913                 image->copy_from_device();
00914                 tmp_data->copy_from_device();
00915         }
00916 #endif
00917         
00918         bool sqrtnorm=params.set_default("sqrtnorm",false);
00919         normalize_threed(sqrtnorm);
00920         
00921 //      tmp_data->write_image("density.mrc");
00922 
00923         // we may as well delete the tmp data now... it saves memory and the calling program might
00924         // need memory after it gets the return volume.
00925         // If this delete didn't happen now, it would happen when the deconstructor was called --david
00926         // no longer a good idea with the new iterative scheme -- steve
00927 //      if ( tmp_data != 0 )
00928 //      {
00929 //              delete tmp_data;
00930 //              tmp_data = 0;
00931 //      }
00932 
00933 /*      image->process_inplace("xform.fourierorigin.tocorner");*/
00934 
00935         if (doift) {
00936                 image->do_ift_inplace();
00937                 image->depad();
00938                 image->process_inplace("xform.phaseorigin.tocenter");
00939         }
00940         // If the image was padded it should be the original size, as the client would expect
00941         //  I blocked the rest, it is almost certainly incorrect  PAP 07/31/08
00942         // No, it's not incorrect. You are wrong. You have the meaning of nx mixed up. DSAW 09/23/cd
00943         // This should now be handled in the calling program --steve 11/03/09
00944 //      bool is_fftodd = (nx % 2 == 1);
00945 //      if ( (nx-2*(!is_fftodd)) != output_x || ny != output_y || nz != output_z )
00946 //      {
00947 //              FloatPoint origin( (nx-output_x)/2, (ny-output_y)/2, (nz-output_z)/2 );
00948 //              FloatSize region_size( output_x, output_y, output_z);
00949 //              Region clip_region( origin, region_size );
00950 //              image->clip_inplace( clip_region );
00951 //      }
00952 
00953         // Should be an "if (verbose)" here or something
00954         //print_stats(quality_scores);
00955 
00956         image->update();
00957         
00958         if (params.has_key("savenorm") && strlen((const char *)params["savenorm"])>0) {
00959                 if (tmp_data->get_ysize()%2==0 && tmp_data->get_zsize()%2==0) tmp_data->process_inplace("xform.fourierorigin.tocenter");
00960                 tmp_data->write_image((const char *)params["savenorm"]);
00961         }
00962 
00963         delete tmp_data;
00964         tmp_data=0;
00965         //Since we give up the ownership of the pointer to-be-returned, it's caller's responsibility to delete the returned image.
00966         //So we wrap this function with return_value_policy< manage_new_object >() in libpyReconstructor2.cpp to hand over ownership to Python.
00967         EMData *ret=image;
00968         image=0;
00969         
00970         return ret;
00971 }
00972 
00973 int WienerFourierReconstructor::insert_slice(const EMData* const input_slice, const Transform & arg, const float weight)
00974 {
00975         // Are these exceptions really necessary? (d.woolford)
00976         if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
00977 
00978         Transform * rotation;
00979 /*      if ( input_slice->has_attr("xform.projection") ) {
00980                 rotation = (Transform*) (input_slice->get_attr("xform.projection")); // assignment operator
00981         } else {*/
00982         rotation = new Transform(arg); // assignment operator
00983 //      }
00984 
00985         if (!input_slice->has_attr("ctf_snr_total")) 
00986                 throw NotExistingObjectException("ctf_snr_total","No SNR information present in class-average. Must use the ctf.auto or ctfw.auto averager.");
00987 
00988         EMData *slice;
00989         if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
00990         else slice = preprocess_slice( input_slice, *rotation);
00991 
00992 
00993         // We must use only the rotational component of the transform, scaling, translation and mirroring
00994         // are not implemented in Fourier space, but are in preprocess_slice
00995         rotation->set_scale(1.0);
00996         rotation->set_mirror(false);
00997         rotation->set_trans(0,0,0);
00998 
00999         // Finally to the pixel wise slice insertion
01000         do_insert_slice_work(slice, *rotation, weight);
01001 
01002         delete rotation; rotation=0;
01003         delete slice;
01004 
01005 //      image->update();
01006         return 0;
01007 }
01008 
01009 void WienerFourierReconstructor::do_insert_slice_work(const EMData* const input_slice, const Transform & arg,const float inweight)
01010 {
01011 
01012         vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
01013 
01014         float inx=(float)(input_slice->get_xsize());            // x/y dimensions of the input image
01015         float iny=(float)(input_slice->get_ysize());
01016         
01017         int undo_wiener=(int)input_slice->get_attr_default("ctf_wiener_filtered",0);    // indicates whether we need to undo a wiener filter before insertion
01018 //      if (undo_wiener) throw UnexpectedBehaviorException("wiener_fourier does not yet accept already Wiener filtered class-averages. Suggest using ctf.auto averager for now.");
01019         
01020         vector<float> snr=input_slice->get_attr("ctf_snr_total");
01021         float sub=1.0;
01022         if (inweight<0) sub=-1.0;
01023         float weight;
01024         
01025         for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
01026                 Transform t3d = arg*(*it);
01027                 for (int y = -iny/2; y < iny/2; y++) {
01028                         for (int x = 0; x <=  inx/2; x++) {
01029 
01030                                 float rx = (float) x/(inx-2.0f);        // coords relative to Nyquist=.5
01031                                 float ry = (float) y/iny;
01032 
01033                                 // This deals with the SNR weight
01034                                 float rn = (float)hypot(rx,ry);
01035                                 if (rn>=.5) continue;           // no SNR in the corners, and we're going to mask them later anyway
01036                                 rn*=snr.size()*2.0f;
01037                                 int rni=(int)floor(rn);
01038                                 if ((unsigned int)rni>=snr.size()-1) weight=snr[snr.size()-1]*sub;
01039                                 else {
01040                                         rn-=rni;
01041                                         weight=Util::linear_interpolate(snr[rni],snr[rni+1],rn);
01042                                 }
01043 //                              if (weight>500.0) printf("%f %d %d %f %f %d %f\n",weight,x,y,rx,ry,rni);
01044                                 
01045                                 Vec3f coord(rx,ry,0);
01046                                 coord = coord*t3d; // transpose multiplication
01047                                 float xx = coord[0]; // transformed coordinates in terms of Nyquist
01048                                 float yy = coord[1];
01049                                 float zz = coord[2];
01050 
01051                                 // Map back to real pixel coordinates in output volume
01052                                 xx=xx*(nx-2);
01053                                 yy=yy*ny;
01054                                 zz=zz*nz;
01055 
01056 //                              printf("%f\n",weight);
01057                                 if (undo_wiener) inserter->insert_pixel(xx,yy,zz,(input_slice->get_complex_at(x,y))*((weight+1.0f)/weight),weight*sub);
01058                                 else inserter->insert_pixel(xx,yy,zz,input_slice->get_complex_at(x,y),weight*sub);
01059                         }
01060                 }
01061         }
01062 }
01063 
01064 int WienerFourierReconstructor::determine_slice_agreement(EMData*  input_slice, const Transform & arg, const float weight,bool sub)
01065 {
01066         // Are these exceptions really necessary? (d.woolford)
01067         if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
01068 
01069         Transform * rotation;
01070         rotation = new Transform(arg); // assignment operator
01071 
01072         EMData *slice;
01073         if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
01074         else slice = preprocess_slice( input_slice, *rotation);
01075 
01076 
01077         // We must use only the rotational component of the transform, scaling, translation and mirroring
01078         // are not implemented in Fourier space, but are in preprocess_slice
01079         rotation->set_scale(1.0);
01080         rotation->set_mirror(false);
01081         rotation->set_trans(0,0,0);
01082 
01083 //      tmp_data->write_image("dbug.hdf",0);
01084         
01085         // Remove the current slice first (not threadsafe, but otherwise performance would be awful)
01086         if (sub) do_insert_slice_work(slice, *rotation, -weight);
01087 
01088         // Compare
01089         do_compare_slice_work(slice, *rotation,weight);
01090 
01091         input_slice->set_attr("reconstruct_norm",slice->get_attr("reconstruct_norm"));
01092         input_slice->set_attr("reconstruct_absqual",slice->get_attr("reconstruct_absqual"));
01093 //      input_slice->set_attr("reconstruct_qual",slice->get_attr("reconstruct_qual"));
01094         input_slice->set_attr("reconstruct_weight",slice->get_attr("reconstruct_weight"));
01095 
01096         // Now put the slice back
01097         if (sub) do_insert_slice_work(slice, *rotation, weight);
01098 
01099 
01100         delete rotation;
01101         delete slice;
01102 
01103 //      image->update();
01104         return 0;
01105 
01106 }
01107 
01108 void WienerFourierReconstructor::do_compare_slice_work(EMData* input_slice, const Transform & arg,float weight)
01109 {
01110 
01111         float dt[3];    // This stores the complex and weight from the volume
01112         float dt2[2];   // This stores the local image complex
01113         float *dat = input_slice->get_data();
01114         vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
01115 
01116         float inx=(float)(input_slice->get_xsize());            // x/y dimensions of the input image
01117         float iny=(float)(input_slice->get_ysize());
01118 
01119         double dot=0;           // summed pixel*weight dot product
01120         double vweight=0;               // sum of weights
01121         double power=0;         // sum of inten*weight from volume
01122         double power2=0;                // sum of inten*weight from image
01123         for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
01124                 Transform t3d = arg*(*it);
01125                 for (int y = -iny/2; y < iny/2; y++) {
01126                         for (int x = 0; x <=  inx/2; x++) {
01127                                 if (x==0 && y==0) continue;             // We don't want to use the Fourier origin
01128 
01129                                 float rx = (float) x/(inx-2);   // coords relative to Nyquist=.5
01130                                 float ry = (float) y/iny;
01131 
01132 //                              if ((rx * rx + Util::square(ry - max_input_dim / 2)) > rl)
01133 //                                      continue;
01134 
01135                                 Vec3f coord(rx,ry,0);
01136                                 coord = coord*t3d; // transpose multiplication
01137                                 float xx = coord[0]; // transformed coordinates in terms of Nyquist
01138                                 float yy = coord[1];
01139                                 float zz = coord[2];
01140 
01141 
01142                                 if (fabs(xx)>0.5 || fabs(yy)>=0.5 || fabs(zz)>=0.5) continue;
01143 
01144                                 // Map back to actual pixel coordinates in output volume
01145                                 xx=xx*(nx-2);
01146                                 yy=yy*ny;
01147                                 zz=zz*nz;
01148 
01149 
01150                                 int idx = (int)(x * 2 + inx*(y<0?iny+y:y));
01151                                 dt2[0] = dat[idx];
01152                                 dt2[1] = dat[idx+1];
01153 
01154                                 // value returned indirectly in dt
01155                                 if (!pixel_at(xx,yy,zz,dt) || dt[2]<=0) continue;
01156 
01157 //                              printf("%f\t%f\t%f\t%f\t%f\n",dt[0],dt[1],dt[2],dt2[0],dt2[1]);
01158                                 dot+=(dt[0]*dt2[0]+dt[1]*dt2[1])*dt[2];
01159                                 vweight+=dt[2];
01160                                 power+=(dt[0]*dt[0]+dt[1]*dt[1])*dt[2];
01161                                 power2+=(dt2[0]*dt2[0]+dt2[1]*dt2[1])*dt[2];
01162                         }
01163                 }
01164         }
01165 
01166         dot/=sqrt(power*power2);                // normalize the dot product
01167 //      input_slice->set_attr("reconstruct_norm",(float)(power2<=0?1.0:sqrt(power/power2)/(inx*iny)));
01168         input_slice->set_attr("reconstruct_norm",(float)(power2<=0?1.0:sqrt(power/power2)));
01169         input_slice->set_attr("reconstruct_absqual",(float)dot);
01170         float rw=weight<=0?1.0f:1.0f/weight;
01171         input_slice->set_attr("reconstruct_qual",(float)(dot*rw/((rw-1.0)*dot+1.0)));   // here weight is a proxy for SNR
01172         input_slice->set_attr("reconstruct_weight",(float)vweight/(float)(subnx*subny*subnz));
01173 //      printf("** %g\t%g\t%g\t%g ##\n",dot,vweight,power,power2);
01174         //printf("** %f %f %f ##\n",(float)(power2<=0?1.0:sqrt(power/power2)/(inx*iny)),(float)dot,(float)(dot*weight/((weight-1.0)*dot+1.0)));
01175 }
01176 
01177 bool WienerFourierReconstructor::pixel_at(const float& xx, const float& yy, const float& zz, float *dt)
01178 {
01179         int x0 = (int) floor(xx);
01180         int y0 = (int) floor(yy);
01181         int z0 = (int) floor(zz);
01182         
01183         float *rdata=image->get_data();
01184         float *norm=tmp_data->get_data();
01185         float normsum=0,normsum2=0;
01186 
01187         dt[0]=dt[1]=dt[2]=0.0;
01188 
01189         if (nx==subnx) {                        // normal full reconstruction
01190                 if (x0<-nx2-1 || y0<-ny2-1 || z0<-nz2-1 || x0>nx2 || y0>ny2 || z0>nz2 ) return false;
01191 
01192                 // no error checking on add_complex_fast, so we need to be careful here
01193                 int x1=x0+1;
01194                 int y1=y0+1;
01195                 int z1=z0+1;
01196                 if (x0<-nx2) x0=-nx2;
01197                 if (x1>nx2) x1=nx2;
01198                 if (y0<-ny2) y0=-ny2;
01199                 if (y1>ny2) y1=ny2;
01200                 if (z0<-nz2) z0=-nz2;
01201                 if (z1>nz2) z1=nz2;
01202                 
01203                 size_t idx=0;
01204                 float r, gg;
01205                 for (int k = z0 ; k <= z1; k++) {
01206                         for (int j = y0 ; j <= y1; j++) {
01207                                 for (int i = x0; i <= x1; i ++) {
01208                                         r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
01209                                         idx=image->get_complex_index_fast(i,j,k);
01210                                         gg = Util::fast_exp(-r / EMConsts::I2G);
01211                                         
01212                                         dt[0]+=gg*rdata[idx];
01213                                         dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1];
01214                                         dt[2]+=norm[idx/2]*gg;
01215                                         normsum2+=gg;
01216                                         normsum+=gg*norm[idx/2];                                
01217                                 }
01218                         }
01219                 }
01220                 if (normsum==0) return false;
01221                 dt[0]/=normsum;
01222                 dt[1]/=normsum;
01223                 dt[2]/=normsum2;
01224 //              printf("%1.2f,%1.2f,%1.2f\t%1.3f\t%1.3f\t%1.3f\t%1.3f\t%1.3f\n",xx,yy,zz,dt[0],dt[1],dt[2],rdata[idx],rdata[idx+1]);
01225                 return true;
01226         } 
01227         else {                                  // for subvolumes, not optimized yet
01228                 size_t idx;
01229                 float r, gg;
01230                 for (int k = z0 ; k <= z0 + 1; k++) {
01231                         for (int j = y0 ; j <= y0 + 1; j++) {
01232                                 for (int i = x0; i <= x0 + 1; i ++) {
01233                                         r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
01234                                         idx=image->get_complex_index(i,j,k,subx0,suby0,subz0,nx,ny,nz);
01235                                         gg = Util::fast_exp(-r / EMConsts::I2G)*norm[idx/2];
01236                                         
01237                                         dt[0]+=gg*rdata[idx];
01238                                         dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1];
01239                                         dt[2]+=norm[idx/2];
01240                                         normsum+=gg;                            
01241                                 }
01242                         }
01243                 }
01244                 
01245                 if (normsum==0)  return false;
01246                 return true;
01247         }
01248 }
01249 
01250 
01251 EMData *WienerFourierReconstructor::finish(bool doift)
01252 {
01253 
01254         bool sqrtnorm=params.set_default("sqrtnorm",false);
01255         normalize_threed(sqrtnorm,true);                // true is the wiener filter
01256 
01257         if (doift) {
01258                 image->do_ift_inplace();
01259                 image->depad();
01260                 image->process_inplace("xform.phaseorigin.tocenter");
01261         }
01262 
01263         image->update();
01264         
01265         if (params.has_key("savenorm") && strlen((const char *)params["savenorm"])>0) {
01266                 if (tmp_data->get_ysize()%2==0 && tmp_data->get_zsize()%2==0) tmp_data->process_inplace("xform.fourierorigin.tocenter");
01267                 tmp_data->write_image((const char *)params["savenorm"]);
01268         }
01269 
01270         delete tmp_data;
01271         tmp_data=0;
01272         EMData *ret=image;
01273         image=0;
01274         
01275         return ret;
01276 }
01277 
01278 /*
01279 void BaldwinWoolfordReconstructor::setup()
01280 {
01281         //This is a bit of a hack - but for now it suffices
01282         params.set_default("mode","nearest_neighbor");
01283         WienerFourierReconstructor::setup();
01284         // Set up the Baldwin Kernel
01285         int P = (int)((1.0+0.25)*max_input_dim+1);
01286         float r = (float)(max_input_dim+1)/(float)P;
01287         dfreq = 0.2f;
01288         if (W != 0) delete [] W;
01289         int maskwidth = params.set_default("maskwidth",2);
01290         W = Util::getBaldwinGridWeights(maskwidth, (float)P, r,dfreq,0.5f,0.2f);
01291 }
01292 
01293 EMData* BaldwinWoolfordReconstructor::finish(bool doift)
01294 {
01295         tmp_data->write_image("density.mrc");
01296         image->process_inplace("xform.fourierorigin.tocorner");
01297         image->do_ift_inplace();
01298         image->depad();
01299         image->process_inplace("xform.phaseorigin.tocenter");
01300 
01301         if ( (bool) params.set_default("postmultiply", false) )
01302         {
01303                 cout << "POST MULTIPLYING" << endl;
01304         // now undo the Fourier convolution with real space division
01305                 float* d = image->get_data();
01306                 float N = (float) image->get_xsize()/2.0f;
01307                 N *= N;
01308                 size_t rnx = image->get_xsize();
01309                 size_t rny = image->get_ysize();
01310                 size_t rnxy = rnx*rny;
01311                 int cx = image->get_xsize()/2;
01312                 int cy = image->get_ysize()/2;
01313                 int cz = image->get_zsize()/2;
01314                 size_t idx;
01315                 for (int k = 0; k < image->get_zsize(); ++k ){
01316                         for (int j = 0; j < image->get_ysize(); ++j ) {
01317                                 for (int i =0; i < image->get_xsize(); ++i ) {
01318                                         float xd = (float)(i-cx); xd *= xd;
01319                                         float yd = (float)(j-cy); yd *= yd;
01320                                         float zd = (float)(k-cz); zd *= zd;
01321                                         float weight = exp((xd+yd+zd)/N);
01322                                         idx = k*rnxy + j*rnx + i;
01323                                         d[idx] *=  weight;
01324                                 }
01325                         }
01326                 }
01327         }
01328         image->update();
01329         return  image;
01330 }
01331 
01332 #include <iomanip>
01333 
01334 // int BaldwinWoolfordReconstructor::insert_slice_weights(const Transform& t3d)
01335 // {
01336 //      bool fftodd = image->is_fftodd();
01337 //      int rnx = nx-2*!fftodd;
01338 //
01339 //      float y_scale = 1.0, x_scale = 1.0;
01340 //
01341 //      if ( ny != rnx  )
01342 //      {
01343 //              if ( rnx > ny ) y_scale = (float) rnx / (float) ny;
01344 //              else x_scale = (float) ny / (float) rnx;
01345 //      }
01346 //
01347 //      int tnx = tmp_data->get_xsize();
01348 //      int tny = tmp_data->get_ysize();
01349 //      int tnz = tmp_data->get_zsize();
01350 //
01351 //      vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
01352 //      for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
01353 //              Transform n3d = t3d*(*it);
01354 //
01355 //              for (int y = 0; y < tny; y++) {
01356 //                      for (int x = 0; x < tnx; x++) {
01357 //
01358 //                              float rx = (float) x;
01359 //                              float ry = (float) y;
01360 //
01361 //                              if ( ny != rnx )
01362 //                              {
01363 //                                      if ( rnx > ny ) ry *= y_scale;
01364 //                                      else rx *= x_scale;
01365 //                              }
01366 // //                           float xx = rx * n3d[0][0] + (ry - tny/2) * n3d[1][0];
01367 // //                           float yy = rx * n3d[0][1] + (ry - tny/2) * n3d[1][1];
01368 // //                           float zz = rx * n3d[0][2] + (ry - tny/2) * n3d[1][2];
01369 //
01370 //                              Vec3f coord(rx,(ry - tny/2),0);
01371 //                              coord = coord*n3d; // transpose multiplication
01372 //                              float xx = coord[0];
01373 //                              float yy = coord[1];
01374 //                              float zz = coord[2];
01375 //
01376 //                              if (xx < 0 ){
01377 //                                      xx = -xx;
01378 //                                      yy = -yy;
01379 //                                      zz = -zz;
01380 //                              }
01381 //
01382 //                              yy += tny/2;
01383 //                              zz += tnz/2;
01384 //                              insert_density_at(xx,yy,zz);
01385 //                      }
01386 //              }
01387 //      }
01388 //
01389 //      return 0;
01390 // }
01391 
01392 void BaldwinWoolfordReconstructor::insert_density_at(const float& x, const float& y, const float& z)
01393 {
01394         int xl = Util::fast_floor(x);
01395         int yl = Util::fast_floor(y);
01396         int zl = Util::fast_floor(z);
01397 
01398         // w is the windowing width
01399         int w = params.set_default("maskwidth",2);
01400         float wsquared = (float) w*w;
01401         float dw = 1.0f/w;
01402         dw *= dw;
01403 
01404         // w minus one - this control the number of
01405         // pixels/voxels to the left of the main pixel
01406         // that will have density
01407         int wmox = w-1;
01408         int wmoy = w-1;
01409         int wmoz = w-1;
01410 
01411         // If any coordinate is incedental with a vertex, then
01412         // make sure there is symmetry in density accruing.
01413         // i.e. the window width must be equal in both directions
01414         if ( ((float) xl) == x ) wmox = w;
01415         if ( ((float) yl) == y ) wmoy = w;
01416         if ( ((float) zl) == z ) wmoz = w;
01417 
01418         float* d = tmp_data->get_data();
01419         int tnx = tmp_data->get_xsize();
01420         int tny = tmp_data->get_ysize();
01421         int tnz = tmp_data->get_zsize();
01422         size_t tnxy = tnx*tny;
01423 
01424         int mode = params.set_default("mode","nearest_neighbor");
01425 
01426         for(int k = zl-wmoz; k <= zl+w; ++k ) {
01427                 for(int j = yl-wmoy; j <= yl+w; ++j) {
01428                         for( int i = xl-wmox; i <= xl+w; ++i) {
01429                                 float fac = 1.0;
01430                                 int ic = i, jc = j, kc = k;
01431 
01432                                 // Fourier space is periodic, which is enforced
01433                                 // by the next 6 if statements. These if statements
01434                                 // assume that the Fourier DC components is at
01435                                 // (0,ny/2,nz/2).
01436                                 if ( i <= 0 ) {
01437 
01438                                         if ( x != 0 && i == 0 ) fac = 1.0;
01439                                         else if ( x == 0 && i < 0) continue;
01440 //                                      if (i < 0 ) ic = -i;
01441                                         if (i < 0 ) {
01442                                                 continue;
01443                                                 ic = -i;
01444                                                 jc = tny-jc;
01445                                                 kc = tnz-kc;
01446                                         }
01447                                 }
01448                                 if ( ic >= tnx ) ic = 2*tnx-ic-1;
01449                                 if ( jc < 0 ) jc = tny+jc;
01450                                 if ( jc >= tny ) jc = jc-tny;
01451                                 if ( kc < 0 ) kc = tnz+kc;
01452                                 if ( kc >= tnz ) kc = kc-tnz;
01453 //                              if ( ic >= tnx ) continue;
01454 //                              if ( jc < 0 ) continue;
01455 //                              if ( jc >= tny ) continue;
01456 //                              if ( kc < 0 ) continue;
01457 //                              if ( kc >= tnz ) continue;
01458                                 // This shouldn't happen
01459                                 // Debug remove later
01460                                 if ( ic < 0 ) { cout << "wo 1" << endl; }
01461                                 if ( ic >= tnx  ){ cout << "wo 2" << endl; }
01462                                 if ( jc < 0 ) { cout << "wo 3" << endl; }
01463                                 if ( jc >= tny ) { cout << "wo 4" << endl; }
01464                                 if ( kc < 0 ) { cout << "wo 5" << endl; }
01465                                 if ( kc >= tnz ) { cout << "wo 6" << endl; }
01466 
01467 
01468                                 float zd = (z-(float)k);
01469                                 float yd = (y-(float)j);
01470                                 float xd = (x-(float)i);
01471                                 zd *= zd; yd *= yd; xd *= xd;
01472                                 float distsquared = xd+yd+zd;
01473                                 // We enforce a spherical kernel
01474                                 if ( mode == 1 && distsquared > wsquared ) continue;
01475 
01476 //                              float f = fac*exp(-dw*(distsquared));
01477                                 float f = fac*exp(-2.467f*(distsquared));
01478                                 // Debug - this error should never occur.
01479                                 if ( (kc*tnxy+jc*tnx+ic) >= tnxy*tnz ) throw OutofRangeException(0,tnxy*tnz,kc*tnxy+jc*tnx+ic, "in density insertion" );
01480                                 d[kc*tnxy+jc*tnx+ic] += f;
01481                         }
01482                 }
01483         }
01484 }
01485 
01486 int BaldwinWoolfordReconstructor::insert_slice(const EMData* const input_slice, const Transform & t, const float weight)
01487 {
01488         Transform * rotation;
01489         if ( input_slice->has_attr("xform.projection") ) {
01490                 rotation = (Transform*) (input_slice->get_attr("xform.projection")); // assignment operator
01491         } else {
01492                 rotation = new Transform(t); // assignment operator
01493         }
01494         Transform tmp(*rotation);
01495         tmp.set_rotation(Dict("type","eman")); // resets the rotation to 0 implicitly
01496 
01497         Vec2f trans = tmp.get_trans_2d();
01498         float scale = tmp.get_scale();
01499         bool mirror = tmp.get_mirror();
01500         EMData* slice = 0;
01501         if (trans[0] != 0 || trans[1] != 0 || scale != 1.0 ) {
01502                 slice = input_slice->process("xform",Dict("transform",&tmp));
01503         } else if ( mirror == true ) {
01504                 slice = input_slice->process("xform.flip",Dict("axis","x"));
01505         }
01506         if ( slice == 0 ) {
01507                 slice = input_slice->process("xform.phaseorigin.tocorner");
01508         } else {
01509                 slice->process_inplace("xform.phaseorigin.tocorner");
01510         }
01511 
01512         slice->do_fft_inplace();
01513         slice->process_inplace("xform.fourierorigin.tocenter");
01514         float *dat = slice->get_data();
01515         float dt[2];
01516 
01517         bool fftodd = image->is_fftodd();
01518         int rnx = nx-2*!fftodd;
01519 
01520         float y_scale = 1.0, x_scale = 1.0;
01521 
01522         if ( ny != rnx  )
01523         {
01524                 if ( rnx > ny ) y_scale = (float) rnx / (float) ny;
01525                 else x_scale = (float) ny / (float) rnx;
01526         }
01527 
01528         int tnx = tmp_data->get_xsize();
01529         int tny = tmp_data->get_ysize();
01530         int tnz = tmp_data->get_zsize();
01531 
01532         vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
01533 //      float weight = params.set_default("weight",1.0f);
01534 
01535         rotation->set_scale(1.0); rotation->set_mirror(false); rotation->set_trans(0,0,0);
01536         for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
01537                 Transform t3d = (*rotation)*(*it);
01538 
01539                 for (int y = 0; y < tny; y++) {
01540                         for (int x = 0; x < tnx; x++) {
01541                                 float rx = (float) x;
01542                                 float ry = (float) y;
01543 
01544                                 if ( ny != rnx )
01545                                 {
01546                                         if ( rnx > ny ) ry *= y_scale;
01547                                         else rx *= x_scale;
01548                                 }
01549 
01550 //                              float xx = rx * n3d[0][0] + (ry - tny/2) * n3d[1][0];
01551 //                              float yy = rx * n3d[0][1] + (ry - tny/2) * n3d[1][1];
01552 //                              float zz = rx * n3d[0][2] + (ry - tny/2) * n3d[1][2];
01553 
01554                                 Vec3f coord(rx,(ry - tny/2),0);
01555                                 coord = coord*t3d; // transpose multiplication
01556                                 float xx = coord[0];
01557                                 float yy = coord[1];
01558                                 float zz = coord[2];
01559 
01560 
01561                                 float cc = 1;
01562                                 if (xx < 0 ){
01563                                         xx = -xx;
01564                                         yy = -yy;
01565                                         zz = -zz;
01566                                         cc = -1;
01567                                 }
01568 
01569                                 yy += tny/2;
01570                                 zz += tnz/2;
01571 
01572                                 int idx = x * 2 + y * (slice->get_xsize());
01573                                 dt[0] = dat[idx];
01574                                 dt[1] = cc * dat[idx+1];
01575 
01576                                 insert_pixel(xx,yy,zz,dt);
01577                         }
01578                 }
01579         }
01580 
01581         if(rotation) {delete rotation; rotation=0;}
01582         delete slice;
01583 
01584         return 0;
01585 }
01586 
01587 void BaldwinWoolfordReconstructor::insert_pixel(const float& x, const float& y, const float& z, const float dt[2])
01588 {
01589         int xl = Util::fast_floor(x);
01590         int yl = Util::fast_floor(y);
01591         int zl = Util::fast_floor(z);
01592 
01593         // w is the windowing width
01594         int w = params.set_default("maskwidth",2);
01595         float wsquared = (float) w*w;
01596         float dw = 1.0f/w;
01597         dw *= dw;
01598 
01599         int wmox = w-1;
01600         int wmoy = w-1;
01601         int wmoz = w-1;
01602 
01603         // If any coordinate is incedental with a vertex, then
01604         // make sure there is symmetry in density accruing.
01605         // i.e. the window width must be equal in both directions
01606         if ( ((float) xl) == x ) wmox = w;
01607         if ( ((float) yl) == y ) wmoy = w;
01608         if ( ((float) zl) == z ) wmoz = w;
01609 
01610         float* we = tmp_data->get_data();
01611         int tnx = tmp_data->get_xsize();
01612         int tny = tmp_data->get_ysize();
01613         int tnz = tmp_data->get_zsize();
01614         int tnxy = tnx*tny;
01615 
01616         int rnx = 2*tnx;
01617         int rnxy = 2*tnxy;
01618 
01619         int mode = params.set_default("mode","nearest_neighbor");
01620 
01621         float* d = image->get_data();
01622         for(int k = zl-wmoz; k <= zl+w; ++k ) {
01623                 for(int j = yl-wmoy; j <= yl+w; ++j) {
01624                         for( int i = xl-wmox; i <= xl+w; ++i) {
01625                                 float fac = 1.0;
01626                                 int ic = i, jc = j, kc = k;
01627 
01628                                 // Fourier space is periodic, which is enforced
01629                                 // by the next 6 if statements. These if statements
01630                                 // assume that the Fourier DC component is at
01631                                 // (0,ny/2,nz/2).
01632                                 float negfac=1.0;
01633                                 if ( i <= 0 ) {
01634                                         if ( x != 0 && i == 0 ) fac = 1.0;
01635                                         else if ( x == 0 && i < 0) continue;
01636                                         if (i < 0 ) {
01637                                                 continue;
01638                                                 ic = -i;
01639                                                 jc = tny-jc;
01640                                                 kc = tnz-kc;
01641                                                 negfac=-1.0;
01642                                         }
01643                                 }
01644                                 if ( ic >= tnx ) ic = 2*tnx-ic-1;
01645                                 if ( jc < 0 ) jc = tny+jc;
01646                                 if ( jc >= tny ) jc = jc-tny;
01647                                 if ( kc < 0 ) kc = tnz+kc;
01648                                 if ( kc >= tnz ) kc = kc-tnz;
01649 //                              if ( ic >= tnx ) continue;
01650 //                              if ( jc < 0 ) continue;
01651 //                              if ( jc >= tny ) continue;
01652 //                              if ( kc < 0 ) continue;
01653 //                              if ( kc >= tnz ) continue;
01654 
01655                                 float zd = (z-(float)k);
01656                                 float yd = (y-(float)j);
01657                                 float xd = (x-(float)i);
01658                                 zd *= zd; yd *= yd; xd *= xd;
01659                                 float distsquared = xd+yd+zd;
01660 //                              float f = fac*exp(-dw*(distsquared));
01661                                 float f = fac*exp(-2.467f*(distsquared));
01662                                 float weight = f/we[kc*tnxy+jc*tnx+ic];
01663                                 // debug - this error should never occur
01664                                 if ( (kc*rnxy+jc*rnx+2*ic+1) >= rnxy*tnz ) throw OutofRangeException(0,rnxy*tnz,kc*rnxy+jc*rnx+2*ic+1, "in pixel insertion" );
01665                                 size_t k = kc*rnxy+jc*rnx+2*ic;
01666 
01667                                 float factor, dist,residual;
01668                                 int sizeW,sizeWmid,idx;
01669                                 switch (mode) {
01670                                         case 0:
01671                                                 d[k] += weight*f*dt[0];
01672                                                 d[k+1] += negfac*weight*f*dt[1];
01673                                                 cout << "hello" << endl;
01674                                         break;
01675 
01676                                         case 1:
01677                                                 // We enforce a spherical kernel
01678                                                 if ( distsquared > wsquared ) continue;
01679 
01680                                                 sizeW = (int)(1+2*w/dfreq);
01681                                                 sizeWmid = sizeW/2;
01682 
01683                                                 dist = sqrtf(distsquared);
01684                                                 idx = (int)(sizeWmid + dist/dfreq);
01685                                                 if (idx >= sizeW) throw InvalidValueException(idx, "idx was greater than or equal to sizeW");
01686                                                 residual = dist/dfreq - (int)(dist/dfreq);
01687                                                 if ( fabs(residual) > 1) throw InvalidValueException(residual, "Residual was too big");
01688 
01689                                                 factor = (W[idx]*(1.0f-residual)+W[idx+1]*residual)*weight;
01690 
01691                                                 d[k] += dt[0]*factor;
01692                                                 d[k+1] += dt[1]*factor;
01693                                         break;
01694 
01695                                         default:
01696                                                 throw InvalidValueException(mode, "The mode was unsupported in BaldwinWoolfordReconstructor::insert_pixel");
01697                                         break;
01698                                 }
01699                         }
01700                 }
01701         }
01702 }
01703 
01704 // void BaldwinWoolfordReconstructor::insert_pixel(const float& x, const float& y, const float& z, const float dt[2])
01705 // {
01706 //      int xl = Util::fast_floor(x);
01707 //      int yl = Util::fast_floor(y);
01708 //      int zl = Util::fast_floor(z);
01709 //
01710 //      // w is the windowing width
01711 //      int w = params.set_default("maskwidth",2);
01712 //      float dw = 1.0/w;
01713 //      dw *= dw;
01714 // //   dw = 2;
01715 // //   cout << w << endl;
01716 //      //      int w = 3;
01717 //      // w minus one - this control the number of
01718 //      // pixels/voxels to the left of the main pixel
01719 //      // that will have density
01720 //      int wmox = w-1;
01721 //      int wmoy = w-1;
01722 //      int wmoz = w-1;
01723 //
01724 //      // If any coordinate is incedental with a vertex, then
01725 //      // make sure there is symmetry in density accruing.
01726 //      // i.e. the window width must be equal in both directions
01727 //      if ( ((float) xl) == x ) wmox = w;
01728 //      if ( ((float) yl) == y ) wmoy = w;
01729 //      if ( ((float) zl) == z ) wmoz = w;
01730 //
01731 //      float* d = tmp_data->get_data();
01732 //      int tnx = tmp_data->get_xsize();
01733 //      int tny = tmp_data->get_ysize();
01734 //      int tnz = tmp_data->get_zsize();
01735 //      int tnxy = tnx*tny;
01736 //
01737 //      float weight = 1.0;
01738 // //
01739 //      for(int k = zl-wmoz; k <= zl+w; ++k ) {
01740 //              for(int j = yl-wmoy; j <= yl+w; ++j) {
01741 //                      for( int i = xl-wmox; i <= xl+w; ++i) {
01742 //                              float fac = 1.0;
01743 //                              int ic = i, jc = j, kc = k;
01744 //
01745 //                              // Fourier space is periodic, which is enforced
01746 //                              // by the next 6 if statements. These if statements
01747 //                              // assume that the Fourier DC components is at
01748 //                              // (0,ny/2,nz/2).
01749 //                              if ( i <= 0 ) {
01750 //                                      if ( x != 0 && i == 0 ) fac = 1.0;
01751 //                                      else if ( x == 0 && i < 0) continue;
01752 // //                                   if (i < 0 ) ic = -i;
01753 //                                      if (i < 0 ) {
01754 //                                              ic = -i;
01755 //                                              jc = tny-jc;
01756 //                                              kc = tnz-kc;
01757 //                                      }
01758 //                              }
01759 //                              if ( ic >= tnx ) ic = 2*tnx-ic-1;
01760 //                              if ( jc < 0 ) jc = tny+jc;
01761 //                              if ( jc >= tny ) jc = jc-tny;
01762 //                              if ( kc < 0 ) kc = tnz+kc;
01763 //                              if ( kc >= tnz ) kc = kc-tnz;
01764 //                              // This shouldn't happen
01765 //                              // Debug remove later
01766 //                              if ( ic < 0 ) { cout << "wo 1" << endl; }
01767 //                              if ( ic >= tnx  ){ cout << "wo 2" << endl; }
01768 //                              if ( jc < 0 ) { cout << "wo 3" << endl; }
01769 //                              if ( jc >= tny ) { cout << "wo 4" << endl; }
01770 //                              if ( kc < 0 ) { cout << "wo 5" << endl; }
01771 //                              if ( kc >= tnz ) { cout << "wo 6" << endl; }
01772 //
01773 //
01774 //                              float zd = (z-(float)k);
01775 //                              float yd = (y-(float)j);
01776 //                              float xd = (x-(float)i);
01777 //                              zd *= zd; yd *= yd; xd *= xd;
01778 //                              // Debug - this error should never occur.
01779 //                              if ( (kc*tnxy+jc*tnx+ic) >= tnxy*tnz ) throw OutofRangeException(0,tnxy*tnz,kc*tnxy+jc*tnx+ic, "in weight determination insertion" );
01780 // //                           float f = fac*exp(-dw*(xd+yd+zd)*0.5);
01781 //                              float f = exp(-2.467*(xd+yd+zd));
01782 //                              weight += f*(d[kc*tnxy+jc*tnx+ic] - f);
01783 //                      }
01784 //              }
01785 //      }
01786 //      weight = 1.0/weight;
01787 //      int rnx = 2*tnx;
01788 //      int rnxy = 2*tnxy;
01789 //      d = image->get_data();
01790 //      for(int k = zl-wmoz; k <= zl+w; ++k ) {
01791 //              for(int j = yl-wmoy; j <= yl+w; ++j) {
01792 //                      for( int i = xl-wmox; i <= xl+w; ++i) {
01793 //                              float fac = 1.0;
01794 //                              int ic = i, jc = j, kc = k;
01795 //
01796 //                              // Fourier space is periodic, which is enforced
01797 //                              // by the next 6 if statements. These if statements
01798 //                              // assume that the Fourier DC components is at
01799 //                              // (0,ny/2,nz/2).
01800 //                              float negfac=1.0;
01801 //                              if ( i <= 0 ) {
01802 //                                      if ( x != 0 && i == 0 ) fac = 1.0;
01803 //                                      else if ( x == 0 && i < 0) continue;
01804 //                                      if (i < 0 ) {
01805 //                                              continue;
01806 //                                              ic = -i;
01807 //                                              jc = tny-jc;
01808 //                                              kc = tnz-kc;
01809 //                                              negfac=-1.0;
01810 //                                      }
01811 //                              }
01812 //                              if ( ic >= tnx ) ic = 2*tnx-ic-1;
01813 //                              if ( jc < 0 ) jc = tny+jc;
01814 //                              if ( jc >= tny ) jc = jc-tny;
01815 //                              if ( kc < 0 ) kc = tnz+kc;
01816 //                              if ( kc >= tnz ) kc = kc-tnz;
01817 //                              // This shouldn't happen
01818 //                              // Debug remove later
01819 //
01820 //
01821 //                              float zd = (z-(float)k);
01822 //                              float yd = (y-(float)j);
01823 //                              float xd = (x-(float)i);
01824 //                              zd *= zd; yd *= yd; xd *= xd;
01825 // //                           float f = fac*exp(-dw*(xd+yd+zd));
01826 //                              float f = exp(-4.934*(xd+yd+zd));
01827 //                              // Debug - this error should never occur.
01828 //                              if ( (kc*rnxy+jc*rnx+2*ic+1) >= rnxy*tnz ) throw OutofRangeException(0,rnxy*tnz,kc*rnxy+jc*rnx+2*ic+1, "in pixel insertion" );
01829 //
01830 //                              d[kc*rnxy+jc*rnx+2*ic] += weight*f*dt[0];
01831 //                              d[kc*rnxy+jc*rnx+2*ic+1] += negfac*weight*f*dt[1];
01832 //                      }
01833 //              }
01834 //      }
01835 // }
01836 */
01837 
01838 
01839 void BackProjectionReconstructor::setup()
01840 {
01841         int size = params["size"];
01842         image = new EMData();
01843         nx = size;
01844         ny = size;
01845         if ( (int) params["zsample"] != 0 ) nz = params["zsample"];
01846         else nz = size;
01847         image->set_size(nx, ny, nz);
01848 }
01849 
01850 EMData* BackProjectionReconstructor::preprocess_slice(const EMData* const slice, const Transform& t)
01851 {
01852 
01853         EMData* return_slice = slice->process("normalize.edgemean");
01854         return_slice->process_inplace("filter.linearfourier");
01855 
01856         Transform tmp(t);
01857         tmp.set_rotation(Dict("type","eman")); // resets the rotation to 0 implicitly
01858         Vec2f trans = tmp.get_trans_2d();
01859         float scale = tmp.get_scale();
01860         bool mirror = tmp.get_mirror();
01861         if (trans[0] != 0 || trans[1] != 0 || scale != 1.0 ) {
01862                 return_slice->transform(tmp);
01863         } else if ( mirror == true ) {
01864                 return_slice = slice->process("xform.flip",Dict("axis","x"));
01865         }
01866 
01867         return return_slice;
01868 }
01869 
01870 int BackProjectionReconstructor::insert_slice(const EMData* const input, const Transform &t, const float)
01871 {
01872         if (!input) {
01873                 LOGERR("try to insert NULL slice");
01874                 return 1;
01875         }
01876 
01877         if (input->get_xsize() != input->get_ysize() || input->get_xsize() != nx) {
01878                 LOGERR("tried to insert image that was not correction dimensions");
01879                 return 1;
01880         }
01881 
01882         Transform * transform;
01883         if ( input->has_attr("xform.projection") ) {
01884                 transform = (Transform*) (input->get_attr("xform.projection")); // assignment operator
01885         } else {
01886                 transform = new Transform(t); // assignment operator
01887         }
01888         EMData* slice = preprocess_slice(input, t);
01889 
01890         float weight = params["weight"];
01891         slice->mult(weight);
01892 
01893         EMData *tmp = new EMData();
01894         tmp->set_size(nx, ny, nz);
01895 
01896         float *slice_data = slice->get_data();
01897         float *tmp_data = tmp->get_data();
01898 
01899         size_t nxy = nx * ny;
01900         size_t nxy_size = nxy * sizeof(float);;
01901         for (int i = 0; i < nz; ++i) {
01902                 memcpy(&tmp_data[nxy * i], slice_data, nxy_size);
01903         }
01904 
01905         transform->set_scale(1.0);
01906         transform->set_mirror(false);
01907         transform->set_trans(0,0,0);
01908         transform->invert();
01909 
01910         tmp->transform(*transform);
01911         image->add(*tmp);
01912 
01913         if(transform) {delete transform; transform=0;}
01914         delete tmp;
01915         delete slice;
01916 
01917         return 0;
01918 }
01919 
01920 EMData *BackProjectionReconstructor::finish(bool)
01921 {
01922 
01923         Symmetry3D* sym = Factory<Symmetry3D>::get((string)params["sym"]);
01924         vector<Transform> syms = sym->get_syms();
01925 
01926         for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
01927 
01928                 EMData tmpcopy(*image);
01929                 tmpcopy.transform(*it);
01930                 image->add(tmpcopy);
01931         }
01932 
01933         image->mult(1.0f/(float)sym->get_nsym());
01934         delete sym;
01935 
01936         EMData *ret = image;
01937         image = 0 ;
01938         return ret;
01939 }
01940 
01941 EMData* EMAN::padfft_slice( const EMData* const slice, const Transform& t, int npad )
01942 {
01943         int nx = slice->get_xsize();
01944         int ny = slice->get_ysize();
01945         int ndim = (ny==1) ? 1 : 2;
01946 
01947         if( ndim==2 && nx!=ny )
01948         {
01949                 // FIXME: What kind of exception should we throw here?
01950                 throw std::runtime_error("Tried to padfft a 2D slice which is not square.");
01951         }
01952 
01953         // process 2D slice or 1D line -- subtract the average outside of the circle, zero-pad, fft extend, and fft
01954         EMData* temp = slice->average_circ_sub();
01955 
01956         Assert( temp != NULL );
01957         EMData* zeropadded = temp->norm_pad( false, npad );
01958         Assert( zeropadded != NULL );
01959         checked_delete( temp );
01960 
01961         zeropadded->do_fft_inplace();
01962         EMData* padfftslice = zeropadded;
01963 
01964         // shift the projection
01965         Vec2f trans = t.get_trans_2d();
01966         float sx = -trans[0];
01967         float sy = -trans[1];
01968         if(sx != 0.0f || sy != 0.0)
01969                 padfftslice->process_inplace("filter.shift", Dict("x_shift", sx, "y_shift", sy, "z_shift", 0.0f));
01970 
01971         int remove = slice->get_attr_default("remove", 0);
01972         padfftslice->set_attr( "remove", remove );
01973 
01974 
01975 
01976         padfftslice->center_origin_fft();
01977         return padfftslice;
01978 }
01979 
01980 nn4Reconstructor::nn4Reconstructor()
01981 {
01982         m_volume = NULL;
01983         m_wptr   = NULL;
01984 }
01985 
01986 nn4Reconstructor::nn4Reconstructor( const string& symmetry, int size, int npad )
01987 {
01988         m_volume = NULL;
01989         m_wptr   = NULL;
01990         setup( symmetry, size, npad );
01991         load_default_settings();
01992         print_params();
01993 }
01994 
01995 nn4Reconstructor::~nn4Reconstructor()
01996 {
01997         //if( m_delete_volume ) checked_delete(m_volume); 
01998 
01999         //if( m_delete_weight ) checked_delete( m_wptr );
02000 
02001         //checked_delete( m_result );
02002 }
02003 
02004 enum weighting_method { NONE, ESTIMATE, VORONOI };
02005 
02006 float max2d( int kc, const vector<float>& pow_a )
02007 {
02008         float max = 0.0;
02009         for( int i=-kc; i <= kc; ++i ) {
02010                 for( int j=-kc; j <= kc; ++j ) {
02011                         if( i==0 && j==0 ) continue;
02012                         {
02013                                 int c = 2*kc+1 - std::abs(i) - std::abs(j);
02014                                 max = max + pow_a[c];
02015                         }
02016                 }
02017         }
02018         return max;
02019 }
02020 
02021 float max3d( int kc, const vector<float>& pow_a )
02022 {
02023         float max = 0.0;
02024         for( int i=-kc; i <= kc; ++i ) {
02025                 for( int j=-kc; j <= kc; ++j ) {
02026                         for( int k=-kc; k <= kc; ++k ) {
02027                                 if( i==0 && j==0 && k==0 ) continue;
02028                                 // if( i!=0 )
02029                                 {
02030                                         int c = 3*kc+1 - std::abs(i) - std::abs(j) - std::abs(k);
02031                                         max = max + pow_a[c];
02032                                         // max = max + c * c;
02033                                         // max = max + c;
02034                                 }
02035                         }
02036                 }
02037         }
02038         return max;
02039 }
02040 
02041 
02042 void nn4Reconstructor::setup()
02043 {
02044         int size = params["size"];
02045         int npad = params["npad"];
02046 
02047 
02048         string symmetry;
02049         if( params.has_key("symmetry") )  symmetry = params["symmetry"].to_str();
02050         else                               symmetry = "c1";
02051 
02052         if( params.has_key("ndim") )  m_ndim = params["ndim"];
02053         else                          m_ndim = 3;
02054 
02055         if( params.has_key( "snr" ) )  m_osnr = 1.0f/float( params["snr"] );
02056         else                           m_osnr = 0.0;
02057 
02058         setup( symmetry, size, npad );
02059 }
02060 
02061 void nn4Reconstructor::setup( const string& symmetry, int size, int npad )
02062 {
02063         m_weighting = ESTIMATE;
02064         m_wghta = 0.2f;
02065 
02066         m_symmetry = symmetry;
02067         m_npad = npad;
02068         m_nsym = Transform::get_nsym(m_symmetry);
02069 
02070         m_vnx = size;
02071         m_vny = size;
02072         m_vnz = (m_ndim==3) ? size : 1;
02073 
02074         m_vnxp = size*npad;
02075         m_vnyp = size*npad;
02076         m_vnzp = (m_ndim==3) ? size*npad : 1;
02077 
02078         m_vnxc = m_vnxp/2;
02079         m_vnyc = m_vnyp/2;
02080         m_vnzc = (m_ndim==3) ? m_vnzp/2 : 1;
02081 
02082         buildFFTVolume();
02083         buildNormVolume();
02084         
02085 }
02086 
02087 
02088 void nn4Reconstructor::buildFFTVolume() {
02089         int offset = 2 - m_vnxp%2;
02090 
02091         m_volume = params["fftvol"];
02092 
02093         if( m_volume->get_xsize() != m_vnxp+offset && m_volume->get_ysize() != m_vnyp && m_volume->get_zsize() != m_vnzp ) {
02094                 m_volume->set_size(m_vnxp+offset,m_vnyp,m_vnzp);
02095                 m_volume->to_zero();
02096         }
02097         // ----------------------------------------------------------------
02098         // Added by Zhengfan Yang on 03/15/07
02099         // Original author: please check whether my revision is correct and
02100         // other Reconstructor need similiar revision.
02101         if ( m_vnxp % 2 == 0 )  m_volume->set_fftodd(0);
02102         else                    m_volume->set_fftodd(1);
02103         // ----------------------------------------------------------------
02104 
02105         m_volume->set_nxc(m_vnxp/2);
02106         m_volume->set_complex(true);
02107         m_volume->set_ri(true);
02108         m_volume->set_fftpad(true);
02109         m_volume->set_attr("npad", m_npad);
02110         m_volume->set_array_offsets(0,1,1);
02111 }
02112 
02113 void nn4Reconstructor::buildNormVolume() {
02114 
02115         m_wptr = params["weight"];
02116 
02117         if( m_wptr->get_xsize() != m_vnxc+1 &&
02118                 m_wptr->get_ysize() != m_vnyp &&
02119                 m_wptr->get_zsize() != m_vnzp ) {
02120                 m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
02121                 m_wptr->to_zero();
02122         }
02123         m_wptr->set_array_offsets(0,1,1);
02124 }
02125 
02126 void printImage( const EMData* line )
02127 {
02128         Assert( line->get_zsize()==1 );
02129 
02130 
02131         int nx = line->get_xsize();
02132         int ny = line->get_ysize();
02133         for( int j=0; j < ny; ++j ) {
02134                 for( int i=0; i < nx; ++i )  printf( "%10.3f ", line->get_value_at(i,j) );
02135                 printf( "\n" );
02136         }
02137 }
02138 
02139 
02140 
02141 int nn4Reconstructor::insert_slice(const EMData* const slice, const Transform& t, const float) {
02142         // sanity checks
02143         if (!slice) {
02144                 LOGERR("try to insert NULL slice");
02145                 return 1;
02146         }
02147 
02148         int padffted= slice->get_attr_default( "padffted", 0 );
02149         if( m_ndim==3 ) {
02150                 if ( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_vnx)  ) {
02151                         // FIXME: Why doesn't this throw an exception?
02152                         LOGERR("Tried to insert a slice that is the wrong size.");
02153                         return 1;
02154                 }
02155         } else {
02156                 Assert( m_ndim==2 );
02157                 if( slice->get_ysize() !=1 ) {
02158                         LOGERR( "for 2D reconstruction, a line is excepted" );
02159                         return 1;
02160                 }
02161         }
02162 
02163         EMData* padfft = NULL;
02164 
02165         if( padffted != 0 ) padfft = new EMData(*slice);
02166         else                padfft = padfft_slice( slice, t,  m_npad );
02167 
02168         int mult= slice->get_attr_default( "mult", 1 );
02169         Assert( mult > 0 );
02170 
02171         if( m_ndim==3 ) {
02172                 insert_padfft_slice( padfft, t, mult );
02173         } else {
02174                 float alpha = padfft->get_attr( "alpha" );
02175                 alpha = alpha/180.0f*M_PI;
02176                 for(int i=0; i < m_vnxc+1; ++i ) {
02177                         float xnew = i*cos(alpha);
02178                         float ynew = -i*sin(alpha);
02179                         float btqr = padfft->get_value_at( 2*i, 0, 0 );
02180                         float btqi = padfft->get_value_at( 2*i+1, 0, 0 );
02181                         if( xnew < 0.0 ) {
02182                                 xnew *= -1;
02183                                 ynew *= -1;
02184                                 btqi *= -1;
02185                         }
02186 
02187                         int ixn = int(xnew+0.5+m_vnxp) - m_vnxp;
02188                         int iyn = int(ynew+0.5+m_vnyp) - m_vnyp;
02189 
02190                         if(iyn < 0 ) iyn += m_vnyp;
02191 
02192                         (*m_volume)( 2*ixn, iyn+1, 1 ) += btqr *float(mult);
02193                         (*m_volume)( 2*ixn+1, iyn+1, 1 ) += btqi * float(mult);
02194                         (*m_wptr)(ixn,iyn+1, 1) += float(mult);
02195                 }
02196 
02197         }
02198         checked_delete( padfft );
02199         return 0;
02200 }
02201 
02202 int nn4Reconstructor::insert_padfft_slice( EMData* padfft, const Transform& t, int mult )
02203 {
02204         Assert( padfft != NULL );
02205 
02206         vector<Transform> tsym = t.get_sym_proj(m_symmetry);
02207         for (unsigned int isym=0; isym < tsym.size(); isym++) {
02208                 m_volume->nn( m_wptr, padfft, tsym[isym], mult);
02209         }
02210         
02211         /*for (int isym=0; isym < m_nsym; isym++) {
02212                 Transform tsym = t.get_sym(m_symmetry, isym);
02213                 m_volume->nn( m_wptr, padfft, tsym, mult);
02214         }*/
02215         return 0;
02216 }
02217 
02218 
02219 #define  tw(i,j,k)      tw[ i-1 + (j-1+(k-1)*iy)*ix ]
02220 
02221 void circumf( EMData* win , int npad)
02222 {
02223         float *tw = win->get_data();
02224         //  correct for the fall-off
02225         //  mask and subtract circumference average
02226         int ix = win->get_xsize();
02227         int iy = win->get_ysize();
02228         int iz = win->get_zsize();
02229         int L2 = (ix/2)*(ix/2);
02230         int L2P = (ix/2-1)*(ix/2-1);
02231 
02232         int IP = ix/2+1;
02233         int JP = iy/2+1;
02234         int KP = iz/2+1;
02235 
02236         //  sinc functions tabulated for fall-off
02237         float* sincx = new float[IP+1];
02238         float* sincy = new float[JP+1];
02239         float* sincz = new float[KP+1];
02240 
02241         sincx[0] = 1.0f;
02242         sincy[0] = 1.0f;
02243         sincz[0] = 1.0f;
02244 
02245         float cdf = M_PI/float(npad*2*ix);
02246         for (int i = 1; i <= IP; ++i)  sincx[i] = sin(i*cdf)/(i*cdf);
02247         cdf = M_PI/float(npad*2*iy);
02248         for (int i = 1; i <= JP; ++i)  sincy[i] = sin(i*cdf)/(i*cdf);
02249         cdf = M_PI/float(npad*2*iz);
02250         for (int i = 1; i <= KP; ++i)  sincz[i] = sin(i*cdf)/(i*cdf);
02251         for (int k = 1; k <= iz; ++k) {
02252                 int kkp = abs(k-KP);
02253                 for (int j = 1; j <= iy; ++j) {
02254                         cdf = sincy[abs(j- JP)]*sincz[kkp];
02255                         for (int i = 1; i <= ix; ++i)  tw(i,j,k) /= (sincx[abs(i-IP)]*cdf);
02256                 }
02257         }
02258 
02259         delete[] sincx;
02260         delete[] sincy;
02261         delete[] sincz;
02262 
02263         float  TNR = 0.0f;
02264         size_t m = 0;
02265         for (int k = 1; k <= iz; ++k) {
02266                 for (int j = 1; j <= iy; ++j) {
02267                         for (int i = 1; i <= ix; ++i) {
02268                                 size_t LR = (k-KP)*(k-KP)+(j-JP)*(j-JP)+(i-IP)*(i-IP);
02269                                 if (LR >= (size_t)L2P && LR<=(size_t)L2) {
02270                                         TNR += tw(i,j,k);
02271                                         ++m;
02272                                 }
02273                         }
02274                 }
02275         }
02276 
02277         TNR /=float(m);
02278         
02279         
02280         for (int k = 1; k <= iz; ++k) {
02281                 for (int j = 1; j <= iy; ++j) {
02282                         for (int i = 1; i <= ix; ++i) {
02283                                 size_t LR = (k-KP)*(k-KP)+(j-JP)*(j-JP)+(i-IP)*(i-IP);
02284                                 if (LR<=(size_t)L2) tw(i,j,k) -= TNR;
02285                                 else                tw(i,j,k) = 0.0f;
02286 
02287                         }
02288                 }
02289         }
02290 
02291 }
02292 
02293 
02294 EMData* nn4Reconstructor::finish(bool) {
02295 
02296         if( m_ndim==3 ) {
02297                 m_volume->symplane0(m_wptr);
02298         } else {
02299                 for( int i=1; i <= m_vnyp; ++i ) {
02300 
02301                         if( (*m_wptr)(0, i, 1)==0.0 ) {
02302                                 int j = m_vnyp + 1 - i;
02303                                 (*m_wptr)(0, i, 1) = (*m_wptr)(0, j, 1);
02304                                 (*m_volume)(0, i, 1) = (*m_volume)(0, j, 1);
02305                                 (*m_volume)(1, i, 1) = (*m_volume)(1, j, 1);
02306                         }
02307                 }
02308         }
02309 
02310 
02311         int box = 7;
02312         int kc = (box-1)/2;
02313         vector< float > pow_a( m_ndim*kc+1, 1.0 );
02314         for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
02315         pow_a.back()=0.0;
02316 
02317         float alpha = 0.0;
02318         if( m_ndim==3) {
02319                 int vol = box*box*box;
02320                 float max = max3d( kc, pow_a );
02321                 alpha = ( 1.0f - 1.0f/(float)vol ) / max;
02322         } else {
02323                 int ara = box*box;
02324                 float max = max2d( kc, pow_a );
02325                 alpha = ( 1.0f - 1.0f/(float)ara ) / max;
02326         }
02327 
02328         int ix,iy,iz;
02329         for (iz = 1; iz <= m_vnzp; iz++) {
02330                 for (iy = 1; iy <= m_vnyp; iy++) {
02331                         for (ix = 0; ix <= m_vnxc; ix++) {
02332                                 if ( (*m_wptr)(ix,iy,iz) > 0) {//(*v) should be treated as complex!!
02333                                         float tmp;
02334                                         tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+m_osnr);
02335 
02336                                         if( m_weighting == ESTIMATE ) {
02337                                                 int cx = ix;
02338                                                 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
02339                                                 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
02340                                                 float sum = 0.0;
02341                                                 for( int ii = -kc; ii <= kc; ++ii ) {
02342                                                         int nbrcx = cx + ii;
02343                                                         if( nbrcx >= m_vnxc ) continue;
02344                                                         for( int jj= -kc; jj <= kc; ++jj ) {
02345                                                                 int nbrcy = cy + jj;
02346                                                                 if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
02347 
02348                                                                 int kcz = (m_ndim==3) ? kc : 0;
02349                                                                 for( int kk = -kcz; kk <= kcz; ++kk ) {
02350                                                                         int nbrcz = cz + kk;
02351                                                                         if( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
02352                                                                         if( nbrcx < 0 ) {
02353                                                                                 nbrcx = -nbrcx;
02354                                                                                 nbrcy = -nbrcy;
02355                                                                                 nbrcz = -nbrcz;
02356                                                                         }
02357                                                                         int nbrix = nbrcx;
02358                                                                         int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
02359                                                                         int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
02360                                                                         if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0 ) {
02361                                                                                 int c = m_ndim*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
02362                                                                                 sum = sum + pow_a[c];
02363                                                                         }
02364                                                                 }
02365                                                         }
02366                                                 }
02367                                                 float wght = 1.0f / ( 1.0f - alpha * sum );
02368                                                 tmp = tmp * wght;
02369                                         }
02370                                         (*m_volume)(2*ix,iy,iz)   *= tmp;
02371                                         (*m_volume)(2*ix+1,iy,iz) *= tmp;
02372                                 }
02373                         }
02374                 }
02375         }
02376 
02377         //if(m_ndim==2) printImage( m_volume );
02378 
02379         // back fft
02380         m_volume->do_ift_inplace();
02381 
02382         // EMData* win = m_volume->window_center(m_vnx);
02383         int npad = m_volume->get_attr("npad");
02384         m_volume->depad();
02385         circumf( m_volume, npad );
02386         m_volume->set_array_offsets( 0, 0, 0 );
02387 
02388         return 0;
02389 }
02390 #undef  tw
02391 
02392 
02393 nn4_rectReconstructor::nn4_rectReconstructor()
02394 {
02395         m_volume = NULL;
02396         m_wptr   = NULL;
02397 }
02398 
02399 nn4_rectReconstructor::nn4_rectReconstructor( const string& symmetry, int size, int npad )
02400 {
02401         m_volume = NULL;
02402         m_wptr   = NULL;
02403         setup( symmetry, size, npad );
02404         load_default_settings();
02405         print_params();
02406 }
02407 
02408 nn4_rectReconstructor::~nn4_rectReconstructor()
02409 {
02410         //if( m_delete_volume ) checked_delete(m_volume);
02411 
02412         //if( m_delete_weight ) checked_delete( m_wptr );
02413 
02414         //checked_delete( m_result );
02415 }
02416 
02417 
02418 void nn4_rectReconstructor::setup()
02419 {
02420         m_sizeofprojection = params["sizeprojection"];
02421         int npad = params["npad"];
02422         m_count=0;
02423 
02424         string symmetry;
02425         if( params.has_key("symmetry") )  symmetry = params["symmetry"].to_str();
02426         else                               symmetry = "c1";
02427 
02428         if( params.has_key("ndim") )  m_ndim = params["ndim"];             
02429         else                            m_ndim = 3;
02430     
02431         if( params.has_key( "snr" ) )  m_osnr = 1.0f/float( params["snr"] );
02432         else                           m_osnr = 0.0;
02433 
02434         setup( symmetry, m_sizeofprojection, npad );
02435 }
02436 
02437 void nn4_rectReconstructor::setup( const string& symmetry, int sizeprojection, int npad )
02438 {
02439         m_weighting = ESTIMATE;
02440         m_wghta = 0.2f;
02441         m_symmetry = symmetry;
02442         m_npad = npad;
02443         m_nsym = Transform::get_nsym(m_symmetry);
02444 
02445         if( params.has_key("sizex") )  m_vnx = params["sizex"];
02446         else if(params.has_key("xratio")) 
02447                 {
02448                 float temp=params["xratio"];
02449                 m_vnx=int(float(sizeprojection)*temp);
02450                 }
02451         else                           m_vnx=sizeprojection;
02452 
02453         if( params.has_key("sizey") )  m_vny = params["sizey"];
02454         else if (params.has_key("yratio"))  
02455                {
02456                 float temp=params["yratio"];
02457                  m_vny=int(float(sizeprojection)*temp);
02458                 }
02459         else m_vny=sizeprojection;
02460 
02461         if( params.has_key("sizez") ) 
02462                 m_vnz = params["sizez"];
02463         else                          
02464                 m_vnz = (m_ndim==3) ? sizeprojection : 1;
02465         
02466         m_xratio=float(m_vnx)/float(sizeprojection);    
02467         m_yratio=float(m_vny)/float(sizeprojection);
02468 
02469         m_vnxp = m_vnx*npad;
02470         m_vnyp = m_vny*npad;
02471         m_vnzp = (m_ndim==3) ? m_vnz*npad : 1;
02472 
02473         m_vnxc = m_vnxp/2;
02474         m_vnyc = m_vnyp/2;
02475         m_vnzc = (m_ndim==3) ? m_vnzp/2 : 1;
02476 
02477         buildFFTVolume();
02478         buildNormVolume();
02479 }
02480 
02481 
02482 void nn4_rectReconstructor::buildFFTVolume() {
02483         int offset = 2 - m_vnxp%2;
02484 
02485         m_volume = params["fftvol"];
02486 
02487         if( m_volume->get_xsize() != m_vnxp+offset && m_volume->get_ysize() != m_vnyp && m_volume->get_zsize() != m_vnzp ) {
02488                 m_volume->set_size(m_vnxp+offset,m_vnyp,m_vnzp);
02489                 m_volume->to_zero();
02490         }
02491         // ----------------------------------------------------------------
02492         // Added by Zhengfan Yang on 03/15/07
02493         // Original author: please check whether my revision is correct and
02494         // other Reconstructor need similiar revision.
02495         if ( m_vnxp % 2 == 0 )  m_volume->set_fftodd(0);
02496         else                    m_volume->set_fftodd(1);
02497         // ----------------------------------------------------------------
02498 
02499         m_volume->set_nxc(m_vnxp/2);
02500         m_volume->set_complex(true);
02501         m_volume->set_ri(true);
02502         m_volume->set_fftpad(true);
02503         m_volume->set_attr("npad", m_npad);
02504         m_volume->set_array_offsets(0,1,1);
02505 }
02506 
02507 void nn4_rectReconstructor::buildNormVolume() {
02508 
02509         m_wptr = params["weight"];
02510 
02511         if( m_wptr->get_xsize() != m_vnxc+1 &&
02512                 m_wptr->get_ysize() != m_vnyp &&
02513                 m_wptr->get_zsize() != m_vnzp ) {
02514                 m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
02515                 m_wptr->to_zero();
02516         }
02517 
02518         m_wptr->set_array_offsets(0,1,1);
02519 }
02520 
02521 int nn4_rectReconstructor::insert_slice(const EMData* const slice, const Transform& t, const float) {
02522         // sanity checks
02523 
02524 
02525         //m_count=m_count+1;
02526         if (!slice) {
02527                 LOGERR("try to insert NULL slice");
02528                 return 1;
02529         }
02530 
02531         int padffted= slice->get_attr_default( "padffted", 0 );
02532         if( m_ndim==3 ) {
02533                 if ( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_sizeofprojection)  ) {
02534                         
02535                         // FIXME: Why doesn't this throw an exception?
02536                         LOGERR("Tried to insert a slice that is the wrong size.");
02537                         return 1;
02538                 }
02539         } 
02540        if (m_ndim==2) {
02541                 Assert( m_ndim==2 );
02542                 if( slice->get_ysize() !=1 ) {
02543                         LOGERR( "for 2D reconstruction, a line is excepted" );
02544                         return 1;
02545                 }
02546         }
02547 
02548         EMData* padfft = NULL;
02549 
02550         if( padffted != 0 ) padfft = new EMData(*slice);
02551         else                padfft = padfft_slice( slice, t,  m_npad );
02552 
02553         int mult= slice->get_attr_default( "mult", 1 );
02554         Assert( mult > 0 );
02555 
02556         if( m_ndim==3 ) {
02557                 insert_padfft_slice( padfft, t, mult );         
02558                 
02559         } else {
02560                 float ellipse_length,ellipse_step,cos_alpha,sin_alpha;
02561                 int ellipse_length_int;
02562                 float alpha = padfft->get_attr( "alpha" );
02563                 alpha = alpha/180.0f*M_PI;
02564                 int loop_range;
02565                 float temp1,temp2;
02566                 temp1=m_xratio*cos(alpha)*float(m_sizeofprojection*m_npad)/2;
02567                 temp2=m_yratio*sin(alpha)*float(m_sizeofprojection*m_npad)/2;
02568                 ellipse_length=sqrt(temp1*temp1+temp2*temp2);
02569                 ellipse_length_int=int(ellipse_length);
02570                 ellipse_step=0.5f*(m_sizeofprojection*m_npad)/float(ellipse_length_int);
02571                 loop_range=ellipse_length_int;
02572                 cos_alpha=temp1/ellipse_length;
02573                 sin_alpha=temp2/ellipse_length;
02574                 if(m_count%100==0)
02575                 {
02576                         std::cout<<"#############################################################"<<std::endl;
02577                         std::cout<<"line insert start=="<<m_count<<std::endl;
02578                         std::cout<<"ellipse lenth=="<<ellipse_length_int<<"ellips step=="<<ellipse_step<<std::endl;
02579                         std::cout<<"loop_range"<<loop_range<<std::endl;
02580                         std::cout<<"x and y ratio=="<<m_xratio<<"  "<<m_yratio<<std::endl;
02581                         std::cout<<"cos sin of alpha=="<<cos(alpha)<<"   "<<sin(alpha)<<std::endl;
02582                         std::cout<<"cos sin of alpha_new==="<<cos_alpha<<sin_alpha<<std::endl;
02583                         std::cout<<"alpah dig==="<<cos_alpha<<sin_alpha<<std::endl;
02584                         std::cout<<"prjection maximum==="<<loop_range*ellipse_step<<"ideal maximum"<<m_sizeofprojection*m_npad/2<<std::endl;
02585                         std::cout<<"x_size=="<<m_volume->get_xsize()<<"y_size=="<<m_volume->get_ysize()<<std::endl;
02586                         std::cout<<"#############################################################"<<std::endl;
02587 
02588 
02589                         
02590 
02591                 }
02592                 for(int i=0; i <=loop_range; ++i ) {
02593                         float xnew = i*cos_alpha;
02594                         float ynew = -i*sin_alpha;
02595                         if(m_count%100==0&&i==loop_range)
02596                                 std::cout<<"x_new=="<<xnew<<"Y_new=="<<ynew<<std::endl;
02597                         float btqr=0,btqi=0;
02598                         float xprj=i*ellipse_step;
02599                         float t=xprj-int(xprj);
02600                         btqr = (1-t)*padfft->get_value_at( 2*int(xprj), 0, 0 )+t*padfft->get_value_at( 2*(1+int(xprj)), 0, 0 );
02601                         btqi = (1-t)*padfft->get_value_at( 2*int(xprj)+1, 0, 0 )+t*padfft->get_value_at( 2*(1+int(xprj))+1, 0, 0 );
02602                         if( xnew < 0.0 ) {
02603                                 xnew *= -1;
02604                                 ynew *= -1;
02605                                 btqi *= -1;
02606                         }
02607 
02608                         int ixn = int(xnew+0.5+m_vnxp) - m_vnxp;
02609                         int iyn = int(ynew+0.5+m_vnyp) - m_vnyp;
02610 
02611                         if(iyn < 0 ) iyn += m_vnyp;
02612                         if(m_count%100==0&&i==loop_range)
02613                                 std::cout<<"xnn=="<<ixn<<"ynn=="<<iyn<<std::endl;
02614                         (*m_volume)( 2*ixn, iyn+1, 1 ) += btqr *float(mult);
02615                         (*m_volume)( 2*ixn+1, iyn+1, 1 ) += btqi * float(mult);
02616                         (*m_wptr)(ixn,iyn+1, 1) += float(mult);
02617                 }
02618 
02619 
02620         }
02621         checked_delete( padfft );
02622         return 0;
02623 }
02624 
02625 
02626 
02627 
02628 int nn4_rectReconstructor::insert_padfft_slice( EMData* padded, const Transform& t, int mult )
02629 {
02630         Assert( padded != NULL );
02631         
02632         vector<Transform> tsym = t.get_sym_proj(m_symmetry);
02633         for (unsigned int isym=0; isym < tsym.size(); isym++) {
02634                 m_volume->insert_rect_slice(m_wptr, padded, tsym[isym], m_sizeofprojection, m_xratio, m_yratio, m_npad, mult);
02635         }
02636         
02637 
02638         return 0;
02639 
02640 }
02641 
02642 #define  tw(i,j,k)      tw[ i-1 + (j-1+(k-1)*iy)*ix ]
02643 void circumf_rect( EMData* win , int npad)
02644 {
02645         float *tw = win->get_data();
02646         //  correct for the fall-off
02647         //  mask and subtract circumference average
02648         int ix = win->get_xsize();
02649         int iy = win->get_ysize();
02650         int iz = win->get_zsize();
02651 
02652         int IP = ix/2+1;
02653         int JP = iy/2+1;
02654         int KP = iz/2+1;
02655         
02656         //  sinc functions tabulated for fall-off
02657         float* sincx = new float[IP+1];
02658         float* sincy = new float[JP+1];
02659         float* sincz = new float[KP+1];
02660 
02661         sincx[0] = 1.0f;
02662         sincy[0] = 1.0f;
02663         sincz[0] = 1.0f;
02664 
02665         float cdf = M_PI/float(npad*2*ix);
02666         for (int i = 1; i <= IP; ++i)  sincx[i] = sin(i*cdf)/(i*cdf);
02667         cdf = M_PI/float(npad*2*iy);
02668         for (int i = 1; i <= JP; ++i)  sincy[i] = sin(i*cdf)/(i*cdf);
02669         cdf = M_PI/float(npad*2*iz);
02670         for (int i = 1; i <= KP; ++i)  sincz[i] = sin(i*cdf)/(i*cdf);
02671         for (int k = 1; k <= iz; ++k) {
02672                 int kkp = abs(k-KP);
02673                 for (int j = 1; j <= iy; ++j) {
02674                         cdf = sincy[abs(j- JP)]*sincz[kkp];
02675                         for (int i = 1; i <= ix; ++i)  tw(i,j,k) /= (sincx[abs(i-IP)]*cdf);
02676                 }
02677         }
02678 
02679         delete[] sincx;
02680         delete[] sincy;
02681         delete[] sincz;
02682         
02683         
02684         
02685         float dxx = 1.0f/float(0.25*ix*ix);
02686         float dyy = 1.0f/float(0.25*iy*iy);
02687         
02688         
02689         
02690         float LR2=(float(ix)/2-1)*(float(ix)/2-1)*dxx;
02691 
02692         float  TNR = 0.0f;
02693         size_t m = 0;
02694         for (int k = 1; k <= iz; ++k) {
02695                 for (int j = 1; j <= iy; ++j) {
02696                         for (int i = 1; i <= ix; ++i) {
02697                                 float LR = (j-JP)*(j-JP)*dyy+(i-IP)*(i-IP)*dxx;
02698                                 if (LR<=1.0f && LR >= LR2) {
02699                                         TNR += tw(i,j,k);
02700                                         ++m;
02701                                 }
02702                         }
02703                 }
02704         }
02705 
02706         TNR /=float(m);
02707         
02708 
02709         for (int k = 1; k <= iz; ++k) {
02710                 for (int j = 1; j <= iy; ++j) {
02711                         for (int i = 1; i <= ix; ++i) {
02712                                 float LR = (j-JP)*(j-JP)*dyy+(i-IP)*(i-IP)*dxx;
02713                                 if (LR<=1.0f)  tw(i,j,k)=tw(i,j,k)-TNR;
02714                                 else            tw(i,j,k) = 0.0f;       
02715                         }
02716                 }
02717         }
02718 
02719 }
02720 #undef tw 
02721 
02722 EMData* nn4_rectReconstructor::finish(bool)
02723 {
02724         
02725         if( m_ndim==3 ) {
02726                 m_volume->symplane0_rect(m_wptr);
02727         } else {
02728                 for( int i=1; i <= m_vnyp; ++i ) {
02729 
02730                         if( (*m_wptr)(0, i, 1)==0.0 ) {
02731                                 int j = m_vnyp + 1 - i;
02732                                 (*m_wptr)(0, i, 1) = (*m_wptr)(0, j, 1);
02733                                 (*m_volume)(0, i, 1) = (*m_volume)(0, j, 1);
02734                                 (*m_volume)(1, i, 1) = (*m_volume)(1, j, 1);
02735                         }
02736                 }
02737         }
02738 
02739 
02740         int box = 7;
02741         int kc = (box-1)/2;
02742         vector< float > pow_a( m_ndim*kc+1, 1.0 );
02743         for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
02744         pow_a.back()=0.0;
02745 
02746         float alpha = 0.0;
02747         if( m_ndim==3) {
02748                 int vol = box*box*box;
02749                 float max = max3d( kc, pow_a );
02750                 alpha = ( 1.0f - 1.0f/(float)vol ) / max;
02751         } else {
02752                 int ara = box*box;
02753                 float max = max2d( kc, pow_a );
02754                 alpha = ( 1.0f - 1.0f/(float)ara ) / max;
02755         }
02756 
02757         int ix,iy,iz;
02758         for (iz = 1; iz <= m_vnzp; iz++) {
02759                 for (iy = 1; iy <= m_vnyp; iy++) {
02760                         for (ix = 0; ix <= m_vnxc; ix++) {
02761                                 if ( (*m_wptr)(ix,iy,iz) > 0) {//(*v) should be treated as complex!!
02762                                         float tmp;
02763                                         tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+m_osnr);
02764                                         
02765                                         if( m_weighting == ESTIMATE ) {
02766                                                 int cx = ix;
02767                                                 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
02768                                                 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
02769                                                 float sum = 0.0;
02770                                                 for( int ii = -kc; ii <= kc; ++ii ) {
02771                                                         int nbrcx = cx + ii;
02772                                                         if( nbrcx >= m_vnxc ) continue;
02773                                                         for( int jj= -kc; jj <= kc; ++jj ) {
02774                                                                 int nbrcy = cy + jj;
02775                                                                 if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
02776 
02777                                                                 int kcz = (m_ndim==3) ? kc : 0;
02778                                                                 for( int kk = -kcz; kk <= kcz; ++kk ) {
02779                                                                         int nbrcz = cz + kk;
02780                                                                         if( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
02781                                                                         if( nbrcx < 0 ) {
02782                                                                                 nbrcx = -nbrcx;
02783                                                                                 nbrcy = -nbrcy;
02784                                                                                 nbrcz = -nbrcz;
02785                                                                         }
02786                                                                         int nbrix = nbrcx;
02787                                                                         int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
02788                                                                         int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
02789                                                                         if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0 ) {
02790                                                                                 int c = m_ndim*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
02791                                                                                 sum = sum + pow_a[c];
02792                                                                         }
02793                                                                 }
02794                                                         }
02795                                                 }
02796                                                 float wght = 1.0f / ( 1.0f - alpha * sum );
02797                                                 tmp = tmp * wght;
02798                                         }
02799                                         (*m_volume)(2*ix,iy,iz)   *= tmp;
02800                                         (*m_volume)(2*ix+1,iy,iz) *= tmp;
02801                                 }
02802                         }
02803                 }
02804         }
02805 
02806         //if(m_ndim==2) printImage( m_volume );
02807 
02808         // back fft
02809         m_volume->do_ift_inplace();
02810 
02811         
02812         int npad = m_volume->get_attr("npad");
02813         m_volume->depad();
02814         circumf_rect( m_volume, npad );
02815         //circumf_rect_new( m_volume, npad,m_xratio,m_yratio);
02816         m_volume->set_array_offsets( 0, 0, 0 );
02817 
02818         return 0;
02819 }
02820 
02821 
02822 // Added By Zhengfan Yang on 03/16/07
02823 // Beginning of the addition
02824 // --------------------------------------------------------------------------------
02825 
02826 nnSSNR_Reconstructor::nnSSNR_Reconstructor()
02827 {
02828         m_volume = NULL;
02829         m_wptr   = NULL;
02830         m_wptr2  = NULL;
02831 }
02832 
02833 nnSSNR_Reconstructor::nnSSNR_Reconstructor( const string& symmetry, int size, int npad)
02834 {
02835         m_volume = NULL;
02836         m_wptr   = NULL;
02837         m_wptr2  = NULL;
02838 
02839         setup( symmetry, size, npad );
02840 }
02841 
02842 nnSSNR_Reconstructor::~nnSSNR_Reconstructor()
02843 {
02844         //if( m_delete_volume ) checked_delete(m_volume);
02845 
02846         //if( m_delete_weight ) checked_delete( m_wptr );
02847 
02848         //if( m_delete_weight2 ) checked_delete( m_wptr2 );
02849 
02850         //checked_delete( m_result );
02851 }
02852 
02853 void nnSSNR_Reconstructor::setup()
02854 {
02855         int size = params["size"];
02856         int npad = params["npad"];
02857 
02858         string symmetry;
02859         if( params.has_key("symmetry") ) symmetry = params["symmetry"].to_str();
02860         else                             symmetry = "c1";
02861 
02862         setup( symmetry, size, npad );
02863 }
02864 
02865 void nnSSNR_Reconstructor::setup( const string& symmetry, int size, int npad )
02866 {
02867 
02868         m_weighting = ESTIMATE;
02869         m_wghta = 0.2f;
02870         m_wghtb = 0.004f;
02871 
02872         m_symmetry = symmetry;
02873         m_npad = npad;
02874         m_nsym = Transform::get_nsym(m_symmetry);
02875 
02876         m_vnx = size;
02877         m_vny = size;
02878         m_vnz = size;
02879 
02880         m_vnxp = size*npad;
02881         m_vnyp = size*npad;
02882         m_vnzp = size*npad;
02883 
02884         m_vnxc = m_vnxp/2;
02885         m_vnyc = m_vnyp/2;
02886         m_vnzc = m_vnzp/2;
02887 
02888         buildFFTVolume();
02889         buildNormVolume();
02890         buildNorm2Volume();
02891 
02892 }
02893 
02894 void nnSSNR_Reconstructor::buildFFTVolume() {
02895 
02896         m_volume = params["fftvol"]; 
02897         m_volume->set_size(m_vnxp+ 2 - m_vnxp%2,m_vnyp,m_vnzp);
02898         m_volume->to_zero();
02899         if ( m_vnxp % 2 == 0 ) m_volume->set_fftodd(0);
02900         else                   m_volume->set_fftodd(1);
02901 
02902         m_volume->set_nxc(m_vnxc);
02903         m_volume->set_complex(true);
02904         m_volume->set_ri(true);
02905         m_volume->set_fftpad(true);
02906         m_volume->set_attr("npad", m_npad);
02907         m_volume->set_array_offsets(0,1,1);
02908 }
02909 
02910 void nnSSNR_Reconstructor::buildNormVolume() {
02911 
02912         m_wptr = params["weight"];
02913         m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
02914         m_wptr->to_zero();
02915         m_wptr->set_array_offsets(0,1,1);
02916 }
02917 
02918 void nnSSNR_Reconstructor::buildNorm2Volume() {
02919 
02920         m_wptr2 = params["weight2"];
02921         m_wptr2->set_size(m_vnxc+1,m_vnyp,m_vnzp);
02922         m_wptr2->to_zero();
02923         m_wptr2->set_array_offsets(0,1,1);
02924 }
02925 
02926 
02927 int nnSSNR_Reconstructor::insert_slice(const EMData* const slice, const Transform& t, const float) {
02928         // sanity checks
02929         if (!slice) {
02930                 LOGERR("try to insert NULL slice");
02931                 return 1;
02932         }
02933 
02934         int padffted=slice->get_attr_default( "padffted", 0 );
02935 
02936         if ( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_vnx)  ) {
02937                 // FIXME: Why doesn't this throw an exception?
02938                 LOGERR("Tried to insert a slice that has wrong size.");
02939                 return 1;
02940         }
02941 
02942         EMData* padfft = NULL;
02943 
02944         if( padffted != 0 ) padfft = new EMData(*slice);
02945         else                padfft = padfft_slice( slice, t, m_npad );
02946 
02947         int mult = slice->get_attr_default("mult", 1);
02948 
02949         Assert( mult > 0 );
02950         insert_padfft_slice( padfft, t, mult );
02951 
02952         checked_delete( padfft );
02953         return 0;
02954 }
02955 
02956 int nnSSNR_Reconstructor::insert_padfft_slice( EMData* padfft, const Transform& t, int mult )
02957 {
02958         Assert( padfft != NULL );
02959         // insert slice for all symmetry related positions
02960         for (int isym=0; isym < m_nsym; isym++) {
02961                 Transform tsym = t.get_sym(m_symmetry, isym);
02962                 m_volume->nn_SSNR( m_wptr, m_wptr2, padfft, tsym, mult);
02963         }
02964         return 0;
02965 }
02966 
02967 
02968 #define  tw(i,j,k)      tw[ i-1 + (j-1+(k-1)*iy)*ix ]
02969 EMData* nnSSNR_Reconstructor::finish(bool)
02970 {
02971 /*
02972   I changed the code on 05/15 so it only returns variance.
02973   Lines commented out are marked by //#
02974   The version prior to the currect changes is r1.190
02975   PAP
02976 */
02977         int kz, ky;
02978         //#int iix, iiy, iiz;
02979         int box = 7;
02980         int kc = (box-1)/2;
02981         float alpha = 0.0;
02982         float argx, argy, argz;
02983         vector< float > pow_a( 3*kc+1, 1.0 );
02984         float w = params["w"];
02985         EMData* SSNR = params["SSNR"];
02986         //#EMData* vol_ssnr = new EMData();
02987         //#vol_ssnr->set_size(m_vnxp, m_vnyp, m_vnzp);
02988         //#vol_ssnr->to_zero();
02989         //#  new line follow
02990         EMData* vol_ssnr = params["vol_ssnr"];
02991         vol_ssnr->set_size(m_vnxp+ 2 - m_vnxp%2, m_vnyp ,m_vnzp);
02992         vol_ssnr->to_zero();
02993         if ( m_vnxp % 2 == 0 ) vol_ssnr->set_fftodd(0);
02994         else                   vol_ssnr->set_fftodd(1);
02995         vol_ssnr->set_nxc(m_vnxc);
02996         vol_ssnr->set_complex(true);
02997         vol_ssnr->set_ri(true);
02998         vol_ssnr->set_fftpad(false);
02999         //#
03000 
03001         float dx2 = 1.0f/float(m_vnxc)/float(m_vnxc);
03002         float dy2 = 1.0f/float(m_vnyc)/float(m_vnyc);
03003 #ifdef _WIN32
03004         float dz2 = 1.0f/_cpp_max(float(m_vnzc),1.0f)/_cpp_max(float(m_vnzc),1.0f);
03005         int  inc = Util::round(float(_cpp_max(_cpp_max(m_vnxc,m_vnyc),m_vnzc))/w);
03006 #else
03007         float dz2 = 1.0f/std::max(float(m_vnzc),1.0f)/std::max(float(m_vnzc),1.0f);
03008         int  inc = Util::round(float(std::max(std::max(m_vnxc,m_vnyc),m_vnzc))/w);
03009 #endif  //_WIN32
03010         SSNR->set_size(inc+1,4,1);
03011 
03012         float *nom    = new float[inc+1];
03013         float *denom  = new float[inc+1];
03014         int  *nn     = new int[inc+1];
03015         int  *ka     = new int[inc+1];
03016         float wght = 1.0f;
03017         for (int i = 0; i <= inc; i++) {
03018                 nom[i] = 0.0f;
03019                 denom[i] = 0.0f;
03020                 nn[i] = 0;
03021                 ka[i] = 0;
03022         }
03023 
03024         m_volume->symplane1(m_wptr, m_wptr2);
03025 
03026         if ( m_weighting == ESTIMATE ) {
03027                 int vol = box*box*box;
03028                 for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
03029                 pow_a[3*kc] = 0.0;
03030                 float max = max3d( kc, pow_a );
03031                 alpha = ( 1.0f - 1.0f/(float)vol ) / max;
03032         }
03033 
03034         for (int iz = 1; iz <= m_vnzp; iz++) {
03035                 if ( iz-1 > m_vnzc ) kz = iz-1-m_vnzp; else kz = iz-1;
03036                 argz = float(kz*kz)*dz2;
03037                 for (int iy = 1; iy <= m_vnyp; iy++) {
03038                         if ( iy-1 > m_vnyc ) ky = iy-1-m_vnyp; else ky = iy-1;
03039                         argy = argz + float(ky*ky)*dy2;
03040                         for (int ix = 0; ix <= m_vnxc; ix++) {
03041                                 float Kn = (*m_wptr)(ix,iy,iz);
03042                                 argx = std::sqrt(argy + float(ix*ix)*dx2);
03043                                 int r = Util::round(float(inc)*argx);
03044                                 if ( r >= 0 && Kn > 4.5f ) {
03045                                         if ( m_weighting == ESTIMATE ) {
03046                                                 int cx = ix;
03047                                                 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
03048                                                 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
03049 
03050                                                 float sum = 0.0;
03051                                                 for( int ii = -kc; ii <= kc; ++ii ) {
03052                                                         int nbrcx = cx + ii;
03053                                                         if( nbrcx >= m_vnxc ) continue;
03054                                                     for ( int jj= -kc; jj <= kc; ++jj ) {
03055                                                                 int nbrcy = cy + jj;
03056                                                                 if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
03057                                                                 for( int kk = -kc; kk <= kc; ++kk ) {
03058                                                                         int nbrcz = cz + jj;
03059                                                                         if ( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
03060                                                                         if( nbrcx < 0 ) {
03061                                                                                 nbrcx = -nbrcx;
03062                                                                                 nbrcy = -nbrcy;
03063                                                                                 nbrcz = -nbrcz;
03064                                                                         }
03065                                                                         int nbrix = nbrcx;
03066                                                                         int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
03067                                                                         int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
03068                                                                         if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0 ) {
03069                                                                                 int c = 3*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
03070                                                                                 sum = sum + pow_a[c];
03071                                                                         }
03072                                                                 }
03073                                                         }
03074                                                 }
03075                                                 wght = 1.0f / ( 1.0f - alpha * sum );
03076                                         } // end of ( m_weighting == ESTIMATE )
03077                                         float nominator = std::norm(m_volume->cmplx(ix,iy,iz))/Kn;
03078                                         float denominator = ((*m_wptr2)(ix,iy,iz)-nominator)/(Kn-1.0f);
03079                                         // Skip Friedel related values
03080                                         if( (ix>0 || (kz>=0 && (ky>=0 || kz!=0)))) {
03081                                                 if ( r <= inc ) {
03082                                                         nom[r]   += nominator*wght;
03083                                                         denom[r] += denominator/Kn*wght;
03084                                                         nn[r]    += 2;
03085                                                         ka[r]    += int(Kn);
03086                                                 }
03087 /*
03088 #ifdef  _WIN32
03089                                                 //#float  tmp = _cpp_max(nominator/denominator/Kn-1.0f,0.0f);
03090 #else
03091                                                 //#float  tmp = std::max(nominator/denominator/Kn-1.0f,0.0f);
03092 #endif  //_WIN32
03093                                                 //  Create SSNR as a 3D array (-n/2:n/2+n%2-1)
03094                                                 iix = m_vnxc + ix; iiy = m_vnyc + ky; iiz = m_vnzc + kz;
03095                                                 if( iix >= 0 && iix < m_vnxp && iiy >= 0 && iiy < m_vnyp && iiz >= 0 && iiz < m_vnzp )
03096                                                         (*vol_ssnr)(iix, iiy, iiz) = tmp;
03097                                                 // Friedel part
03098                                                 iix = m_vnxc - ix; iiy = m_vnyc - ky; iiz = m_vnzc - kz;
03099                                                 if( iix >= 0 && iix < m_vnxp && iiy >= 0 && iiy < m_vnyp && iiz >= 0 && iiz < m_vnzp )
03100                                                         (*vol_ssnr)(iix, iiy, iiz) = tmp;
03101 */
03102 
03103                                         }
03104                                         (*vol_ssnr)(2*ix, iy-1, iz-1) = denominator*wght;
03105                                         //(*vol_ssnr)(2*ix, iy-1, iz-1) =  real(m_volume->cmplx(ix,iy,iz))*wght/Kn;
03106                                         //(*vol_ssnr)(2*ix+1, iy-1, iz-1) = imag(m_volume->cmplx(ix,iy,iz))*wght/Kn;
03107                                 } // end of Kn>4.5
03108                         }
03109                 }
03110         }
03111 
03112         for (int i = 0; i <= inc; i++)  {
03113                 (*SSNR)(i,0,0) = nom[i];  
03114                 (*SSNR)(i,1,0) = denom[i];    // variance
03115                 (*SSNR)(i,2,0) = static_cast<float>(nn[i]);
03116                 (*SSNR)(i,3,0) = static_cast<float>(ka[i]);
03117         }
03118         vol_ssnr->update();
03119         
03120         delete[] nom;
03121         delete[] denom;
03122         delete[] nn;
03123         delete[] ka;
03124 
03125         return 0;
03126 }
03127 #undef  tw
03128 
03129 // -----------------------------------------------------------------------------------
03130 // End of this addition
03131 
03132 //####################################################################################
03133 //** nn4 ctf reconstructor
03134 
03135 nn4_ctfReconstructor::nn4_ctfReconstructor()
03136 {
03137         m_volume  = NULL;
03138         m_wptr    = NULL;
03139 }
03140 
03141 nn4_ctfReconstructor::nn4_ctfReconstructor( const string& symmetry, int size, int npad, float snr, int sign )
03142 {
03143         setup( symmetry, size, npad, snr, sign );
03144 }
03145 
03146 nn4_ctfReconstructor::~nn4_ctfReconstructor()
03147 {
03148         //if( m_delete_volume ) checked_delete(m_volume);
03149 
03150         //if( m_delete_weight ) checked_delete( m_wptr );
03151 
03152         //checked_delete( m_result );
03153 }
03154 
03155 void nn4_ctfReconstructor::setup()
03156 {
03157         if( ! params.has_key("size") ) throw std::logic_error("Error: image size is not given");
03158 
03159         int size = params["size"];
03160         int npad = params.has_key("npad") ? int(params["npad"]) : 4;
03161         // int sign = params.has_key("sign") ? int(params["sign"]) : 1;
03162         int sign = 1;
03163         string symmetry = params.has_key("symmetry")? params["symmetry"].to_str() : "c1";
03164 
03165         float snr = params["snr"];
03166 
03167         m_varsnr = params.has_key("varsnr") ? int(params["varsnr"]) : 0;
03168         setup( symmetry, size, npad, snr, sign );
03169 
03170 }
03171 
03172 void nn4_ctfReconstructor::setup( const string& symmetry, int size, int npad, float snr, int sign )
03173 {
03174         m_weighting = ESTIMATE;
03175         if( params.has_key("weighting") ) {
03176                 int tmp = int( params["weighting"] );
03177                 if( tmp==0 ) m_weighting = NONE;
03178         }
03179 
03180 
03181 
03182         m_wghta = 0.2f;
03183         m_wghtb = 0.004f;
03184 
03185         m_symmetry = symmetry;
03186         m_npad = npad;
03187         m_sign = sign;
03188         m_nsym = Transform::get_nsym(m_symmetry);
03189 
03190         m_snr = snr;
03191 
03192         m_vnx = size;
03193         m_vny = size;
03194         m_vnz = size;
03195 
03196         m_vnxp = size*npad;
03197         m_vnyp = size*npad;
03198         m_vnzp = size*npad;
03199 
03200         m_vnxc = m_vnxp/2;
03201         m_vnyc = m_vnyp/2;
03202         m_vnzc = m_vnzp/2;
03203 
03204         buildFFTVolume();
03205         buildNormVolume();
03206 }
03207 
03208 void nn4_ctfReconstructor::buildFFTVolume() {
03209         int offset = 2 - m_vnxp%2;
03210 
03211         m_volume = params["fftvol"];
03212 
03213         if( m_volume->get_xsize() != m_vnxp+offset && m_volume->get_ysize() != m_vnyp && m_volume->get_zsize() != m_vnzp ) {
03214                 m_volume->set_size(m_vnxp+offset,m_vnyp,m_vnzp);
03215                 m_volume->to_zero();
03216         }
03217 
03218         m_volume->set_nxc(m_vnxp/2);
03219         m_volume->set_complex(true);
03220         m_volume->set_ri(true);
03221         m_volume->set_fftpad(true);
03222         m_volume->set_attr("npad", m_npad);
03223         m_volume->set_array_offsets(0,1,1);
03224 }
03225 
03226 void nn4_ctfReconstructor::buildNormVolume()
03227 {
03228         m_wptr = params["weight"];
03229 
03230         if( m_wptr->get_xsize() != m_vnxc+1 && m_wptr->get_ysize() != m_vnyp && m_wptr->get_zsize() != m_vnzp ) {
03231                m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
03232                m_wptr->to_zero();
03233         }
03234 
03235         m_wptr->set_array_offsets(0,1,1);
03236 
03237 }
03238 
03239 int nn4_ctfReconstructor::insert_slice(const EMData* const slice, const Transform& t, const float)
03240 {
03241         // sanity checks
03242         if (!slice) {
03243                 LOGERR("try to insert NULL slice");
03244                 return 1;
03245         }
03246 
03247         int buffed = slice->get_attr_default( "buffed", 0 );
03248         if( buffed > 0 ) {
03249                 int mult = slice->get_attr_default( "mult", 1 );
03250                 insert_buffed_slice( slice, mult );
03251                 return 0;
03252         }
03253 
03254         int padffted= slice->get_attr_default("padffted", 0);
03255         if( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_vnx)  )
03256         {
03257                 // FIXME: Why doesn't this throw an exception?
03258                 LOGERR("Tried to insert a slice that is the wrong size.");
03259                 return 1;
03260         }
03261 
03262         EMData* padfft = NULL;
03263 
03264         if( padffted != 0 ) padfft = new EMData(*slice);
03265         else                padfft = padfft_slice( slice, t, m_npad );
03266 
03267         int mult= slice->get_attr_default("mult", 1);
03268 
03269         Assert( mult > 0 );
03270         insert_padfft_slice( padfft, t, mult );
03271 
03272         checked_delete( padfft );
03273 
03274         return 0;
03275 }
03276 
03277 int nn4_ctfReconstructor::insert_buffed_slice( const EMData* buffed, int mult )
03278 {
03279         const float* bufdata = buffed->get_data();
03280         float* cdata = m_volume->get_data();
03281         float* wdata = m_wptr->get_data();
03282 
03283         int npoint = buffed->get_xsize()/4;
03284         for( int i=0; i < npoint; ++i ) {
03285 
03286                 int pos2 = int( bufdata[4*i] );
03287                 int pos1 = pos2 * 2;
03288                 cdata[pos1  ] += bufdata[4*i+1]*mult;
03289                 cdata[pos1+1] += bufdata[4*i+2]*mult;
03290                 wdata[pos2  ] += bufdata[4*i+3]*mult;
03291 /*
03292         std::cout << "pos1, pos2, ctfv1, ctfv2, ctf2: ";
03293         std::cout << pos1 << " " << bufdata[5*i+1] << " " << bufdata[5*i+2] << " ";
03294         std::cout << pos2 << " " << bufdata[5*i+4] << std::endl;
03295  */
03296         }
03297         return 0;
03298 }
03299 
03300 int nn4_ctfReconstructor::insert_padfft_slice( EMData* padfft, const Transform& t, int mult )
03301 {
03302         Assert( padfft != NULL );
03303         float tmp = padfft->get_attr("ctf_applied");
03304         int   ctf_applied = (int) tmp;
03305         
03306         vector<Transform> tsym = t.get_sym_proj(m_symmetry);
03307         for (unsigned int isym=0; isym < tsym.size(); isym++) {
03308                 if(ctf_applied) m_volume->nn_ctf_applied(m_wptr, padfft, tsym[isym], mult);
03309                                 
03310                 else            m_volume->nn_ctf(m_wptr, padfft, tsym[isym], mult);
03311                 
03312         }
03313 
03314 
03315         return 0;
03316 
03317 }
03318 
03319 #define  tw(i,j,k)      tw[ i-1 + (j-1+(k-1)*iy)*ix ]
03320 EMData* nn4_ctfReconstructor::finish(bool)
03321 {
03322         m_volume->set_array_offsets(0, 1, 1);
03323         m_wptr->set_array_offsets(0, 1, 1);
03324         m_volume->symplane0_ctf(m_wptr);
03325 
03326         int box = 7;
03327         int vol = box*box*box;
03328         int kc = (box-1)/2;
03329         vector< float > pow_a( 3*kc+1, 1.0 );
03330         for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
03331         pow_a[3*kc]=0.0;
03332 
03333 
03334         float max = max3d( kc, pow_a );
03335         float alpha = ( 1.0f - 1.0f/(float)vol ) / max;
03336         float osnr = 1.0f/m_snr;
03337 
03338         // normalize
03339         int ix,iy,iz;
03340         for (iz = 1; iz <= m_vnzp; iz++) {
03341                 for (iy = 1; iy <= m_vnyp; iy++) {
03342                         for (ix = 0; ix <= m_vnxc; ix++) {
03343                                 if ( (*m_wptr)(ix,iy,iz) > 0.0f) {//(*v) should be treated as complex!!
03344                     int iyp = (iy<=m_vnyc) ? iy - 1 : iy-m_vnyp-1;
03345                     int izp = (iz<=m_vnzc) ? iz - 1 : iz-m_vnzp-1;
03346                     float tmp=0.0;
03347                     if( m_varsnr )
03348                     {
03349                                             float freq = sqrt( (float)(ix*ix+iyp*iyp+izp*izp) );
03350                         tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+freq*osnr)*m_sign;
03351                     }
03352                     else
03353                     {
03354                         tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+osnr)*m_sign;
03355                     }
03356 
03357                                         if( m_weighting == ESTIMATE ) {
03358                                                 int cx = ix;
03359                                                 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
03360                                                 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
03361                                                 float sum = 0.0;
03362                                                 for( int ii = -kc; ii <= kc; ++ii ) {
03363                                                         int nbrcx = cx + ii;
03364                                                         if( nbrcx >= m_vnxc ) continue;
03365                                                         for( int jj= -kc; jj <= kc; ++jj ) {
03366                                                                 int nbrcy = cy + jj;
03367                                                                 if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
03368                                                                 for( int kk = -kc; kk <= kc; ++kk ) {
03369                                                                         int nbrcz = cz + jj;
03370                                                                         if( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
03371                                                                         if( nbrcx < 0 ) {
03372                                                                                 nbrcx = -nbrcx;
03373                                                                                 nbrcy = -nbrcy;
03374                                                                                 nbrcz = -nbrcz;
03375                                                                         }
03376 
03377                                                                         int nbrix = nbrcx;
03378                                                                         int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
03379                                                                         int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
03380                                                                         if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0.0 ) {
03381                                                                                 int c = 3*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
03382                                                                                 sum = sum + pow_a[c];
03383                                                                                   // if(ix%20==0 && iy%20==0 && iz%20==0)
03384                                                                                  //   std::cout << boost::format( "%4d %4d %4d %4d %10.3f" ) % nbrix % nbriy % nbriz % c % sum << std::endl;
03385                                                                         }
03386                                                                 }
03387                                                         }
03388                                                 }
03389                                                 float wght = 1.0f / ( 1.0f - alpha * sum );
03390 /*
03391                         if(ix%10==0 && iy%10==0)
03392                         {
03393                             std::cout << boost::format( "%4d %4d %4d " ) % ix % iy %iz;
03394                             std::cout << boost::format( "%10.3f %10.3f %10.3f " )  % tmp % wght % sum;
03395                             std::  << boost::format( "%10.3f %10.3e " ) % pow_b[r] % alpha;
03396                             std::cout << std::endl;
03397                         }
03398  */
03399                                                 tmp = tmp * wght;
03400                                         }
03401                                         (*m_volume)(2*ix,iy,iz) *= tmp;
03402                                         (*m_volume)(2*ix+1,iy,iz) *= tmp;
03403                                 }
03404                         }
03405                 }
03406         }
03407 
03408         // back fft
03409         m_volume->do_ift_inplace();
03410         int npad = m_volume->get_attr("npad");
03411         m_volume->depad();
03412         circumf( m_volume, npad );
03413         m_volume->set_array_offsets( 0, 0, 0 );
03414 
03415         return 0;
03416 }
03417 #undef  tw
03418 
03419 
03420 
03421 //####################################################################################
03422 //** nn4 ctf rect reconstructor
03423 
03424 nn4_ctf_rectReconstructor::nn4_ctf_rectReconstructor()
03425 {
03426         m_volume  = NULL;
03427         m_wptr    = NULL;
03428 }
03429 
03430 nn4_ctf_rectReconstructor::nn4_ctf_rectReconstructor( const string& symmetry, int size, int npad, float snr, int sign )
03431 {
03432         setup( symmetry, size, npad, snr, sign );
03433 }
03434 
03435 nn4_ctf_rectReconstructor::~nn4_ctf_rectReconstructor()
03436 {
03437         //if( m_delete_volume ) checked_delete(m_volume);
03438 
03439         //if( m_delete_weight ) checked_delete( m_wptr );
03440 
03441         //checked_delete( m_result );
03442 }
03443 
03444 void nn4_ctf_rectReconstructor::setup()
03445 {
03446         if( ! params.has_key("sizeprojection") ) throw std::logic_error("Error: projection size is not given");
03447         m_sizeofprojection = params["sizeprojection"];
03448         int npad = params.has_key("npad") ? int(params["npad"]) : 4;
03449         // int sign = params.has_key("sign") ? int(params["sign"]) : 1;
03450         int sign = 1;
03451         string symmetry = params.has_key("symmetry")? params["symmetry"].to_str() : "c1";
03452 
03453         float snr = params["snr"];
03454 
03455         m_varsnr = params.has_key("varsnr") ? int(params["varsnr"]) : 0;
03456         setup( symmetry, m_sizeofprojection, npad, snr, sign );
03457 
03458 }
03459 
03460 void nn4_ctf_rectReconstructor::setup( const string& symmetry, int sizeprojection, int npad, float snr, int sign )
03461 {
03462         m_weighting = ESTIMATE;
03463         if( params.has_key("weighting") ) {
03464                 int tmp = int( params["weighting"] );
03465                 if( tmp==0 ) m_weighting = NONE;
03466         }
03467 
03468         m_wghta = 0.2f;
03469         m_wghtb = 0.004f;
03470 
03471         m_symmetry = symmetry;
03472         m_npad = npad;
03473         m_sign = sign;
03474         m_nsym = Transform::get_nsym(m_symmetry);
03475 
03476         m_snr = snr;
03477         if (params.has_key("sizex"))  m_vnx = params["sizex"];
03478         else if (params.has_key("xratio")) {
03479                 float temp=params["xratio"];
03480                 m_vnx=int(float(sizeprojection)*temp);
03481         } else  m_vnx=sizeprojection;
03482 
03483         if (params.has_key("sizey"))  m_vny = params["sizey"];
03484         else if (params.has_key("yratio")) {
03485                 float temp=params["yratio"];
03486                 m_vny=int(float(sizeprojection)*temp);
03487         }
03488         else m_vny=sizeprojection;
03489 
03490         if( params.has_key("sizez") ) 
03491                 m_vnz = params["sizez"];
03492         else                          
03493                 m_vnz =sizeprojection;
03494         
03495         m_xratio=float(m_vnx)/float(sizeprojection);    
03496         m_yratio=float(m_vny)/float(sizeprojection);
03497 
03498         //std::cout<<"xratio=="<<m_xratio<<"yratio=="<<m_yratio<<std::endl;
03499         //std::cout<<"sx=="<<m_vnx<<"sy=="<<m_vny<<"sz=="<<m_vnz<<std::endl;
03500 
03501         m_vnxp = m_vnx*npad;
03502         m_vnyp = m_vny*npad;
03503         m_vnzp = m_vnz*npad;
03504 
03505         m_vnxc = m_vnxp/2;
03506         m_vnyc = m_vnyp/2;
03507         m_vnzc = m_vnzp/2;
03508 
03509         buildFFTVolume();
03510         buildNormVolume();
03511 }
03512 
03513 void nn4_ctf_rectReconstructor::buildFFTVolume() {
03514         int offset = 2 - m_vnxp%2;
03515 
03516         m_volume = params["fftvol"];
03517 
03518         if( m_volume->get_xsize() != m_vnxp+offset && m_volume->get_ysize() != m_vnyp && m_volume->get_zsize() != m_vnzp ) {
03519                 m_volume->set_size(m_vnxp+offset,m_vnyp,m_vnzp);
03520                 m_volume->to_zero();
03521         }
03522 
03523         m_volume->set_nxc(m_vnxp/2);
03524         m_volume->set_complex(true);
03525         m_volume->set_ri(true);
03526         m_volume->set_fftpad(true);
03527         m_volume->set_attr("npad", m_npad);
03528         m_volume->set_array_offsets(0,1,1);
03529 }
03530 
03531 void nn4_ctf_rectReconstructor::buildNormVolume()
03532 {
03533         m_wptr = params["weight"];
03534 
03535         if( m_wptr->get_xsize() != m_vnxc+1 && m_wptr->get_ysize() != m_vnyp && m_wptr->get_zsize() != m_vnzp ) {
03536                m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
03537                m_wptr->to_zero();
03538         }
03539 
03540         m_wptr->set_array_offsets(0,1,1);
03541 
03542 }
03543 
03544 int nn4_ctf_rectReconstructor::insert_slice(const EMData* const slice, const Transform& t, const float)
03545 {
03546         // sanity checks
03547         if (!slice) {
03548                 LOGERR("try to insert NULL slice");
03549                 return 1;
03550         }
03551 
03552         int buffed = slice->get_attr_default( "buffed", 0 );
03553         if( buffed > 0 ) {
03554                 int mult = slice->get_attr_default( "mult", 1 );
03555                 insert_buffed_slice( slice, mult );
03556                 return 0;
03557         }
03558 
03559         int padffted= slice->get_attr_default("padffted", 0);
03560         //if( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_vnx)  )
03561         if( padffted==0 && (slice->get_xsize()!=slice->get_ysize())  )
03562         {
03563                 // FIXME: Why doesn't this throw an exception?
03564                 LOGERR("Tried to insert a slice that is the wrong size.");
03565                 return 1;
03566         }
03567 
03568         EMData* padfft = NULL;
03569 
03570         if( padffted != 0 ) padfft = new EMData(*slice);
03571         else                padfft = padfft_slice( slice, t, m_npad );
03572 
03573         int mult= slice->get_attr_default("mult", 1);
03574 
03575         Assert( mult > 0 );
03576         insert_padfft_slice( padfft, t, mult );
03577 
03578         checked_delete( padfft );
03579 
03580         return 0;
03581 }
03582 
03583 int nn4_ctf_rectReconstructor::insert_buffed_slice( const EMData* buffed, int mult )
03584 {
03585         const float* bufdata = buffed->get_data();
03586         float* cdata = m_volume->get_data();
03587         float* wdata = m_wptr->get_data();
03588 
03589         int npoint = buffed->get_xsize()/4;
03590         for( int i=0; i < npoint; ++i ) {
03591 
03592                 int pos2 = int( bufdata[4*i] );
03593                 int pos1 = pos2 * 2;
03594                 cdata[pos1  ] += bufdata[4*i+1]*mult;
03595                 cdata[pos1+1] += bufdata[4*i+2]*mult;
03596                 wdata[pos2  ] += bufdata[4*i+3]*mult;
03597 /*
03598         std::cout << "pos1, pos2, ctfv1, ctfv2, ctf2: ";
03599         std::cout << pos1 << " " << bufdata[5*i+1] << " " << bufdata[5*i+2] << " ";
03600         std::cout << pos2 << " " << bufdata[5*i+4] << std::endl;
03601  */
03602         }
03603         return 0;
03604 }
03605 
03606 
03607 int nn4_ctf_rectReconstructor::insert_padfft_slice( EMData* padfft, const Transform& t, int mult )
03608 {
03609         Assert( padfft != NULL );
03610         float tmp = padfft->get_attr("ctf_applied");
03611         int   ctf_applied = (int) tmp;
03612         vector<Transform> tsym = t.get_sym_proj(m_symmetry);
03613         for (unsigned int isym=0; isym < tsym.size(); isym++) {
03614                 if(ctf_applied) m_volume->insert_rect_slice_ctf_applied(m_wptr, padfft, tsym[isym], m_sizeofprojection, m_xratio,m_yratio, m_npad, mult);
03615                                 
03616                 else            m_volume->insert_rect_slice_ctf(m_wptr, padfft, tsym[isym], m_sizeofprojection, m_xratio, m_yratio, m_npad, mult);
03617                 
03618         }
03619 
03620 
03621         return 0;
03622 
03623 }
03624 
03625 #define  tw(i,j,k)      tw[ i-1 + (j-1+(k-1)*iy)*ix ]
03626 EMData* nn4_ctf_rectReconstructor::finish(bool)
03627 {
03628         m_volume->set_array_offsets(0, 1, 1);
03629         m_wptr->set_array_offsets(0, 1, 1);
03630         m_volume->symplane0_rect(m_wptr);
03631 
03632         int box = 7;
03633         int vol = box*box*box;
03634         int kc = (box-1)/2;
03635         vector< float > pow_a( 3*kc+1, 1.0 );
03636         for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
03637         pow_a[3*kc]=0.0;
03638 
03639 
03640         float max = max3d( kc, pow_a );
03641         float alpha = ( 1.0f - 1.0f/(float)vol ) / max;
03642         float osnr = 1.0f/m_snr;
03643 
03644         // normalize
03645         int ix,iy,iz;
03646         for (iz = 1; iz <= m_vnzp; iz++) {
03647                 for (iy = 1; iy <= m_vnyp; iy++) {
03648                         for (ix = 0; ix <= m_vnxc; ix++) {
03649                                 if ( (*m_wptr)(ix,iy,iz) > 0.0f) {//(*v) should be treated as complex!!
03650                     int iyp = (iy<=m_vnyc) ? iy - 1 : iy-m_vnyp-1;
03651                     int izp = (iz<=m_vnzc) ? iz - 1 : iz-m_vnzp-1;
03652                     float tmp=0.0;
03653                     if( m_varsnr )
03654                     {
03655                         float freq = sqrt( (float)(ix*ix/(m_xratio*m_xratio)+iyp*iyp/(m_zratio*m_yratio)+izp*izp) );
03656                         tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+freq*osnr)*m_sign;
03657                     }
03658                     else
03659                     {
03660                         tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+osnr)*m_sign;
03661                     }
03662 
03663                                         if( m_weighting == ESTIMATE ) {
03664                                                 int cx = ix;
03665                                                 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
03666                                                 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
03667                                                 float sum = 0.0;
03668                                                 for( int ii = -kc; ii <= kc; ++ii ) {
03669                                                         int nbrcx = cx + ii;
03670                                                         if( nbrcx >= m_vnxc ) continue;
03671                                                         for( int jj= -kc; jj <= kc; ++jj ) {
03672                                                                 int nbrcy = cy + jj;
03673                                                                 if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
03674                                                                 for( int kk = -kc; kk <= kc; ++kk ) {
03675                                                                         int nbrcz = cz + jj;
03676                                                                         if( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
03677                                                                         if( nbrcx < 0 ) {
03678                                                                                 nbrcx = -nbrcx;
03679                                                                                 nbrcy = -nbrcy;
03680                                                                                 nbrcz = -nbrcz;
03681                                                                         }
03682 
03683                                                                         int nbrix = nbrcx;
03684                                                                         int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
03685                                                                         int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
03686                                                                         if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0.0 ) {
03687                                                                                 int c = 3*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
03688                                                                                 sum = sum + pow_a[c];
03689                                                                                   // if(ix%20==0 && iy%20==0 && iz%20==0)
03690                                                                                  //   std::cout << boost::format( "%4d %4d %4d %4d %10.3f" ) % nbrix % nbriy % nbriz % c % sum << std::endl;
03691                                                                         }
03692                                                                 }
03693                                                         }
03694                                                 }
03695                                                 float wght = 1.0f / ( 1.0f - alpha * sum );
03696 /*
03697                         if(ix%10==0 && iy%10==0)
03698                         {
03699                             std::cout << boost::format( "%4d %4d %4d " ) % ix % iy %iz;
03700                             std::cout << boost::format( "%10.3f %10.3f %10.3f " )  % tmp % wght % sum;
03701                             std::  << boost::format( "%10.3f %10.3e " ) % pow_b[r] % alpha;
03702                             std::cout << std::endl;
03703                         }
03704  */
03705                                                 tmp = tmp * wght;
03706                                         }
03707                                         (*m_volume)(2*ix,iy,iz) *= tmp;
03708                                         (*m_volume)(2*ix+1,iy,iz) *= tmp;
03709                                 }
03710                         }
03711                 }
03712         }
03713 
03714         // back fft
03715         m_volume->do_ift_inplace();
03716         int npad = m_volume->get_attr("npad");
03717         m_volume->depad();
03718         circumf_rect( m_volume, npad );
03719         m_volume->set_array_offsets( 0, 0, 0 );
03720         return 0;
03721 }
03722 #undef  tw
03723 
03724 
03725 
03726 // Added By Zhengfan Yang on 04/11/07
03727 // Beginning of the addition
03728 // --------------------------------------------------------------------------------
03729 
03730 nnSSNR_ctfReconstructor::nnSSNR_ctfReconstructor()
03731 {
03732         m_volume  = NULL;
03733         m_wptr    = NULL;
03734         m_wptr2   = NULL;
03735         m_wptr3   = NULL;
03736 }
03737 
03738 nnSSNR_ctfReconstructor::nnSSNR_ctfReconstructor( const string& symmetry, int size, int npad, float snr, int sign)
03739 {
03740         m_volume  = NULL;
03741         m_wptr    = NULL;
03742         m_wptr2   = NULL;
03743         m_wptr3   = NULL;
03744 
03745         setup( symmetry, size, npad, snr, sign );
03746 }
03747 
03748 nnSSNR_ctfReconstructor::~nnSSNR_ctfReconstructor()
03749 {
03750 
03751         //if( m_delete_volume )  checked_delete(m_volume);
03752         //if( m_delete_weight )  checked_delete( m_wptr );
03753         //if( m_delete_weight2 ) checked_delete( m_wptr2 );
03754         //if( m_delete_weight3 ) checked_delete( m_wptr3 );
03755         //checked_delete( m_result );
03756 }
03757 
03758 void nnSSNR_ctfReconstructor::setup()
03759 {
03760         int  size = params["size"];
03761         int  npad = params["npad"];
03762         int  sign = params["sign"];
03763         float snr = params["snr"];
03764         string symmetry;
03765         if( params.has_key("symmetry") )  symmetry = params["symmetry"].to_str();
03766         else                              symmetry = "c1";
03767         
03768         setup( symmetry, size, npad, snr, sign );
03769 }
03770 void nnSSNR_ctfReconstructor::setup( const string& symmetry, int size, int npad, float snr, int sign )
03771 {
03772 
03773         m_weighting = ESTIMATE;
03774         m_wghta     = 0.2f;
03775         m_wghtb     = 0.004f;
03776         wiener      = 1;
03777 
03778         m_symmetry  = symmetry;
03779         m_npad      = npad;
03780         m_nsym      = Transform::get_nsym(m_symmetry);
03781 
03782         m_sign      = sign;
03783         m_snr       = snr;
03784 
03785         m_vnx       = size;
03786         m_vny       = size;
03787         m_vnz       = size;
03788 
03789         m_vnxp      = size*npad;
03790         m_vnyp      = size*npad;
03791         m_vnzp      = size*npad;
03792 
03793         m_vnxc      = m_vnxp/2;
03794         m_vnyc      = m_vnyp/2;
03795         m_vnzc      = m_vnzp/2;
03796 
03797         buildFFTVolume();
03798         buildNormVolume();
03799         buildNorm2Volume();
03800         buildNorm3Volume();
03801 }
03802 
03803 void nnSSNR_ctfReconstructor::buildFFTVolume() {
03804 
03805         int offset = 2 - m_vnxp%2;
03806         m_volume = params["fftvol"];
03807 
03808         m_volume->set_size(m_vnxp+offset,m_vnyp,m_vnzp);
03809         m_volume->to_zero();
03810 
03811         m_volume->set_fftodd(m_vnxp % 2);
03812 
03813         m_volume->set_nxc(m_vnxc);
03814         m_volume->set_complex(true);
03815         m_volume->set_ri(true); //(real, imaginary) instead of polar coordinate
03816         m_volume->set_fftpad(true);
03817         m_volume->set_attr("npad", m_npad);
03818         m_volume->set_array_offsets(0,1,1);
03819 }
03820 
03821 
03822 void nnSSNR_ctfReconstructor::buildNormVolume()
03823 {
03824         m_wptr = params["weight"];
03825         m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
03826         m_wptr->to_zero();
03827         m_wptr->set_array_offsets(0,1,1);
03828 }
03829 
03830 void nnSSNR_ctfReconstructor::buildNorm2Volume() {
03831 
03832         m_wptr2 = params["weight2"];
03833         m_wptr2->set_size(m_vnxc+1,m_vnyp,m_vnzp);
03834         m_wptr2->to_zero();
03835         m_wptr2->set_array_offsets(0,1,1);
03836 }
03837 
03838 void nnSSNR_ctfReconstructor::buildNorm3Volume() {
03839 
03840         m_wptr3 = params["weight3"];
03841         m_wptr3->set_size(m_vnxc+1,m_vnyp,m_vnzp);
03842         m_wptr3->to_zero();
03843         m_wptr3->set_array_offsets(0,1,1);
03844 }
03845 
03846 int nnSSNR_ctfReconstructor::insert_slice(const EMData *const  slice, const Transform& t, const float) {
03847         // sanity checks
03848         if (!slice) {
03849                 LOGERR("try to insert NULL slice");
03850                 return 1;
03851         }
03852         int padffted= slice->get_attr_default("padffted", 0);
03853         if ( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_vnx)  ) {
03854                 // FIXME: Why doesn't this throw an exception?
03855                 LOGERR("Tried to insert a slice that is the wrong size.");
03856                 return 1;
03857         }
03858         EMData* padfft = NULL;
03859 
03860         if( padffted != 0 ) padfft = new EMData(*slice);
03861         else                 padfft = padfft_slice( slice, t, m_npad );
03862 
03863         int mult= slice->get_attr_default("mult", 1);
03864 
03865         Assert( mult > 0 );
03866         insert_padfft_slice( padfft, t, mult );
03867 
03868         checked_delete( padfft );
03869         return 0;
03870 }
03871 int nnSSNR_ctfReconstructor::insert_padfft_slice( EMData* padfft, const Transform& t, int mult )
03872 {
03873         Assert( padfft != NULL );
03874 
03875         // insert slice for all symmetry related positions
03876         for (int isym=0; isym < m_nsym; isym++) {
03877                 Transform tsym = t.get_sym(m_symmetry, isym);
03878                 m_volume->nn_SSNR_ctf(m_wptr, m_wptr2, m_wptr3, padfft, tsym, mult);
03879         }
03880 
03881         return 0;
03882 }
03883 
03884 #define  tw(i,j,k)      tw[ i-1 + (j-1+(k-1)*iy)*ix ]
03885 EMData* nnSSNR_ctfReconstructor::finish(bool)
03886 {
03887 /*
03888   I changed the code on 05/15 so it only returns variance.
03889   Lines commented out are marked by //#
03890   The version prior to the currect changes is r1.190
03891   PAP
03892 */
03893         /***
03894             m_volume ctf*(P^2D->3D(F^3D))
03895             m_wptr   ctf^2
03896             m_wptr2  |P^2D->3D(F^3D)|^2
03897             m_wptr3  Kn
03898             nominator = sum_rot [ wght*signal ]
03899             denominator  = sum_rot[ wght*variance ]
03900                                                       ***/
03901         int kz, ky;
03902         int box = 7;
03903         int kc  = (box-1)/2;
03904         float alpha = 0.0;
03905         float argx, argy, argz;
03906         vector< float > pow_a( 3*kc+1, 1.0 );
03907         float w = params["w"];
03908         float dx2 = 1.0f/float(m_vnxc)/float(m_vnxc);
03909         float dy2 = 1.0f/float(m_vnyc)/float(m_vnyc);
03910 #ifdef  _WIN32
03911         float dz2 = 1.0f/_cpp_max(float(m_vnzc),1.0f)/_cpp_max(float(m_vnzc),1.0f);
03912         int inc = Util::round(float(_cpp_max(_cpp_max(m_vnxc,m_vnyc),m_vnzc))/w);
03913 #else
03914         float dz2 = 1.0f/std::max(float(m_vnzc),1.0f)/std::max(float(m_vnzc),1.0f);
03915         int inc = Util::round(float(std::max(std::max(m_vnxc,m_vnyc),m_vnzc))/w);
03916 #endif  //_WIN32
03917 
03918         EMData* SSNR = params["SSNR"];
03919         SSNR->set_size(inc+1,4,1);
03920         //#EMData* vol_ssnr = new EMData();
03921         //#vol_ssnr->set_size(m_vnxp, m_vnyp, m_vnzp);
03922         //#vol_ssnr->to_zero();
03923         //#  new linea follow
03924         EMData* vol_ssnr = params["vol_ssnr"];
03925         vol_ssnr->set_size(m_vnxp+ 2 - m_vnxp%2, m_vnyp ,m_vnzp);
03926         vol_ssnr->to_zero();
03927         if ( m_vnxp % 2 == 0 ) vol_ssnr->set_fftodd(0);
03928         else                   vol_ssnr->set_fftodd(1);
03929         vol_ssnr->set_nxc(m_vnxc);
03930         vol_ssnr->set_complex(true);
03931         vol_ssnr->set_ri(true);
03932         vol_ssnr->set_fftpad(false);
03933         //#
03934         float *nom    = new float[inc+1];
03935         float *denom  = new float[inc+1];
03936         int  *ka     = new int[inc+1];
03937         int  *nn     = new int[inc+1];
03938         float wght = 1.f;
03939         for (int i = 0; i <= inc; i++) {
03940                 nom[i]   = 0.0f;
03941                 denom[i] = 0.0f;
03942                 nn[i]    = 0;
03943                 ka[i]    = 0;
03944         }
03945         m_volume->symplane2(m_wptr, m_wptr2, m_wptr3);
03946         if ( m_weighting == ESTIMATE ) {
03947                 int vol = box*box*box;
03948                 for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
03949                 pow_a[3*kc] = 0.0;
03950                 float max = max3d( kc, pow_a );
03951                 alpha = ( 1.0f - 1.0f/(float)vol ) / max;
03952         }
03953         for (int iz = 1; iz <= m_vnzp; iz++) {
03954                 if ( iz-1 > m_vnzc ) kz = iz-1-m_vnzp; else kz = iz-1;
03955                 argz = float(kz*kz)*dz2;
03956                 for (int iy = 1; iy <= m_vnyp; iy++) {
03957                         if ( iy-1 > m_vnyc ) ky = iy-1-m_vnyp; else ky = iy-1;
03958                         argy = argz + float(ky*ky)*dy2;
03959                         for (int ix = 0; ix <= m_vnxc; ix++) {
03960                                 float Kn = (*m_wptr3)(ix,iy,iz);
03961                                 argx = std::sqrt(argy + float(ix*ix)*dx2);
03962                                 int r = Util::round(float(inc)*argx);
03963                                 if ( r >= 0 && Kn > 4.5f ) {
03964                                         if ( m_weighting == ESTIMATE ) {
03965                                                 int cx = ix;
03966                                                 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
03967                                                 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
03968                                                 float sum = 0.0;
03969                                                 for( int ii = -kc; ii <= kc; ++ii ) {
03970                                                         int nbrcx = cx + ii;
03971                                                         if( nbrcx >= m_vnxc ) continue;
03972                                                         for ( int jj= -kc; jj <= kc; ++jj ) {
03973                                                                 int nbrcy = cy + jj;
03974                                                                 if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
03975                                                                 for( int kk = -kc; kk <= kc; ++kk ) {
03976                                                                         int nbrcz = cz + jj;
03977                                                                         if ( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
03978                                                                         if( nbrcx < 0 ) {
03979                                                                                 nbrcx = -nbrcx;
03980                                                                                 nbrcy = -nbrcy;
03981                                                                                 nbrcz = -nbrcz;
03982                                                                         }
03983                                                                         int nbrix = nbrcx;
03984                                                                         int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
03985                                                                         int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
03986                                                                         if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0 ) {
03987                                                                                 int c = 3*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
03988                                                                                 sum = sum + pow_a[c];
03989                                                                         }
03990                                                                 }
03991                                                         }
03992                                                 }
03993 //                                              int r = std::abs(cx) + std::abs(cy) + std::abs(cz);
03994                                                 wght = 1.0f / ( 1.0f - alpha * sum );
03995                                         } // end of ( m_weighting == ESTIMATE )
03996                                         float nominator   = std::norm(m_volume->cmplx(ix,iy,iz))/(*m_wptr)(ix,iy,iz);
03997                                         float denominator = ((*m_wptr2)(ix,iy,iz)-std::norm(m_volume->cmplx(ix,iy,iz))/(*m_wptr)(ix,iy,iz))/(Kn-1.0f);
03998                                         // Skip Friedel related values
03999                                         if( (ix>0 || (kz>=0 && (ky>=0 || kz!=0)))) {
04000                                                 if ( r <= inc ) {
04001                                                         nom[r]   += nominator*wght;
04002                                                         denom[r] += denominator/(*m_wptr)(ix,iy,iz)*wght;
04003                                                         nn[r]    += 2;
04004                                                         ka[r]    += int(Kn);
04005                                                 }
04006 /*
04007 #ifdef  _WIN32
04008                                                 float  tmp = _cpp_max(nominator/denominator/(*m_wptr)(ix,iy,iz)-1.0f,0.0f);
04009 #else
04010                                                 float  tmp = std::max(nominator/denominator/(*m_wptr)(ix,iy,iz)-1.0f,0.0f);
04011 #endif  //_WIN32
04012                                                 //  Create SSNR as a 3D array (-n/2:n/2+n%2-1)
04013                                                 int iix = m_vnxc + ix; int iiy = m_vnyc + ky; int iiz = m_vnzc + kz;
04014                                                 if( iix >= 0 && iix < m_vnxp && iiy >= 0 && iiy < m_vnyp && iiz >= 0 && iiz < m_vnzp )
04015                                                         (*vol_ssnr)(iix, iiy, iiz) = tmp;
04016                                                 // Friedel part
04017                                                 iix = m_vnxc - ix; iiy = m_vnyc - ky; iiz = m_vnzc - kz;
04018                                                 if( iix >= 0 && iix < m_vnxp && iiy >= 0 && iiy < m_vnyp && iiz >= 0 && iiz < m_vnzp )
04019                                                         (*vol_ssnr)(iix, iiy, iiz) = tmp;
04020 */
04021                                         }
04022                                         (*vol_ssnr)(2*ix, iy-1, iz-1) = denominator*wght;
04023                                 } // end of Kn>4.5 or whatever
04024                         }
04025                 }
04026         }
04027         for (int i = 0; i <= inc; i++) {
04028                 (*SSNR)(i,0,0) = nom[i];
04029                 (*SSNR)(i,1,0) = denom[i];
04030                 (*SSNR)(i,2,0) = static_cast<float>(nn[i]);
04031                 (*SSNR)(i,3,0) = static_cast<float>(ka[i]);
04032         }
04033         vol_ssnr->update();
04034 
04035         delete[] nom;
04036         delete[] denom;
04037         delete[] nn;
04038         delete[] ka;
04039 
04040         return 0;
04041 }
04042 #undef  tw
04043 // -----------------------------------------------------------------------------------
04044 // End of this addition
04045 
04046 
04047 void EMAN::dump_reconstructors()
04048 {
04049         dump_factory < Reconstructor > ();
04050 }
04051 
04052 map<string, vector<string> > EMAN::dump_reconstructors_list()
04053 {
04054         return dump_factory_list < Reconstructor > ();
04055 }
04056 
04057 
04058 
04059 
04060 
04061 using std::ofstream;
04062 using std::ifstream;
04063 
04064 
04065 newfile_store::newfile_store( const string& filename, int npad, bool ctf )
04066     : m_bin_file( filename + ".bin" ),
04067       m_txt_file( filename + ".txt" )
04068 {
04069         m_npad = npad;
04070         m_ctf = ctf;
04071 }
04072 
04073 newfile_store::~newfile_store( )
04074 {
04075 }
04076 
04077 void newfile_store::add_image( EMData* emdata, const Transform& tf )
04078 {
04079         if( m_bin_of == NULL ) {
04080             m_bin_of = shared_ptr<ofstream>( new ofstream(m_bin_file.c_str(), std::ios::out|std::ios::binary) );
04081             m_txt_of = shared_ptr<ofstream>( new ofstream(m_txt_file.c_str()) );
04082         }
04083 
04084 
04085         EMData* padfft = padfft_slice( emdata, tf, m_npad );
04086 
04087         int nx = padfft->get_xsize();
04088         int ny = padfft->get_ysize();
04089         int n2 = ny / 2;
04090         int n = ny;
04091 
04092         float voltage=0.0f, pixel=0.0f, Cs=0.0f, ampcont=0.0f, bfactor=0.0f, defocus=0.0f;
04093 
04094         if( m_ctf ) {
04095                 Ctf* ctf = emdata->get_attr( "ctf" );
04096                 Dict params = ctf->to_dict();
04097                 voltage = params["voltage"];
04098                 pixel   = params["apix"];
04099                 Cs      = params["cs"];
04100                 ampcont = params["ampcont"];
04101                 bfactor = params["bfactor"];
04102                 defocus = params["defocus"];
04103                 if(ctf) {delete ctf; ctf=0;}
04104         }
04105 
04106         vector<point_t> points;
04107         for( int j=-ny/2+1; j <= ny/2; j++ ) {
04108                 int jp = (j>=0) ? j+1 : ny+j+1;
04109                 for( int i=0; i <= n2; ++i ) {
04110                         int r2 = i*i + j*j;
04111                         if( (r2<ny*ny/4) && !( (i==0) && (j<0) ) ) {
04112                                 float ctf;
04113                                 if( m_ctf ) {
04114                                         float ak = std::sqrt( r2/float(ny*ny) )/pixel;
04115                                         ctf = Util::tf( defocus, ak, voltage, Cs, ampcont, bfactor, 1);
04116                                 } else {
04117                                         ctf = 1.0;
04118                                 }
04119 
04120                                 float xnew = i*tf[0][0] + j*tf[1][0];
04121                                 float ynew = i*tf[0][1] + j*tf[1][1];
04122                                 float znew = i*tf[0][2] + j*tf[1][2];
04123                                 std::complex<float> btq;
04124                                 if (xnew < 0.) {
04125                                         xnew = -xnew;
04126                                         ynew = -ynew;
04127                                         znew = -znew;
04128                                         btq = conj(padfft->cmplx(i,jp-1));
04129                                 } else {
04130                                         btq = padfft->cmplx(i,jp-1);
04131                                 }
04132 
04133                                 int ixn = int(xnew + 0.5 + n) - n;
04134                                 int iyn = int(ynew + 0.5 + n) - n;
04135                                 int izn = int(znew + 0.5 + n) - n;
04136                                 if ((ixn <= n2) && (iyn >= -n2) && (iyn <= n2) && (izn >= -n2) && (izn <= n2)) {
04137                                         int ixf, iyf, izf;
04138                                         if (ixn >= 0) {
04139                                                 int iza, iya;
04140                                                 if (izn >= 0)
04141                                                     iza = izn + 1;
04142                                                 else
04143                                                     iza = n + izn + 1;
04144 
04145                                                 if (iyn >= 0)
04146                                                     iya = iyn + 1;
04147                                                 else
04148                                                     iya = n + iyn + 1;
04149 
04150                                                 ixf = ixn;
04151                                                 iyf = iya;
04152                                                 izf = iza;
04153                                         } else {
04154                                                 int izt, iyt;
04155                                                 if (izn > 0)
04156                                                     izt = n - izn + 1;
04157                                                 else
04158                                                     izt = -izn + 1;
04159 
04160                                                 if (iyn > 0)
04161                                                     iyt = n - iyn + 1;
04162                                                 else
04163                                                     iyt = -iyn + 1;
04164 
04165                                                 ixf = -ixn;
04166                                                 iyf = iyt;
04167                                                 izf = izt;
04168                                         }
04169 
04170 
04171                                         int pos2 = ixf + (iyf-1)*nx/2 + (izf-1)*ny*nx/2;
04172                                         float ctfv1 = btq.real() * ctf;
04173                                         float ctfv2 = btq.imag() * ctf;
04174                                         float ctf2 = ctf*ctf;
04175 
04176                                         point_t p;
04177                                         p.pos2 = pos2;
04178                                         p.real = ctfv1;
04179                                         p.imag = ctfv2;
04180                                         p.ctf2 = ctf2;
04181 
04182                                         points.push_back( p );
04183                                 }
04184                         }
04185                 }
04186         }
04187 
04188 
04189         int npoint = points.size();
04190         std::istream::off_type offset = (m_offsets.size()==0) ? 0 : m_offsets.back();
04191         offset += npoint*sizeof(point_t);
04192         m_offsets.push_back( offset );
04193 
04194         *m_txt_of << m_offsets.back() << std::endl;
04195         m_bin_of->write( (char*)(&points[0]), sizeof(point_t)*npoint );
04196         checked_delete( padfft );
04197 }
04198 
04199 void newfile_store::get_image( int id, EMData* buf )
04200 {
04201         if( m_offsets.size()==0 ) {
04202                 ifstream is( m_txt_file.c_str() );
04203                 std::istream::off_type off;
04204                 while( is >> off ) {
04205                     m_offsets.push_back( off );
04206                 }
04207 
04208                 m_bin_if = shared_ptr<std::ifstream>( new ifstream(m_bin_file.c_str(), std::ios::in|std::ios::binary) );
04209         }
04210 
04211         Assert( m_bin_if != NULL );
04212 
04213         std::istream::off_type offset = (id==0) ? 0 : m_offsets[id-1];
04214         Assert( offset >= 0 );
04215         m_bin_if->seekg( offset, std::ios::beg );
04216 
04217 
04218         if( m_bin_if->bad() || m_bin_if->fail() || m_bin_if->eof() ) {
04219                 std::cout << "bad or fail or eof while fetching id, offset: " << id << " " << offset << std::endl;
04220                 throw std::logic_error( "bad happen" );
04221         }
04222 
04223         int bufsize = (m_offsets[id] - offset)/sizeof(float);
04224         if( buf->get_xsize() != bufsize ) {
04225                 buf->set_size( bufsize, 1, 1 );
04226         }
04227 
04228         char* data = (char*)(buf->get_data());
04229         m_bin_if->read( data, sizeof(float)*bufsize );
04230         buf->update();
04231 }
04232 
04233 void newfile_store::read( int nprj )
04234 {
04235         if( m_offsets.size()==0 ) {
04236             ifstream is( m_txt_file.c_str() );
04237             std::istream::off_type off;
04238             while( is >> off ) {
04239                 m_offsets.push_back( off );
04240             }
04241         }
04242 
04243         if( m_bin_if==NULL ) {
04244             m_bin_if = shared_ptr< ifstream>( new ifstream(m_bin_file.c_str(), std::ios::in|std::ios::binary) );
04245         }
04246 
04247 
04248         int npoint = m_offsets[0]/sizeof(point_t);
04249         std::ios::off_type prjsize = m_offsets[0];
04250 
04251         try {
04252             m_points.resize(nprj * npoint);
04253         }
04254         catch( std::exception& e ) {
04255             std::cout << "Error: " << e.what() << std::endl;
04256         }
04257 
04258         int ip = 0;
04259         for( int i=0; i < nprj; ++i ) {
04260                 m_bin_if->read( (char*)(&m_points[ip]), prjsize );
04261                 if( m_bin_if->bad() || m_bin_if->fail() || m_bin_if->eof() )  {
04262                         std::cout << "Error: file hander bad or fail or eof" << std::endl;
04263                         return;
04264                 }
04265                 ip += npoint;
04266         }
04267 }
04268 
04269 void newfile_store::add_tovol( EMData* fftvol, EMData* wgtvol, const vector<int>& mults, int pbegin, int pend )
04270 {
04271         float* vdata = fftvol->get_data();
04272         float* wdata = wgtvol->get_data();
04273 
04274         int npoint = m_offsets[0]/sizeof(point_t);
04275 //    Assert( int(mults.size())==nprj );
04276         Assert( int(m_points.size())== (pend - pbegin)*npoint );
04277 
04278         for( int iprj=pbegin; iprj < pend; ++iprj ) {
04279                 int m = mults[iprj];
04280                 if( m==0 ) continue;
04281 
04282                 int ipt = (iprj-pbegin)*npoint;
04283                 for( int i=0; i < npoint; ++i )  {
04284                     int pos2 = m_points[ipt].pos2;
04285                     int pos1 = pos2*2;
04286 
04287                     wdata[pos2]   += m_points[ipt].ctf2*m;
04288                     vdata[pos1]   += m_points[ipt].real*m;
04289                     vdata[pos1+1] += m_points[ipt].imag*m;
04290                     ++ipt;
04291                 }
04292         }
04293 }
04294 
04295 void newfile_store::restart()
04296 {
04297     m_bin_if = shared_ptr< ifstream>( new ifstream(m_bin_file.c_str(), std::ios::in|std::ios::binary) );
04298 }
04299 
04300 file_store::file_store(const string& filename, int npad, int write, bool ctf)
04301     : m_bin_file(filename + ".bin"),
04302       m_txt_file(filename + ".txt")
04303 {
04304     m_ctf = ctf;
04305     m_prev = -1;
04306     m_npad = npad;
04307     m_write = write;
04308 }
04309 
04310 file_store::~file_store()
04311 {
04312 }
04313 
04314 void file_store::add_image( EMData* emdata, const Transform& tf )
04315 {
04316 
04317     EMData* padfft = padfft_slice( emdata, tf, m_npad );
04318 
04319     float* data = padfft->get_data();
04320 
04321     if( m_write && m_bin_ohandle == NULL )
04322     {
04323         m_bin_ohandle = shared_ptr< ofstream >( new ofstream(m_bin_file.c_str(), std::ios::out | std::ios::binary) );
04324         m_txt_ohandle = shared_ptr< ofstream >( new ofstream(m_txt_file.c_str() ) );
04325         if( m_ctf )
04326                  *m_txt_ohandle << "Cs pixel voltage ctf_applied amp_contrast defocus ";
04327 
04328         *m_txt_ohandle << "phi theta psi" << std::endl;
04329     }
04330 
04331     m_x_out = padfft->get_xsize();
04332     m_y_out = padfft->get_ysize();
04333     m_z_out = padfft->get_zsize();
04334     m_totsize = m_x_out*m_y_out*m_z_out;
04335 
04336     if( m_ctf )
04337     {
04338         Ctf* ctf = padfft->get_attr( "ctf" );
04339         Dict ctf_params = ctf->to_dict();
04340 
04341         m_ctf_applied = padfft->get_attr( "ctf_applied" );
04342 
04343         m_Cs = ctf_params["cs"];
04344         m_pixel = ctf_params["apix"];
04345         m_voltage = ctf_params["voltage"];
04346         m_amp_contrast = ctf_params["ampcont"];
04347         m_defocuses.push_back( ctf_params["defocus"] );
04348         if(ctf) {delete ctf; ctf=0;}
04349     }
04350 
04351     Dict params = tf.get_rotation( "spider" );
04352     float phi = params.get( "phi" );
04353     float tht = params.get( "theta" );
04354     float psi = params.get( "psi" );
04355 
04356 
04357     m_phis.push_back( phi );
04358     m_thetas.push_back( tht );
04359     m_psis.push_back( psi );
04360 
04361     if( m_write )
04362     {
04363         m_bin_ohandle->write( (char*)data, sizeof(float)*m_totsize );
04364 
04365         if( m_ctf )
04366         {
04367             *m_txt_ohandle << m_Cs << " ";
04368             *m_txt_ohandle << m_pixel << " ";
04369             *m_txt_ohandle << m_voltage << " ";
04370             *m_txt_ohandle << m_ctf_applied << " ";
04371             *m_txt_ohandle << m_amp_contrast << " ";
04372             *m_txt_ohandle << m_defocuses.back() << " ";
04373         }
04374         *m_txt_ohandle << m_phis.back() << " ";
04375         *m_txt_ohandle << m_thetas.back() << " ";
04376         *m_txt_ohandle << m_psis.back() << " ";
04377         *m_txt_ohandle << m_x_out << " ";
04378         *m_txt_ohandle << m_y_out << " ";
04379         *m_txt_ohandle << m_z_out << " ";
04380         *m_txt_ohandle << m_totsize << std::endl;
04381     }
04382 
04383     checked_delete(padfft);
04384 
04385 }
04386 
04387 void file_store::get_image( int id, EMData* padfft )
04388 {
04389 
04390     if( m_phis.size() == 0 ) {
04391         ifstream m_txt_ifs( m_txt_file.c_str() );
04392 
04393         if( !m_txt_ifs )
04394         {
04395             std::cerr << "Error: file " << m_txt_file << " does not exist" << std::endl;
04396         }
04397 
04398         string line;
04399         std::getline( m_txt_ifs, line );
04400 
04401         float first, defocus, phi, theta, psi;
04402 
04403 
04404 
04405         while( m_txt_ifs >> first ) {
04406 
04407             if( m_ctf )
04408             {
04409                 m_Cs = first;
04410                 m_txt_ifs >> m_pixel >> m_voltage;
04411                 m_txt_ifs >> m_ctf_applied >> m_amp_contrast;
04412                 m_txt_ifs >> defocus >> phi >> theta >> psi;
04413                 m_defocuses.push_back( defocus );
04414             }
04415             else
04416             {
04417                 phi = first;
04418                 m_txt_ifs >> theta >> psi;
04419             }
04420 
04421             m_txt_ifs >> m_x_out >> m_y_out >> m_z_out >> m_totsize;
04422             m_phis.push_back( phi );
04423             m_thetas.push_back( theta );
04424             m_psis.push_back( psi );
04425         }
04426     }
04427 
04428     Assert( m_ihandle != NULL );
04429 
04430     std::istream::off_type offset = id*sizeof(float)*m_totsize;
04431     Assert( offset >= 0 );
04432 
04433     if( offset > 0 )
04434     {
04435         m_ihandle->seekg(offset, std::ios::beg);
04436     }
04437 
04438     if( m_ihandle->bad() )
04439     {
04440         std::cout << "bad while fetching id, offset: " << id << " " << offset << std::endl;
04441         throw std::logic_error( "bad happen" );
04442     }
04443 
04444     if( m_ihandle->fail() )
04445     {
04446         std::cout << "fail while fetching id, offset, curoff: " << id << " " << offset << std::endl;
04447         throw std::logic_error( "fail happen" );
04448     }
04449 
04450     if( m_ihandle->eof() )
04451     {
04452         std::cout << "eof while fetching id, offset: " << id << " " << offset << std::endl;
04453         throw std::logic_error( "eof happen" );
04454     }
04455 
04456     if( padfft->get_xsize() != m_x_out ||
04457         padfft->get_ysize() != m_y_out ||
04458         padfft->get_zsize() != m_z_out )
04459     {
04460         padfft->set_size(m_x_out, m_y_out, m_z_out);
04461     }
04462 
04463     char* data = (char*)(padfft->get_data());
04464     m_ihandle->read( data, sizeof(float)*m_totsize );
04465     padfft->update();
04466 
04467     if( m_ctf )
04468     {
04469         padfft->set_attr( "Cs", m_Cs );
04470         padfft->set_attr( "Pixel_size", m_pixel );
04471         padfft->set_attr( "voltage", m_voltage );
04472         padfft->set_attr( "ctf_applied", m_ctf_applied );
04473         padfft->set_attr( "amp_contrast", m_amp_contrast );
04474         padfft->set_attr( "defocus", m_defocuses[id] );
04475     }
04476 
04477     padfft->set_attr( "padffted", 1 );
04478     padfft->set_attr( "phi", m_phis[id] );
04479     padfft->set_attr( "theta", m_thetas[id] );
04480     padfft->set_attr( "psi", m_psis[id] );
04481 
04482 }
04483 
04484 void file_store::restart( )
04485 {
04486     if( m_ihandle == NULL )
04487     {
04488         m_ihandle = shared_ptr< ifstream >( new ifstream(m_bin_file.c_str(), std::ios::in | std::ios::binary) );
04489     }
04490 
04491     if( m_ihandle->bad() || m_ihandle->fail() || m_ihandle->eof() )
04492     {
04493         m_ihandle->open( m_bin_file.c_str(), std::ios::binary );
04494     }
04495 
04496     m_ihandle->seekg( 0, std::ios::beg );
04497 }
04498 
04499 /* vim: set ts=4 noet: */

Generated on Tue Jul 12 13:45:51 2011 for EMAN2 by  doxygen 1.4.7