00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
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
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
00056 force_add<FourierInserter3DMode5>();
00057 force_add<FourierInserter3DMode6>();
00058 force_add<FourierInserter3DMode7>();
00059
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) {
00125 if (x0<-nx2-1 || y0<-ny2-1 || z0<-nz2-1 || x0>nx2 || y0>ny2 || z0>nz2 ) return false;
00126
00127
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
00139 float h=1.0f/EMConsts::I2G;
00140
00141 float r, gg;
00142
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
00148 gg = Util::fast_exp(-r *h)*weight;
00149
00150
00151
00152 size_t off;
00153 off=data->add_complex_at_fast(i,j,k,dt*gg);
00154
00155 norm[off/2]+=gg;
00156 }
00157 }
00158 }
00159 return true;
00160 }
00161 else {
00162
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) {
00192 if (x0<-nx2-2 || y0<-ny2-2 || z0<-nz2-2 || x0>nx2+1 || y0>ny2+1 || z0>nz2+1 ) return false;
00193
00194
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
00206
00207 float h=32.0f/((8.0f+Util::hypot3(xx,yy,zz))*EMConsts::I3G);
00208
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));
00210
00211 float r, gg;
00212
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
00218 gg = Util::fast_exp(-r *h)*w;
00219
00220
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) {
00242 if (x0<-nx2-4 || y0<-ny2-4 || z0<-nz2-4 || x0>nx2+3 || y0>ny2+3 || z0>nz2+3 ) return false;
00243
00244
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
00256
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));
00261
00262 float r, gg;
00263
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
00269 gg = Util::fast_exp(-r *h)*w;
00270
00271
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
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
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
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
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429 return true;
00430 }