EMAN2
Classes | Macros | Enumerations | Functions
reconstructor.cpp File Reference
#include "reconstructor.h"
#include "plugins/reconstructor_template.h"
#include "ctf.h"
#include "emassert.h"
#include "symmetry.h"
#include <cstring>
#include <fstream>
#include <iomanip>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_fit.h>
#include <iostream>
#include <algorithm>
#include <sstream>
Include dependency graph for reconstructor.cpp:

Go to the source code of this file.

Classes

class  ctf_store_real
 

Macros

#define tw(i, j, k)   tw[ i-1 + (j-1+(k-1)*iy)*ix ]
 
#define tw(i, j, k)   tw[ i-1 + (j-1+(k-1)*iy)*ix ]
 

Enumerations

enum  weighting_method { NONE , ESTIMATE , VORONOI }
 

Functions

template<typename T >
void checked_delete (T *&x)
 
float max2d (int kc, const vector< float > &pow_a)
 
float max3d (int kc, const vector< float > &pow_a)
 
void circumfnn (EMData *win, int npad)
 
void circumftrl (EMData *win, int npad)
 
void printImage (const EMData *line)
 
void circumfnn_rect (EMData *win, int npad)
 

Macro Definition Documentation

◆ tw [1/2]

#define tw (   i,
  j,
 
)    tw[ i-1 + (j-1+(k-1)*iy)*ix ]

Definition at line 3694 of file reconstructor.cpp.

◆ tw [2/2]

#define tw (   i,
  j,
 
)    tw[ i-1 + (j-1+(k-1)*iy)*ix ]

Definition at line 3694 of file reconstructor.cpp.

Enumeration Type Documentation

◆ weighting_method

Enumerator
NONE 
ESTIMATE 
VORONOI 

Definition at line 2967 of file reconstructor.cpp.

2967{ NONE, ESTIMATE, VORONOI };
@ ESTIMATE
@ NONE
@ VORONOI

Function Documentation

◆ checked_delete()

template<typename T >
void checked_delete ( T *&  x)

◆ circumfnn()

void circumfnn ( EMData win,
int  npad 
)

Definition at line 3007 of file reconstructor.cpp.

3008{
3009 float *tw = win->get_data();
3010 // correct for the fall-off of NN interpolation using sinc functions
3011 // mask and subtract circumference average
3012 int ix = win->get_xsize();
3013 int iy = win->get_ysize();
3014 int iz = win->get_zsize();
3015 int L2 = (ix/2)*(ix/2);
3016 int L2P = (ix/2-1)*(ix/2-1);
3017
3018 int IP = ix/2+1;
3019 int JP = iy/2+1;
3020 int KP = iz/2+1;
3021
3022 // sinc functions tabulated for fall-off
3023 float* sincx = new float[IP+1];
3024 float* sincy = new float[JP+1];
3025 float* sincz = new float[KP+1];
3026
3027 sincx[0] = 1.0f;
3028 sincy[0] = 1.0f;
3029 sincz[0] = 1.0f;
3030
3031 float cor;
3032 if( npad == 1 ) cor = 1.0;
3033 else cor = 4.0;
3034
3035 float cdf = M_PI/(cor*ix);
3036 for (int i = 1; i <= IP; ++i) sincx[i] = sin(i*cdf)/(i*cdf);
3037 cdf = M_PI/(cor*iy);
3038 for (int i = 1; i <= JP; ++i) sincy[i] = sin(i*cdf)/(i*cdf);
3039 cdf = M_PI/(cor*iz);
3040 for (int i = 1; i <= KP; ++i) sincz[i] = sin(i*cdf)/(i*cdf);
3041 for (int k = 1; k <= iz; ++k) {
3042 int kkp = abs(k-KP);
3043 for (int j = 1; j <= iy; ++j) {
3044 cdf = sincy[abs(j- JP)]*sincz[kkp];
3045 for (int i = 1; i <= ix; ++i) tw(i,j,k) /= (sincx[abs(i-IP)]*cdf);
3046 }
3047 }
3048
3049 delete[] sincx;
3050 delete[] sincy;
3051 delete[] sincz;
3052
3053 float TNR = 0.0f;
3054 size_t m = 0;
3055 for (int k = 1; k <= iz; ++k) {
3056 for (int j = 1; j <= iy; ++j) {
3057 for (int i = 1; i <= ix; ++i) {
3058 size_t LR = (k-KP)*(k-KP)+(j-JP)*(j-JP)+(i-IP)*(i-IP);
3059 if (LR >= (size_t)L2P && LR<=(size_t)L2) {
3060 TNR += tw(i,j,k);
3061 ++m;
3062 }
3063 }
3064 }
3065 }
3066
3067 TNR /=float(m);
3068
3069
3070 for (int k = 1; k <= iz; ++k) {
3071 for (int j = 1; j <= iy; ++j) {
3072 for (int i = 1; i <= ix; ++i) {
3073 size_t LR = (k-KP)*(k-KP)+(j-JP)*(j-JP)+(i-IP)*(i-IP);
3074 if (LR<=(size_t)L2) tw(i,j,k) -= TNR;
3075 else tw(i,j,k) = 0.0f;
3076 }
3077 }
3078 }
3079
3080}
#define tw(i, j, k)

References tw.

Referenced by EMAN::nn4Reconstructor::finish(), and EMAN::nn4_ctfReconstructor::finish().

◆ circumfnn_rect()

void circumfnn_rect ( EMData win,
int  npad 
)

Definition at line 3695 of file reconstructor.cpp.

3696{
3697 float *tw = win->get_data();
3698 // correct for the fall-off
3699 // mask and subtract circumfnnerence average
3700 int ix = win->get_xsize();
3701 int iy = win->get_ysize();
3702 int iz = win->get_zsize();
3703
3704 int IP = ix/2+1;
3705 int JP = iy/2+1;
3706 int KP = iz/2+1;
3707
3708 // sinc functions tabulated for fall-off
3709 float* sincx = new float[IP+1];
3710 float* sincy = new float[JP+1];
3711 float* sincz = new float[KP+1];
3712
3713 sincx[0] = 1.0f;
3714 sincy[0] = 1.0f;
3715 sincz[0] = 1.0f;
3716
3717 float cdf = M_PI/float(npad*2*ix);
3718 for (int i = 1; i <= IP; ++i) sincx[i] = sin(i*cdf)/(i*cdf);
3719 cdf = M_PI/float(npad*2*iy);
3720 for (int i = 1; i <= JP; ++i) sincy[i] = sin(i*cdf)/(i*cdf);
3721 cdf = M_PI/float(npad*2*iz);
3722 for (int i = 1; i <= KP; ++i) sincz[i] = sin(i*cdf)/(i*cdf);
3723 for (int k = 1; k <= iz; ++k) {
3724 int kkp = abs(k-KP);
3725 for (int j = 1; j <= iy; ++j) {
3726 cdf = sincy[abs(j- JP)]*sincz[kkp];
3727 for (int i = 1; i <= ix; ++i) tw(i,j,k) /= (sincx[abs(i-IP)]*cdf);
3728 }
3729 }
3730
3731 delete[] sincx;
3732 delete[] sincy;
3733 delete[] sincz;
3734
3735
3736
3737 float dxx = 1.0f/float(0.25*ix*ix);
3738 float dyy = 1.0f/float(0.25*iy*iy);
3739
3740
3741
3742 float LR2=(float(ix)/2-1)*(float(ix)/2-1)*dxx;
3743
3744 float TNR = 0.0f;
3745 size_t m = 0;
3746 for (int k = 1; k <= iz; ++k) {
3747 for (int j = 1; j <= iy; ++j) {
3748 for (int i = 1; i <= ix; ++i) {
3749 float LR = (j-JP)*(j-JP)*dyy+(i-IP)*(i-IP)*dxx;
3750 if (LR<=1.0f && LR >= LR2) {
3751 TNR += tw(i,j,k);
3752 ++m;
3753 }
3754 }
3755 }
3756 }
3757
3758 TNR /=float(m);
3759
3760
3761 for (int k = 1; k <= iz; ++k) {
3762 for (int j = 1; j <= iy; ++j) {
3763 for (int i = 1; i <= ix; ++i) {
3764 float LR = (j-JP)*(j-JP)*dyy+(i-IP)*(i-IP)*dxx;
3765 if (LR<=1.0f) tw(i,j,k)=tw(i,j,k)-TNR;
3766 else tw(i,j,k) = 0.0f;
3767 }
3768 }
3769 }
3770
3771}

References tw.

Referenced by EMAN::nn4_rectReconstructor::finish(), and EMAN::nn4_ctf_rectReconstructor::finish().

◆ circumftrl()

void circumftrl ( EMData win,
int  npad 
)

Definition at line 3083 of file reconstructor.cpp.

3084{
3085 float *tw = win->get_data();
3086 // correct for the fall-off of tri-linear interpolation using sinc^2 functions
3087 // mask and subtract circumference average
3088 int ix = win->get_xsize();
3089 int iy = win->get_ysize();
3090 int iz = win->get_zsize();
3091 int L2 = (ix/2)*(ix/2);
3092 int L2P = (ix/2-1)*(ix/2-1);
3093
3094 int IP = ix/2+1;
3095 int JP = iy/2+1;
3096 int KP = iz/2+1;
3097
3098 // sinc functions tabulated for fall-off
3099 float* sincx = new float[IP+1];
3100 float* sincy = new float[JP+1];
3101 float* sincz = new float[KP+1];
3102
3103 sincx[0] = 1.0f;
3104 sincy[0] = 1.0f;
3105 sincz[0] = 1.0f;
3106
3107 float cor;
3108 if( npad == 1 ) cor = 1.0;
3109 else cor = 4.0;
3110
3111 float cdf = M_PI/(cor*ix);
3112 for (int i = 1; i <= IP; ++i) sincx[i] = pow(sin(i*cdf)/(i*cdf),2);
3113 cdf = M_PI/(cor*iy);
3114 for (int i = 1; i <= JP; ++i) sincy[i] = pow(sin(i*cdf)/(i*cdf),2);
3115 cdf = M_PI/(cor*iz);
3116 for (int i = 1; i <= KP; ++i) sincz[i] = pow(sin(i*cdf)/(i*cdf),2);
3117 for (int k = 1; k <= iz; ++k) {
3118 int kkp = abs(k-KP);
3119 for (int j = 1; j <= iy; ++j) {
3120 cdf = sincy[abs(j- JP)]*sincz[kkp];
3121 for (int i = 1; i <= ix; ++i) tw(i,j,k) /= (sincx[abs(i-IP)]*cdf);
3122 }
3123 }
3124
3125 delete[] sincx;
3126 delete[] sincy;
3127 delete[] sincz;
3128
3129 float TNR = 0.0f;
3130 size_t m = 0;
3131 for (int k = 1; k <= iz; ++k) {
3132 for (int j = 1; j <= iy; ++j) {
3133 for (int i = 1; i <= ix; ++i) {
3134 size_t LR = (k-KP)*(k-KP)+(j-JP)*(j-JP)+(i-IP)*(i-IP);
3135 if (LR >= (size_t)L2P && LR<=(size_t)L2) {
3136 TNR += tw(i,j,k);
3137 ++m;
3138 }
3139 }
3140 }
3141 }
3142
3143 TNR /=float(m);
3144
3145
3146 for (int k = 1; k <= iz; ++k) {
3147 for (int j = 1; j <= iy; ++j) {
3148 for (int i = 1; i <= ix; ++i) {
3149 size_t LR = (k-KP)*(k-KP)+(j-JP)*(j-JP)+(i-IP)*(i-IP);
3150 if (LR<=(size_t)L2) tw(i,j,k) -= TNR;
3151 else tw(i,j,k) = 0.0f;
3152 }
3153 }
3154 }
3155
3156}

References tw.

Referenced by EMAN::nn4_ctfwReconstructor::finish(), and EMAN::nn4_ctfwsReconstructor::finish().

◆ max2d()

float max2d ( int  kc,
const vector< float > &  pow_a 
)

Definition at line 2969 of file reconstructor.cpp.

2970{
2971 float max = 0.0;
2972 for( int i=-kc; i <= kc; ++i ) {
2973 for( int j=-kc; j <= kc; ++j ) {
2974 if( i==0 && j==0 ) continue;
2975 {
2976 int c = 2*kc+1 - std::abs(i) - std::abs(j);
2977 max = max + pow_a[c];
2978 }
2979 }
2980 }
2981 return max;
2982}

Referenced by EMAN::nn4Reconstructor::finish(), and EMAN::nn4_rectReconstructor::finish().

◆ max3d()

float max3d ( int  kc,
const vector< float > &  pow_a 
)

Definition at line 2984 of file reconstructor.cpp.

2985{
2986 float max = 0.0;
2987 for( int i=-kc; i <= kc; ++i ) {
2988 for( int j=-kc; j <= kc; ++j ) {
2989 for( int k=-kc; k <= kc; ++k ) {
2990 if( i==0 && j==0 && k==0 ) continue;
2991 // if( i!=0 )
2992 {
2993 int c = 3*kc+1 - std::abs(i) - std::abs(j) - std::abs(k);
2994 max = max + pow_a[c];
2995 // max = max + c * c;
2996 // max = max + c;
2997 }
2998 }
2999 }
3000 }
3001 return max;
3002}

Referenced by EMAN::nn4Reconstructor::finish(), EMAN::nn4_rectReconstructor::finish(), EMAN::nnSSNR_Reconstructor::finish(), EMAN::nn4_ctfReconstructor::finish(), EMAN::nn4_ctf_rectReconstructor::finish(), and EMAN::nnSSNR_ctfReconstructor::finish().

◆ printImage()

void printImage ( const EMData line)

Definition at line 3271 of file reconstructor.cpp.

3272{
3273 Assert( line->get_zsize()==1 );
3274
3275
3276 int nx = line->get_xsize();
3277 int ny = line->get_ysize();
3278 for( int j=0; j < ny; ++j ) {
3279 for( int i=0; i < nx; ++i ) printf( "%10.3f ", line->get_value_at(i,j) );
3280 printf( "\n" );
3281 }
3282}
#define Assert(s)
Define Assert() function that is effective only when -DDEBUG is used.
Definition: emassert.h:42

References Assert.