Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

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

Generated on Tue Jun 11 13:40:43 2013 for EMAN2 by  doxygen 1.3.9.1