EMAN2
Public Member Functions | Static Public Member Functions | Static Public Attributes | List of all members
EMAN::HarmonicProcessor Class Reference

This processor computes what I've dubbed the 'harmonic power spectrum'. More...

#include <processor.h>

Inheritance diagram for EMAN::HarmonicProcessor:
Inheritance graph
[legend]
Collaboration diagram for EMAN::HarmonicProcessor:
Collaboration graph
[legend]

Public Member Functions

 HarmonicProcessor ()
 
string get_name () const
 Get the processor's name. More...
 
void process_inplace (EMData *image)
 To process an image in-place. More...
 
virtual EMData * process (const EMData *const image)
 To proccess an image out-of-place. More...
 
string get_desc () const
 Get the descrition of this specific processor. More...
 
TypeDict get_param_types () const
 Get processor parameter information in a dictionary. More...
 
- Public Member Functions inherited from EMAN::Processor
virtual ~Processor ()
 
virtual void process_list_inplace (vector< EMData * > &images)
 To process multiple images using the same algorithm. More...
 
virtual Dict get_params () const
 Get the processor parameters in a key/value dictionary. More...
 
virtual void set_params (const Dict &new_params)
 Set the processor parameters using a key/value dictionary. More...
 

Static Public Member Functions

static Processor * NEW ()
 
- Static Public Member Functions inherited from EMAN::Processor
static string get_group_desc ()
 Get the description of this group of processors. More...
 
static void EMFourierFilterInPlace (EMData *fimage, Dict params)
 Compute a Fourier-filter processed image in place. More...
 
static EMData * EMFourierFilter (EMData *fimage, Dict params)
 Compute a Fourier-processor processed image without altering the original image. More...
 

Static Public Attributes

static const string NAME = "math.harmonic"
 

Additional Inherited Members

- Public Types inherited from EMAN::Processor
enum  fourier_filter_types {
  TOP_HAT_LOW_PASS , TOP_HAT_HIGH_PASS , TOP_HAT_BAND_PASS , TOP_HOMOMORPHIC ,
  GAUSS_LOW_PASS , GAUSS_HIGH_PASS , GAUSS_BAND_PASS , GAUSS_INVERSE ,
  GAUSS_HOMOMORPHIC , BUTTERWORTH_LOW_PASS , BUTTERWORTH_HIGH_PASS , BUTTERWORTH_HOMOMORPHIC ,
  KAISER_I0 , KAISER_SINH , KAISER_I0_INVERSE , KAISER_SINH_INVERSE ,
  SHIFT , TANH_LOW_PASS , TANH_HIGH_PASS , TANH_HOMOMORPHIC ,
  TANH_BAND_PASS , RADIAL_TABLE , CTF_
}
 Fourier filter Processor type enum. More...
 
- Protected Attributes inherited from EMAN::Processor
Dict params
 

Detailed Description

This processor computes what I've dubbed the 'harmonic power spectrum'.

It is a 2 point invariant method based on the product of a Fourier location and the complex conjugate of another Fourier location at r=1/N raised to the N power to cancel phases Since raising complex values to fracional powers causes root uncertainty issues, this insures that phase information is preserved Effectively this is giving a relative phase shift between different harmonics of each specific frequency. It also avoids problems with multiplying noisy components into the invariant as the bispectrum does. Anyway... that's the idea. Let's see if it works ;^)

Author
Steve Ludtke
Date
2019/03/17

Definition at line 681 of file processor.h.

Constructor & Destructor Documentation

◆ HarmonicProcessor()

EMAN::HarmonicProcessor::HarmonicProcessor ( )
inline

Definition at line 684 of file processor.h.

684{}

Referenced by NEW().

Member Function Documentation

◆ get_desc()

string EMAN::HarmonicProcessor::get_desc ( ) const
inlinevirtual

Get the descrition of this specific processor.

This function must be overwritten by a subclass.

Returns
The description of this processor.

Implements EMAN::Processor.

Definition at line 700 of file processor.h.

701 {
702 return "Computes invariants including relative phase in harmonic series";
703 }

◆ get_name()

string EMAN::HarmonicProcessor::get_name ( ) const
inlinevirtual

Get the processor's name.

Each processor is identified by a unique name.

Returns
The processor's name.

Implements EMAN::Processor.

Definition at line 686 of file processor.h.

687 {
688 return NAME;
689 }
EMAN::HarmonicProcessor::NAME
static const string NAME
Definition: processor.h:717

References NAME.

◆ get_param_types()

TypeDict EMAN::HarmonicProcessor::get_param_types ( ) const
inlinevirtual

Get processor parameter information in a dictionary.

Each parameter has one record in the dictionary. Each record contains its name, data-type, and description.

Returns
A dictionary containing the parameter info.

Reimplemented from EMAN::Processor.

Definition at line 705 of file processor.h.

706 {
707 TypeDict d;
708 d.put("hn", EMObject::INT, "Computes a single translational invariant for the nth harmonic, 1 is a normal power spectrum");
709 d.put("rn", EMObject::INT, "Computes a single rot/trans invariant for the nth rotational harmonic, requires hn to be non zero");
710 d.put("rfp", EMObject::INT, "Returns a non square 2-D image with translational invariants, y=radial, x=aziumth. Used for rotational alignment.");
711 d.put("fp", EMObject::INT, "Returns a non-square 2-D image containing n harmonics for each R&T component. R&T invariant. Min 2, default 4");
712// d.put("fb", EMObject::INT, "Fourier Bessel");
713 d.put("size", EMObject::INT, "If specified, will determine the number of rotational samples in the bispectrum. If not set, a size is selected automatically");
714 return d;
715 }
EMAN::TypeDict
TypeDict is a dictionary to store <string, EMObject::ObjectType> pair.
Definition: emobject.h:305
EMAN::TypeDict::put
void put(const string &key, EMObject::ObjectType o, const string &desc="")
Definition: emobject.h:330

References EMAN::EMObject::INT, and EMAN::TypeDict::put().

◆ NEW()

static Processor * EMAN::HarmonicProcessor::NEW ( )
inlinestatic

Definition at line 695 of file processor.h.

696 {
697 return new HarmonicProcessor();
698 }

References HarmonicProcessor().

◆ process()

EMData * HarmonicProcessor::process ( const EMData *const  image)
virtual

To proccess an image out-of-place.

For those processors which can only be processed out-of-place, override this function to give the right behavior.

Parameters
imageThe image will be copied, actual process happen on copy of image.
Returns
the image processing result, may or may not be the same size of the input image

Reimplemented from EMAN::Processor.

Definition at line 13898 of file processor.cpp.

13898 {
13899 if (image->get_zsize()!=1 || image->is_complex()) throw ImageDimensionException("Only 2-D real images supported");
13900
13901 EMData *cimage = NULL;
13902// if (image->is_complex()) cimage = image->copy();
13903// else cimage = image->do_fft();
13904 int s1=image->get_ysize();
13905 EMData *im2=image->get_clip(Region(-s1*3/2,-s1*3/2,s1*4,s1*4));
13906 cimage=im2->do_fft();
13907 delete im2;
13908 cimage->process_inplace("xform.phaseorigin.tocorner");
13909
13910 // Decide how large the harmonic invariant will be
13911// int nx=cimage->get_xsize();
13912 int ny=cimage->get_ysize()/8;
13913 int naz=Util::calc_best_fft_size((int)params.set_default("size",ny));
13914
13915 // Compute a translational invariant for a single harmonic
13916 EMData* trns=new EMData(naz*2,ny/2,1);
13917 trns->to_zero();
13918 trns->set_complex(1);
13919 trns->set_ri(1);
13920 trns->set_fftpad(1); // this is kind of meaningless in this context
13921 size_t xyz=trns->get_size();
13922// printf("apix %f\tnx %d\n",(float)image->get_attr("apix_x"),(int)image->get_xsize());
13923
13924 // Compute a translational invariant for a single harmonic
13925 if (params.has_key("hn")) {
13926 int hn=(int)params.get("hn");
13927
13928 // FFT Debugging code with hn<0
13929 if (hn<0) {
13930 complex<float> *tmp = (complex<float>*)EMfft::fftmalloc(naz*2);
13931 for (int jy=0; jy<ny/4; jy+=2) {
13932 for (int x=0; x<naz; x++) trns->set_complex_at_idx(x,jy,0,(complex<float>)std::polar((float)sin(1.0f+jy/2.0f*x*M_PI*2.0f/naz),-hn/100.0f));
13933 memcpy((void*)tmp,(void*)(trns->get_data()+jy*naz*2),naz*2*sizeof(float));
13934 EMfft::complex_to_complex_1d_inplace(tmp,naz*2);
13935 memcpy((void*)(trns->get_data()+(jy+1)*naz*2),(void*)tmp,naz*2*sizeof(float));
13936 }
13937 EMfft::fftfree((float *)tmp);
13938
13939 delete cimage;
13940 EMData *ret=trns->get_clip(Region(0,0,naz*2,ny/4));
13941 delete trns;
13942// ret->set_complex(0);
13943 return(ret);
13944 }
13945
13946 // Here is the actual hn code
13947 if (hn<1) throw InvalidParameterException("Invalid parameter, hn<1");
13948 int rn = 0;
13949 // rotational/translational single. If rn==0, special case where polar coordinate version of translational invariant is generated
13950 if (params.has_key("rn")) {
13951 rn=params.get("rn");
13952
13953 // Start with the translational invariant in Fourier space in a radial coordinate system
13954 for (int ja=0; ja<naz; ja++) {
13955 float si=sin(float(2.0*M_PI*(ja+0.5)/naz));
13956 float co=cos(float(2.0*M_PI*(ja+0.5)/naz));
13957 for (int jr=3; jr<ny/2 && jr*hn<ny*2; jr++) { // This is cryoEM specific, we have bad values near the origin
13958 float jx=co*jr;
13959 float jy=si*jr;
13960 complex<double> v2 = (complex<double>)cimage->get_complex_at_interp(jx,jy);
13961 complex<double> v1 = (complex<double>)cimage->get_complex_at_interp(jx*hn,jy*hn);
13962// int jx=co*jr;
13963// int jy=si*jr;
13964// complex<double> v2 = (complex<double>)cimage->get_complex_at(jx,jy);
13965// complex<double> v1 = (complex<double>)cimage->get_complex_at(jx*hn,jy*hn);
13966 trns->set_complex_at_idx(ja,jr,0,(complex<float>)(v1*std::pow(std::conj(v2),(float)hn)));
13967 }
13968 }
13969 // rescale components to have linear amplitude WRT the original FFT, without changing phase
13970 trns->ri2ap();
13971
13972 for (size_t i=0; i<xyz; i+=2) {
13973 trns->set_value_at_index(i,pow(trns->get_value_at_index(i),float(1.0/(hn+1))));
13974 }
13975 trns->ap2ri();
13976
13977 // Only if rn is defined
13978 if (rn==0) {
13979 complex<float> *tmp = (complex<float>*)EMfft::fftmalloc(naz*2);
13980 for (int jy=3; jy<ny/2; jy++) {
13981 memcpy((void*)tmp,(void*)(trns->get_data()+jy*naz*2),naz*2*sizeof(float));
13982 EMfft::complex_to_complex_1d_inplace(tmp,naz*2);
13983 memcpy((void*)(trns->get_data()+jy*naz*2),(void*)tmp,naz*2*sizeof(float));
13984 }
13985 EMfft::fftfree((float *)tmp);
13986 }
13987 else if (rn>=0) {
13988 // Now we do the 1-D FFTs on the lines of the translational invariant
13989 complex<float> *tmp = (complex<float>*)EMfft::fftmalloc(naz*2);
13990 for (int jy=3; jy<ny/2; jy++) {
13991 // While it might seem a good idea to do inplace 1D transforms for each row, the potential memory
13992 // alignment change for each row could cause bad things to happen
13993 memcpy((void*)tmp,(void*)(trns->get_data()+jy*naz*2),naz*2*sizeof(float));
13994 EMfft::complex_to_complex_1d_inplace(tmp,naz*2);
13995 for (int jx=0; jx<naz; jx++) {
13996// float l=jx/float(rn);
13997// float frc=l-floor(l);
13998// int li=(int)l;
13999// complex<float> v2 = Util::linear_interpolate_cmplx(tmp[li],tmp[li+1],frc);
14000 if (jx*rn>naz) {
14001 trns->set_complex_at_idx(jx,jy,0,0.0f);
14002 continue;
14003 }
14004 complex<float> v2 = tmp[jx];
14005 complex<float> v1 = tmp[jx*rn];
14006 trns->set_complex_at_idx(jx,jy,0,v1*std::pow(std::conj(v2),(float)rn));
14007 }
14008 // memcpy((void*)(trns->get_data()+jy*naz*2),(void*)tmp,naz*2*sizeof(float));
14009 }
14010 EMfft::fftfree((float *)tmp);
14011
14012 // rescale components to have linear amplitude WRT the original FFT, without changing phase
14013 trns->ri2ap();
14014 for (size_t i=0; i<xyz; i+=2) {
14015 trns->set_value_at_index(i,pow(trns->get_value_at_index(i),float(1.0/(rn+1))));
14016 }
14017 trns->ap2ri();
14018 }
14019 trns->set_attr("is_harmonic_rn",(int)rn);
14020 }
14021 // With no rotational component
14022 else {
14023 trns->set_size(ny,ny,1);
14024 xyz=trns->get_size();
14025 // translational only single
14026 for (int jx=0; jx<ny/2 && jx*hn<ny*2; jx++) {
14027 for (int jy=max(-ny/2,-ny*2/hn); jy<ny/2 && jy*hn<ny*2; jy++) {
14028 if (Util::hypot_fast(jx,jy)<2.5f) {
14029 trns->set_complex_at(jx,jy,0,(complex<float>)0);
14030 continue;
14031 }
14032 complex<double> v2 = (complex<double>)cimage->get_complex_at(jx,jy);
14033 complex<double> v1 = (complex<double>)cimage->get_complex_at(jx*hn,jy*hn);
14034 trns->set_complex_at(jx,jy,0,(complex<float>)(v1*std::pow(std::conj(v2),(float)hn)));
14035 }
14036 }
14037 // rescale components to have linear amplitude WRT the original FFT, without changing phase
14038 trns->ri2ap();
14039 for (size_t i=0; i<xyz; i+=2) {
14040 trns->set_value_at_index(i,pow(trns->get_value_at_index(i),float(1.0/(hn+1))));
14041 }
14042 trns->ap2ri();
14043 // EMData *ret = trns->do_ift();
14044 }
14045 delete cimage;
14046 trns->set_attr("is_harmonic_hn",(int)hn);
14047
14048 return(trns);
14049 }
14050 if (params.has_key("rfp")) {
14051 int rfp=(int)params.get("rfp");
14052 for (int ja=0; ja<naz; ja++) {
14053 float si=sin(float(2.0*M_PI*(ja+0.5)/naz));
14054 float co=cos(float(2.0*M_PI*(ja+0.5)/naz));
14055 int y=0; // This is where the results go
14056 // Start at r=5 due to bad CryoEM values near the origin.
14057 // Go to 1/2 Nyquist because high resolution values are less invariant due to interpolaton
14058 for (int jr=5; jr<ny/4; jr++) {
14059 // Innermost loop is hn (radial harmonic coefficient) to group similar values together
14060 for (int hn=2; hn<=rfp; hn++) {
14061 float jx=co*jr;
14062 float jy=si*jr;
14063 complex<double> v2 = (complex<double>)cimage->get_complex_at_interp(jx,jy);
14064 complex<double> v1 = (complex<double>)cimage->get_complex_at_interp(jx*hn,jy*hn);
14065 v1*=std::pow(std::conj(v2),(double)hn);
14066 v1=std::polar(std::pow(std::abs(v1),1.0/(hn+1.0))/ny,std::arg(v1));
14067 trns->set_complex_at_idx(ja,y,0,(complex<float>)v1);
14068 y++;
14069 if (y>=ny/2) break;
14070 }
14071 if (y>=ny/2) break;
14072 }
14073 // rescale components to have linear amplitude WRT the original FFT, without changing phase
14074 }
14075 delete cimage;
14076 trns->set_attr("is_harmonic_rfp",(int)rfp);
14077// trns->set_complex(0);
14078 return(trns);
14079 }
14080
14081 // Rotational & Translational invariant, fp specifies the maximum harmonic (R&T) to include
14082 if (params.has_key("fp")) {
14083 int fp=(int)params.set_default("fp",4);
14084 if (fp<2) fp=2;
14085 // Start with the translational invariant in Fourier space in a radial coordinate system
14086 for (int ja=0; ja<naz; ja++) {
14087 float si=sin(float(2.0*M_PI*(ja+0.5)/naz));
14088 float co=cos(float(2.0*M_PI*(ja+0.5)/naz));
14089 int y=0; // This is where the results go
14090 // Start at r=3 due to bad CryoEM values near the origin.
14091 // Go to 1/2 Nyquist because high resolution values are less invariant due to interpolaton
14092 for (int jr=3; jr<ny/4; jr++) {
14093 // Innermost loop is hn (radial harmonic coefficient) to group similar values together, skip the phaseless hn=1
14094 for (int hn=2; hn<=fp; hn++) {
14095 float jx=co*jr;
14096 float jy=si*jr;
14097 complex<double> v2 = (complex<double>)cimage->get_complex_at_interp(jx,jy);
14098 complex<double> v1 = (complex<double>)cimage->get_complex_at_interp(jx*hn,jy*hn);
14099 v1*=std::pow(std::conj(v2),(double)hn);
14100 v1=std::polar(std::pow(std::abs(v1),1.0/(hn+1.0))/ny,std::arg(v1));
14101 trns->set_complex_at_idx(ja,y,0,(complex<float>)v1);
14102 y++;
14103 if (y>=ny/2) break;
14104 }
14105 if (y>=ny/2) break;
14106 }
14107 // rescale components to have linear amplitude WRT the original FFT, without changing phase
14108 }
14109
14110 // This holds a copy for doing the 1D FFTs efficiently
14111 complex<float> *tmp = (complex<float>*)EMfft::fftmalloc(naz*2);
14112 // One horizontal line at a time
14113 for (int jy=0; jy<ny/2; jy++) {
14114 // Now we do the 1-D FFTs on the lines of the translational invariant
14115 memcpy((void*)tmp,(void*)(trns->get_data()+jy*naz*2),naz*2*sizeof(float));
14116 EMfft::complex_to_complex_1d_inplace(tmp,naz*2);
14117 int x=0; // output x coordinate
14118 // outer loop over base rotational frequency, skip phaseless jx=0
14119 for (int jx=1; jx<naz/2; jx++) {
14120 // inner loop over rotational harmonic coefficients, skip the phaseless rn=1
14121 complex<double> v2 = tmp[jx];
14122 for (int rn=2; rn<=fp; rn++) {
14123 if (jx*rn>=naz) break;
14124 complex<double> v1 = tmp[jx*rn];
14125 v1*=std::pow(std::conj(v2),(double)rn);
14126 v1=std::polar(std::pow(std::abs(v1),1.0/(rn+1.0))/naz,std::arg(v1));
14127 trns->set_complex_at_idx(x,jy,0,(complex<float>)v1);
14128 x++;
14129 if (x>=naz/2) break;
14130 }
14131 if (x>=naz/2) break;
14132 }
14133 }
14134 EMfft::fftfree((float *)tmp);
14135
14136
14137 delete cimage;
14138 // the /4 in the next line is arbitrary to remove regions which empirically
14139 // aren't useful
14140 EMData *ret=trns->get_clip(Region(0,0,min(naz,fp*naz/4)/2,min(ny/2,fp*(ny/4-3))));
14141 delete trns;
14142 ret->set_attr("is_harmonic_fp",(int)fp);
14143 ret->set_complex(0);
14144 return(ret);
14145 }
14146
14147}
EMAN::Dict::set_default
type set_default(const string &key, type val)
Default setting behavior This can be achieved using a template - d.woolford Jan 2008 (before there wa...
Definition: emobject.h:569
EMAN::Dict::has_key
bool has_key(const string &key) const
Ask the Dictionary if it as a particular key.
Definition: emobject.h:511
EMAN::Dict::get
EMObject get(const string &key) const
Get the EMObject corresponding to the particular key Probably better to just use operator[].
Definition: emobject.h:527
EMAN::EMData
EMData stores an image's data and defines core image processing routines.
Definition: emdata.h:82
EMAN::EMData::get_clip
EMData * get_clip(const Region &area, const float fill=0) const
Get an inclusive clip.
Definition: emdata.cpp:592
EMAN::Region
Region defines a 2D or 3D rectangular region specified by its origin coordinates and all edges' sizes...
Definition: geometry.h:497
EMAN::Util::calc_best_fft_size
static int calc_best_fft_size(int low)
Search the best FFT size with good primes.
Definition: util.cpp:1021
EMAN::Util::hypot_fast
static float hypot_fast(int x, int y)
Euclidean distance in 2D for integers computed fast using a cached lookup table.
Definition: util.cpp:742
InvalidParameterException
#define InvalidParameterException(desc)
Definition: exception.h:361
ImageDimensionException
#define ImageDimensionException(desc)
Definition: exception.h:166
y
#define y(i, j)
Definition: projector.cpp:1516
x
#define x(i)
Definition: projector.cpp:1517

References EMAN::Util::calc_best_fft_size(), EMAN::Dict::get(), EMAN::EMData::get_clip(), EMAN::Dict::has_key(), EMAN::Util::hypot_fast(), ImageDimensionException, InvalidParameterException, EMAN::Processor::params, EMAN::Dict::set_default(), x, and y.

◆ process_inplace()

void EMAN::HarmonicProcessor::process_inplace ( EMData *  image)
inlinevirtual

To process an image in-place.

For those processors which can only be processed out-of-place, override this function to just print out some error message to remind user call the out-of-place version.

Parameters
imageThe image to be processed.

Implements EMAN::Processor.

Definition at line 691 of file processor.h.

691{ throw InvalidCallException("inplace not supported"); }
InvalidCallException
#define InvalidCallException(desc)
Definition: exception.h:348

References InvalidCallException.

Member Data Documentation

◆ NAME

const string HarmonicProcessor::NAME = "math.harmonic"
static

Definition at line 717 of file processor.h.

Referenced by get_name().

The documentation for this class was generated from the following files:
Generated by doxygen 1.9.2