EMAN2
reconstructor_tools.cpp
Go to the documentation of this file.
00001 
00005 /*
00006  * Author: David Woolford, 07/25/2007 (woolford@bcm.edu)
00007  * Copyright (c) 2000-2007 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 <cstring>
00037 #include "reconstructor_tools.h"
00038 
00039 using namespace EMAN;
00040 
00041 const string FourierInserter3DMode1::NAME = "nearest_neighbor";
00042 const string FourierInserter3DMode2::NAME = "gauss_2";
00043 const string FourierInserter3DMode3::NAME = "gauss_3";
00044 //const string FourierInserter3DMode4::NAME = "gauss_4";
00045 const string FourierInserter3DMode5::NAME = "gauss_5";
00046 const string FourierInserter3DMode6::NAME = "gauss_5_slow";
00047 const string FourierInserter3DMode7::NAME = "hypergeom_5";
00048 const string FourierInserter3DMode8::NAME = "experimental";
00049 
00050 template <> Factory < FourierPixelInserter3D >::Factory()
00051 {
00052         force_add<FourierInserter3DMode1>();
00053         force_add<FourierInserter3DMode2>();
00054         force_add<FourierInserter3DMode3>();
00055 //      force_add<FourierInserter3DMode4>();
00056         force_add<FourierInserter3DMode5>();
00057         force_add<FourierInserter3DMode6>();
00058         force_add<FourierInserter3DMode7>();
00059 //      force_add(&FourierInserter3DMode8::NEW);
00060 }
00061 
00062 
00063 
00064 void FourierPixelInserter3D::init()
00065 {
00066         if ( params.has_key("data") )
00067         {
00068                 data = params["data"];
00069                 if ( data == 0 )
00070                         throw NotExistingObjectException("data", "error the data pointer was 0 in FourierPixelInserter3D::init");
00071         }
00072         else throw NotExistingObjectException("data", "the data pointer was not defined in FourierPixelInserter3D::init");
00073 
00074         if ( params.has_key("norm"))
00075         {
00076                 norm = params["norm"];
00077                 if ( norm == 0 )
00078                         throw NotExistingObjectException("norm", "error the norm pointer was 0 in FourierPixelInserter3D::init");
00079         }
00080         else throw NotExistingObjectException("norm", "the norm pointer was not defined in FourierPixelInserter3D::init");
00081 
00082         nx=data->get_xsize();
00083         ny=data->get_ysize();
00084         nz=data->get_zsize();
00085         nxyz=(size_t)nx*ny*nz;
00086         nx2=nx/2-1;
00087         ny2=ny/2;
00088         nz2=nz/2;
00089         
00090         if (data->has_attr("subvolume_x0") && data->has_attr("subvolume_full_nx")) {
00091                 subx0=data->get_attr("subvolume_x0");
00092                 suby0=data->get_attr("subvolume_y0");
00093                 subz0=data->get_attr("subvolume_z0");
00094                 fullnx=data->get_attr("subvolume_full_nx");
00095                 fullny=data->get_attr("subvolume_full_ny");
00096                 fullnz=data->get_attr("subvolume_full_nz");
00097         }
00098         else {
00099                 subx0=suby0=subz0=-1;
00100         }
00101 }
00102 
00103 bool FourierInserter3DMode1::insert_pixel(const float& xx, const float& yy, const float& zz, const std::complex<float> dt, const float& weight)
00104 {
00105         int x0 = (int) floor(xx + 0.5f);
00106         int y0 = (int) floor(yy + 0.5f);
00107         int z0 = (int) floor(zz + 0.5f);
00108 
00109         size_t off;
00110         if (subx0<0) off=data->add_complex_at(x0,y0,z0,dt*weight);
00111         else off=data->add_complex_at(x0,y0,z0,subx0,suby0,subz0,fullnx,fullny,fullnz,dt*weight);
00112         if (static_cast<int>(off)!=nxyz) norm[off/2]+=weight;
00113         else return false;
00114         
00115         return true;
00116 }
00117 
00118 bool FourierInserter3DMode2::insert_pixel(const float& xx, const float& yy, const float& zz, const std::complex<float> dt,const float& weight)
00119 {
00120         int x0 = (int) floor(xx);
00121         int y0 = (int) floor(yy);
00122         int z0 = (int) floor(zz);
00123         
00124         if (subx0<0) {                  // normal full reconstruction
00125                 if (x0<-nx2-1 || y0<-ny2-1 || z0<-nz2-1 || x0>nx2 || y0>ny2 || z0>nz2 ) return false;
00126 
00127                 // no error checking on add_complex_fast, so we need to be careful here
00128                 int x1=x0+1;
00129                 int y1=y0+1;
00130                 int z1=z0+1;
00131                 if (x0<-nx2) x0=-nx2;
00132                 if (x1>nx2) x1=nx2;
00133                 if (y0<-ny2) y0=-ny2;
00134                 if (y1>ny2) y1=ny2;
00135                 if (z0<-nz2) z0=-nz2;
00136                 if (z1>nz2) z1=nz2;
00137                 
00138 //              float h=2.0/((1.0+pow(Util::hypot3sq(xx,yy,zz),.5))*EMConsts::I2G);
00139                 float h=1.0f/EMConsts::I2G;
00140                 //size_t idx;
00141                 float r, gg;
00142 //              int pc=0;
00143                 for (int k = z0 ; k <= z1; k++) {
00144                         for (int j = y0 ; j <= y1; j++) {
00145                                 for (int i = x0; i <= x1; i ++) {
00146                                         r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
00147 //                                      gg=weight;
00148                                         gg = Util::fast_exp(-r *h)*weight;
00149 //                                      gg = Util::fast_exp(-r / EMConsts::I2G)*weight;
00150 //                                      gg = sqrt(Util::fast_exp(-r / EMConsts::I2G))*weight;
00151                                         
00152                                         size_t off;
00153                                         off=data->add_complex_at_fast(i,j,k,dt*gg);
00154 //                                      off=data->add_complex_at(i,j,k,dt*gg);
00155                                         norm[off/2]+=gg;
00156                                 }
00157                         }
00158                 }
00159                 return true;
00160         } 
00161         else {                                  // for subvolumes, not optimized yet
00162                 //size_t idx;
00163                 float r, gg;
00164                 int pc=0;
00165                 for (int k = z0 ; k <= z0 + 1; k++) {
00166                         for (int j = y0 ; j <= y0 + 1; j++) {
00167                                 for (int i = x0; i <= x0 + 1; i ++) {
00168                                         r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
00169                                         gg = Util::fast_exp(-r / EMConsts::I2G)*weight;
00170 
00171                                         size_t off;
00172                                         if (subx0<0) off=data->add_complex_at(i,j,k,dt*gg);
00173                                         else off=data->add_complex_at(i,j,k,subx0,suby0,subz0,fullnx,fullny,fullnz,dt*gg);
00174                                         if (static_cast<int>(off)!=nxyz) { norm[off/2]+=gg; pc+=1; }
00175                                 }
00176                         }
00177                 }
00178                 
00179                 if (pc>0)  return true;
00180                 return false;
00181         }
00182 }
00183 
00184 
00185 bool FourierInserter3DMode3::insert_pixel(const float& xx, const float& yy, const float& zz, const std::complex<float> dt,const float& weight)
00186 {
00187         int x0 = (int) floor(xx-.5);
00188         int y0 = (int) floor(yy-.5);
00189         int z0 = (int) floor(zz-.5);
00190         
00191         if (subx0<0) {                  // normal full reconstruction
00192                 if (x0<-nx2-2 || y0<-ny2-2 || z0<-nz2-2 || x0>nx2+1 || y0>ny2+1 || z0>nz2+1 ) return false;
00193 
00194                 // no error checking on add_complex_fast, so we need to be careful here
00195                 int x1=x0+2;
00196                 int y1=y0+2;
00197                 int z1=z0+2;
00198                 if (x0<-nx2) x0=-nx2;
00199                 if (x1>nx2) x1=nx2;
00200                 if (y0<-ny2) y0=-ny2;
00201                 if (y1>ny2) y1=ny2;
00202                 if (z0<-nz2) z0=-nz2;
00203                 if (z1>nz2) z1=nz2;
00204                 
00205 //              float h=2.0/((1.0+pow(Util::hypot3sq(xx,yy,zz),.5))*EMConsts::I2G);
00206 //              float h=2.0/EMConsts::I3G;
00207                 float h=32.0f/((8.0f+Util::hypot3(xx,yy,zz))*EMConsts::I3G);
00208 //              float w=weight;
00209                 float w=weight/(1.0f+6.0f*Util::fast_exp(-h)+12*Util::fast_exp(-h*2.0f)+8*Util::fast_exp(-h*3.0f));     // approx normalization so higer radii aren't upweighted relative to lower due to wider Gaussian
00210                 //size_t idx;
00211                 float r, gg;
00212 //              int pc=0;
00213                 for (int k = z0 ; k <= z1; k++) {
00214                         for (int j = y0 ; j <= y1; j++) {
00215                                 for (int i = x0; i <= x1; i ++) {
00216                                         r = Util::hypot3sq((float) i - xx, (float)j - yy, (float)k - zz);
00217 //                                      gg=weight;
00218                                         gg = Util::fast_exp(-r *h)*w;
00219 //                                      gg = Util::fast_exp(-r / EMConsts::I2G)*weight;
00220 //                                      gg = sqrt(Util::fast_exp(-r / EMConsts::I2G))*weight;
00221                                         
00222                                         size_t off;
00223                                         off=data->add_complex_at_fast(i,j,k,dt*gg);
00224                                         norm[off/2]+=gg;
00225                                 }
00226                         }
00227                 }
00228                 return true;
00229         }
00230         printf("region writing not supported in mode 3\n");
00231         return false;
00232 }
00233 
00234 
00235 bool FourierInserter3DMode5::insert_pixel(const float& xx, const float& yy, const float& zz, const std::complex<float> dt,const float& weight)
00236 {
00237         int x0 = (int) floor(xx-1.5);
00238         int y0 = (int) floor(yy-1.5);
00239         int z0 = (int) floor(zz-1.5);
00240         
00241         if (subx0<0) {                  // normal full reconstruction
00242                 if (x0<-nx2-4 || y0<-ny2-4 || z0<-nz2-4 || x0>nx2+3 || y0>ny2+3 || z0>nz2+3 ) return false;
00243 
00244                 // no error checking on add_complex_fast, so we need to be careful here
00245                 int x1=x0+4;
00246                 int y1=y0+4;
00247                 int z1=z0+4;
00248                 if (x0<-nx2) x0=-nx2;
00249                 if (x1>nx2) x1=nx2;
00250                 if (y0<-ny2) y0=-ny2;
00251                 if (y1>ny2) y1=ny2;
00252                 if (z0<-nz2) z0=-nz2;
00253                 if (z1>nz2) z1=nz2;
00254                 
00255 //              float h=2.0/((1.0+pow(Util::hypot3sq(xx,yy,zz),.5))*EMConsts::I2G);
00256 //              float h=2.0/EMConsts::I3G;
00257                 float h=32.0f/((8.0f+Util::hypot3(xx,yy,zz))*EMConsts::I3G); 
00258                 float w=weight/(1.0f+6.0f*Util::fast_exp(-h)+12*Util::fast_exp(-h*2.0f)+8*Util::fast_exp(-h*3.0f)+
00259                         6.0f*Util::fast_exp(-h*4.0f)+24.0f*Util::fast_exp(-h*5.0f)+24.0f*Util::fast_exp(-h*6.0f)+12.0f*Util::fast_exp(-h*8.0f)+
00260                         24.0f*Util::fast_exp(-h*9.0f)+8.0f*Util::fast_exp(-h*12.0f));   // approx normalization so higer radii aren't upweighted relative to lower due to wider Gaussian
00261                 //size_t idx;
00262                 float r, gg;
00263 //              int pc=0;
00264                 for (int k = z0 ; k <= z1; k++) {
00265                         for (int j = y0 ; j <= y1; j++) {
00266                                 for (int i = x0; i <= x1; i ++) {
00267                                         r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
00268 //                                      gg=weight;
00269                                         gg = Util::fast_exp(-r *h)*w;
00270 //                                      gg = Util::fast_exp(-r / EMConsts::I2G)*weight;
00271 //                                      gg = sqrt(Util::fast_exp(-r / EMConsts::I2G))*weight;
00272                                         
00273                                         size_t off;
00274                                         off=data->add_complex_at_fast(i,j,k,dt*gg);
00275                                         norm[off/2]+=gg;
00276                                 }
00277                         }
00278                 }
00279                 return true;
00280         }
00281         printf("region writing not supported in mode 3\n");
00282         return false;
00283 }
00284 
00285 
00286 bool FourierInserter3DMode6::insert_pixel(const float& xx, const float& yy, const float& zz, const std::complex<float> dt,const float& weight)
00287 {
00288         int x0 = 2 * (int) floor(xx + 0.5f);
00289         int y0 = (int) floor(yy + 0.5f);
00290         int z0 = (int) floor(zz + 0.5f);
00291 
00292         if (x0<-6||y0<-3||z0<-3||x0>nx+6||y0>ny+3||z0>nz+3) return false;
00293 
00294         //size_t idx;
00295         float r, gg;
00296         for (int k = z0 - 1; k <= z0 + 1; k++) {
00297                 for (int j = y0 - 1; j <= y0 + 1; j++) {
00298                         for (int i = x0 -2; i <= x0 + 2; i += 2) {
00299                                 if (k<0 || j<0 || i<0 || k>=nz || j>=ny || i>=nx) continue;
00300                                 r = Util::hypot3sq((float) i / 2 - xx, j - yy, k - zz);
00301                                 gg = exp(-r / EMConsts::I5G)*weight;
00302 
00303                                 size_t off;
00304                                 if (subx0<0) off=data->add_complex_at(i,j,k,dt*gg);
00305                                 else off=data->add_complex_at(i,j,k,subx0,suby0,subz0,fullnx,fullny,fullnz,dt*gg);
00306                                 if (static_cast<int>(off)!=nxyz) norm[off/2]+=gg;
00307 
00308                         }
00309                 }
00310         }
00311         return true;
00312  }
00313 
00314 
00315 bool FourierInserter3DMode7::insert_pixel(const float& xx, const float& yy, const float& zz, const std::complex<float> dt,const float& weight)
00316 {
00317         int x0 = 2 * (int) floor(xx + 0.5f);
00318         int y0 = (int) floor(yy + 0.5f);
00319         int z0 = (int) floor(zz + 0.5f);
00320 
00321         if (x0 >= nx - 4 || y0 > ny - 3 || z0 > nz - 3 || y0 < 2 || z0 < 2)
00322                 return false;
00323 
00324         int l = x0 - 4;
00325         if (x0 <= 2)
00326                 l = 0;
00327 
00328         //size_t ii;
00329         float r, gg;
00330         for (int k = z0 - 2; k <= z0 + 2; k++) {
00331                 for (int j = y0 - 2; j <= y0 + 2; j++) {
00332                         for (int i = l; i <= x0 + 4; i += 2) {
00333                                 r =     sqrt(Util::hypot3sq((float) i / 2 - xx, (float) j - yy, (float) k - zz));
00334                                 gg = Interp::hyperg(r)*weight;
00335 
00336                                 size_t off;
00337                                 if (subx0<0) off=data->add_complex_at(i,j,k,dt*gg);
00338                                 else off=data->add_complex_at(i,j,k,subx0,suby0,subz0,fullnx,fullny,fullnz,dt*gg);
00339                                 if (static_cast<int>(off)!=nxyz) norm[off/2]+=gg;
00340 
00341                         }
00342                 }
00343         }
00344 
00345         if (x0 <= 2) {
00346                 float xx_b = -xx;
00347                 float yy_b = -(yy - ny / 2) + ny / 2;
00348                 float zz_b = -(zz - nz / 2) + nz / 2;
00349                 x0 = 2 * (int) floor(xx_b + 0.5f);
00350                 y0 = (int) floor(yy_b + 0.5f);
00351                 z0 = (int) floor(zz_b + 0.5f);
00352 
00353                 if (y0 > ny - 3 || z0 > nz - 3 || y0 < 2 || z0 < 2)
00354                         return false;
00355 
00356                 //size_t ii;
00357                 float r, gg;
00358                 for (int k = z0 - 2; k <= z0 + 2; k++) {
00359                         for (int j = y0 - 2; j <= y0 + 2; j++) {
00360                                 for (int i = 0; i <= x0 + 4; i += 2) {
00361                                         r = sqrt(Util::hypot3sq((float) i / 2 - xx_b, (float) j - yy_b,
00362                                                                    (float) k - zz_b));
00363                                         gg = Interp::hyperg(r)*weight;
00364 
00365                                         size_t off;
00366                                         if (subx0<0) off=data->add_complex_at(i,j,k,dt*gg);
00367                                         else off=data->add_complex_at(i,j,k,subx0,suby0,subz0,fullnx,fullny,fullnz,dt*gg);
00368                                         if (static_cast<int>(off)!=nxyz) norm[off/2]+=gg;
00369 
00370                                 }
00371                         }
00372                 }
00373         }
00374 
00375         return true;
00376 }
00377 
00378 
00379 void FourierInserter3DMode8::init()
00380 {
00381         FourierPixelInserter3D::init();
00382         int P = (int)((1.0+0.25)*nx+1);
00383         float r = (float)(nx+1)/(float)P;
00384         mFreqCutoff = 2;
00385         mDFreq = 0.2f;
00386         if (W != 0) delete [] W;
00387         W = Util::getBaldwinGridWeights(mFreqCutoff, (float)P, r,mDFreq,0.5f,0.2f);
00388 
00389 }
00390 bool FourierInserter3DMode8::insert_pixel(const float&, const float&, const float&, const std::complex<float>, const float&)
00391 {
00392 //      int x0 = (int) floor(qx);
00393 //      int y0 = (int) floor(qy);
00394 //      int z0 = (int) floor(qz);
00395 
00396 //      int sizeW = (int)(1+2*mFreqCutoff/mDFreq);
00397 //      int sizeWmid = sizeW/2;
00398 
00399 //      for (int z = z0-mFreqCutoff; z < z0+mFreqCutoff; ++z){
00400 //              for (int y = y0-mFreqCutoff; y < y0+mFreqCutoff; ++y){
00401 //                      for (int x = x0-mFreqCutoff; x < x0+mFreqCutoff; ++x){
00402 //                              if ( x < 0 || x >= nx ) continue;
00403 //                              if ( y < 0 || y >= ny ) continue;
00404 //                              if ( z < 0 || z >= nz ) continue;
00405 //                              float dist = (float)((x-x0)*(x-x0)+(y-y0)*(y-y0)+(z-z0)*(z-z0));
00406 //                              dist = sqrtf(dist);
00407 //                              // We enforce a spherical kernel
00408 //                              if ( dist > mFreqCutoff ) continue;
00409 //                              int idx = (int)(sizeWmid + dist/mDFreq);
00410 //                              if (idx >= sizeW) throw;
00411 //                              float residual = dist/mDFreq - (int)(dist/mDFreq);
00412 //                              if ( fabs(residual) > 1) throw;
00413 // 
00414 //                              float factor = W[idx]*(1.0f-residual)+W[idx+1]*residual*weight;
00415 // 
00416 //                              size_t k = z*nx*ny + y*nx + 2*x;
00417 // 
00418 // //                           float c = Util::agauss(1, x-x0,y-y0,z-z0, EMConsts::I2G);
00419 //                              rdata[k] += fq[0]*factor;
00420 //                              rdata[k+1] += fq[1]*factor;
00421 // 
00422 // 
00423 //                              norm[k/2] += weight;
00424 // 
00425 //                      }
00426 //              }
00427 //      }
00428 
00429         return true;
00430 }