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