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

This processor computes 2-D slices of the 4-D bispectrum of a 2-D image. More...

#include <processor.h>

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

Public Member Functions

 BispecSliceProcessor ()
 
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.bispectrum.slice"
 

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 2-D slices of the 4-D bispectrum of a 2-D image.

It can also integrate over image rotation to produce a set of rotationally/translationally invariant slices

Author
Steve Ludtke
Date
2017/07/04
Parameters
kxComplex x coordinate of the slice
kyComplex y coordinate of the slice
jkxInstead of kx,ky, specify the sum Jx+Kx for the slice
jkyInstead of kx,ky, specify the sum Jy+Ky for the slice
fpCompute a set of k=n rotational invariants and pack them into a single 2-D image stacked along y
kComplex radial coordinate of the slice, integrates over angle at this radius
sizespecifies the x.y dimension of the bispectrum image. Works by truncating high frequencies

Definition at line 734 of file processor.h.

Constructor & Destructor Documentation

◆ BispecSliceProcessor()

EMAN::BispecSliceProcessor::BispecSliceProcessor ( )
inline

Definition at line 737 of file processor.h.

737{}

Referenced by NEW().

Member Function Documentation

◆ get_desc()

string EMAN::BispecSliceProcessor::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 753 of file processor.h.

754 {
755 return "Computes a 2-D slice of the 4-D bispectrum of a 2-D image. Returns zero outside of source image. Input may be real or complex, but output is always complex. kx,ky OR jkx,jky OR k OR fp OR ffp";
756 }

◆ get_name()

string EMAN::BispecSliceProcessor::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 739 of file processor.h.

740 {
741 return NAME;
742 }
EMAN::BispecSliceProcessor::NAME
static const string NAME
Definition: processor.h:776

References NAME.

◆ get_param_types()

TypeDict EMAN::BispecSliceProcessor::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 758 of file processor.h.

759 {
760 TypeDict d;
761 d.put("kx", EMObject::INT, "Kx location of the slice in Fourier pixels");
762 d.put("ky", EMObject::INT, "Ky location of the slice in Fourier pixels");
763 d.put("jkx", EMObject::INT, "Jx+Kx location of the slice in Fourier pixels");
764 d.put("jky", EMObject::INT, "Jy+Ky location of the slice in Fourier pixels");
765 d.put("k", EMObject::FLOAT, "Radius of slice in Fourier pixels, integrates over angle.");
766 d.put("prbk", EMObject::FLOAT, "Radius of slice in Fourier pixels, integrates over angle.");
767 d.put("prbkv2", EMObject::FLOAT, "Radius of slice in Fourier pixels, integrates over angle. k>|q+k|>q");
768 d.put("prb3D", EMObject::FLOAT, "Radius of maximum slic. Returns volume, integrates over angle. k>|q+k|>q");
769 d.put("rfp", EMObject::INT, "Returns a non square 2-D image containing translatinal invariants organized such that X=azimuth. Used for rotational alignment.");
770 d.put("fp", EMObject::INT, "Returns a non-square 2-D image containing n rotationally integrated planes. R&T invariant.");
771 d.put("ffp", EMObject::INT, "Returns a 3-D volume containing n rotationally integrated planes. R&T invariant. This is normally further processed with CTF info to produce fp equivalent.");
772 d.put("size", EMObject::INT, "If specified, will determine the size (x/y) of the real-space bispectrum image. If not set, a size is selected automatically");
773 return d;
774 }
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::FLOAT, EMAN::EMObject::INT, and EMAN::TypeDict::put().

◆ NEW()

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

Definition at line 748 of file processor.h.

749 {
750 return new BispecSliceProcessor();
751 }

References BispecSliceProcessor().

◆ process()

EMData * BispecSliceProcessor::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 13324 of file processor.cpp.

13324 {
13325 if (image->get_zsize()!=1) throw ImageDimensionException("Only 2-D images supported");
13326
13327 EMData *cimage = NULL;
13328 if (image->is_complex()) cimage = image->copy();
13329 else cimage = image->do_fft();
13330 cimage->process_inplace("xform.phaseorigin.tocorner");
13331
13332 // Decide how large the bispectrum will be
13333 int nky=params.set_default("size",0)/2;
13334
13335 int nkx=nky+1;
13336 if (nky<4 || nky>=cimage->get_ysize()/2) nky=cimage->get_ysize()/8;
13337 if (nkx<5 || nky>=cimage->get_xsize()/2) nkx=cimage->get_xsize()/8+1;
13338 //if (image->get_xsize()<nky*2) throw ImageDimensionException("Image size smaller than requested footprint size, this is invalid.");
13339
13340 EMData* ret=new EMData(nkx*2,nky*2,1);
13341 ret->set_complex(1);
13342 ret->set_fftpad(1);
13343 ret->to_zero();
13344// printf("apix %f\tnx %d\n",(float)image->get_attr("apix_x"),(int)image->get_xsize());
13345
13346 // Fourier footprint mode, produces a 3-D image using n rotational invariants, real values from Fourier space
13347 if (params.has_key("ffp")) {
13348 // We need to save this so CTF can be applied accurately later on.
13349 float dsbg=1.0/(float(image->get_attr("apix_x"))*image->get_xsize()*4.0);
13350 int fp=(int)params.set_default("ffp",8);
13351 EMData *ret2=new EMData(nkx*2,nky*2,fp);
13352// ret2->to_zero();
13353 ret2->set_complex(1);
13354 ret2->set_ri(1);
13355 ret2->set_fftpad(1);
13356
13357 // We are doing a lot of rotations, so we make the image as small as possible first
13358 EMData *tmp=cimage;
13359 size_t laysize=nkx*2*nky*2;
13360 cimage=new EMData((nkx+fp+2)*2,(nky+fp+2)*2,1);
13361 cimage->set_complex(1);
13362 cimage->set_ri(1);
13363 cimage->set_fftpad(1);
13364
13365 for (int k=-cimage->get_ysize()/2; k<cimage->get_ysize()/2; k++) {
13366 for (int j=0; j<cimage->get_xsize()/2; j++) {
13367 cimage->set_complex_at(j,k,tmp->get_complex_at(j,k));
13368 }
13369 }
13370 delete tmp;
13371
13372 // now we compute the actual rotationally integrated "footprint"
13373 for (int k=2; k<2+fp; k++) {
13374// int jkx=k;
13375// int jky=0;
13376 int kx=k;
13377 int ky=0;
13378 ret->to_zero();
13379
13380 for (float ang=0; ang<360.0; ang+=360.0/(nky*M_PI) ) {
13381 EMData *cimage2=cimage->process("xform",Dict("alpha",ang));
13382
13383 for (int jy=-nky; jy<nky; jy++) {
13384 for (int jx=0; jx<nkx; jx++) {
13385 if (jx==0 && jy<0) continue; // this is critical! it avoids double-computation along kx=0
13386// int kx=jkx-jx;
13387// int ky=jky-jy;
13388 int jkx=jx+kx;
13389 int jky=jy+ky;
13390
13391// if (abs(jkx)>nkx || abs(jky)>nky) continue; // we go outside this range
13392 complex<double> v1 = (complex<double>)cimage2->get_complex_at(jx,jy);
13393 complex<double> v2 = (complex<double>)cimage2->get_complex_at(kx,ky);
13394 complex<double> v3 = (complex<double>)cimage2->get_complex_at(jkx,jky);
13395 ret->add_complex_at(jx,jy,0,(complex<float>)(v1*v2*std::conj(v3)));
13396 }
13397 }
13398 delete cimage2;
13399 }
13400
13401 // this fixes an issue with adding in the "special" Fourier locations ... sort of
13402// for (int jy=-nky; jy<nky; jy++) {
13403// ret->set_complex_at(0,jy,ret->get_complex_at(0,jy)/sqrt(2.0f));
13404// ret->set_complex_at(nkx-1,jy,ret->get_complex_at(nkx-1,jy)/sqrt(2.0f));
13405// }
13406 // simple fixed high-pass filter to get rid of gradient effects
13407 for (int jy=-2; jy<=2; jy++) {
13408 for (int jx=0; jx<=2; jx++) {
13409 ret->set_complex_at(jx,jy,0.0f);
13410 }
13411 }
13412 memcpy(ret2->get_data()+laysize*(k-2),ret->get_data(),laysize*sizeof(float));
13413 }
13414 delete ret;
13415 delete cimage;
13416 ret2->set_attr("dsbg",dsbg);
13417 return ret2;
13418 }
13419 // footprint mode, produces a 2-D image containing n veritcally arranged rotational invariants
13420 else if (params.has_key("fp")) {
13421 int fp=(int)params.set_default("fp",8);
13422 //EMData *ret2=new EMData(nkx*2-2,nky*2*fp,1);
13423 EMData *ret2=new EMData(nkx-1,nky*2*fp,1);
13424
13425 // We are doing a lot of rotations, so we make the image as small as possible first
13426 // note that 'step' is actually downsampling the image in Fourier space based
13427 // on the size of the desired bispectrum 9/1/18
13428 EMData *tmp=cimage;
13429// cimage=new EMData((nkx*2+fp+2)*2,(nky*2+fp+2)*2,1);
13430 cimage=new EMData((nkx*2+fp)*2-2,(nky*2+fp)*2,1);
13431 cimage->set_complex(1);
13432 cimage->set_ri(1);
13433 cimage->set_fftpad(1);
13434 int step=int(floor(tmp->get_ysize()/(2.0f*cimage->get_ysize())));
13435// int step=int(floor(tmp->get_ysize()/(cimage->get_ysize())));
13436 if (step==0) step=1;
13437// printf("%d\n",step);
13438
13439 for (int k=-cimage->get_ysize()/2; k<cimage->get_ysize()/2; k++) {
13440 for (int j=0; j<cimage->get_xsize()/2; j++) {
13441 cimage->set_complex_at(j,k,tmp->get_complex_at(j*step,k*step));
13442 }
13443 }
13444 delete tmp;
13445
13446// int mjkx=0,mjx=0,mkx=0;
13447 // now we compute the actual rotationally integrated "footprint"
13448 for (int dk=0; dk<fp; dk++) {
13449// int jkx=k;
13450// int jky=0;
13451// int kx=k;
13452// int ky=0;
13453 ret->to_zero();
13454
13455 for (float ang=0; ang<360.0; ang+=360.0/(nky*M_PI) ) {
13456 EMData *cimage2=cimage->process("xform",Dict("alpha",ang,"zerocorners",1));
13457
13458 for (int jy=-nky; jy<nky; jy++) {
13459 for (int jx=0; jx<nkx; jx++) {
13460 if ((jx==0 && jy<0) || (jx+dk<4)) continue; // first avoids double-computation along kx=0, second eliminates low resolution artifacts
13461// int kx=jkx-jx;
13462// int ky=jky-jy;
13463 int kx=jx+dk;
13464 int ky=0;
13465 int jkx=jx+kx;
13466 int jky=jy+ky;
13467// if (jkx>mjkx) { mjkx=jkx; mjx=jx; mkx=kx; }
13468 //int jkx2=jx-kx;
13469 //int jky2=jy-ky;
13470
13471// if (abs(jkx)>nkx || abs(jky)>nky) continue;
13472 complex<double> v1 = (complex<double>)cimage2->get_complex_at(jx,jy);
13473 complex<double> v2 = (complex<double>)cimage2->get_complex_at(kx,ky);
13474 complex<double> v3 = (complex<double>)cimage2->get_complex_at(jkx,jky);
13475// if (jx<4&&jy<4&&Util::hypot_fast(jx,jy)<4.0) v1=0.0; // we "filter" out the very low resolution
13476 ret->add_complex_at(jx,jy,0,(complex<float>)(v1*v2*std::conj(v3)));
13477
13478 //complex<double> v4 = (complex<double>)cimage2->get_complex_at(jkx2,jky2); // This is based on Phil's point that the v1,v2,v3 computation lacks Friedel symmetry
13479 //ret->add_complex_at(jx,jy,0,(complex<float>)(v1*(v2*std::conj(v3)+std::conj(v2)*std::conj(v4))));
13480
13481 }
13482 }
13483 delete cimage2;
13484 }
13485// printf("%d\t%d\t%d\t%d\n",dk,mjkx,mjx,mkx);
13486
13487 // this fixes an issue with adding in the "special" Fourier locations ... sort of
13488// for (int jy=-nky; jy<nky; jy++) {
13489// ret->set_complex_at(0,jy,ret->get_complex_at(0,jy)/sqrt(2.0f));
13490// ret->set_complex_at(nkx-1,jy,ret->get_complex_at(nkx-1,jy)/sqrt(2.0f));
13491// }
13492 // simple fixed high-pass filter to get rid of gradient effects
13493// for (int jy=-2; jy<=2; jy++) {
13494// for (int jx=0; jx<=2; jx++) {
13495// ret->set_complex_at(jx,jy,0.0f);
13496// }
13497// }
13498// ret->write_image("tst2.hdf",-1);
13499// EMData *pln=ret->do_ift();
13500// pln->process_inplace("xform.phaseorigin.tocenter");
13501// pln->process_inplace("normalize");
13502// ret2->insert_clip(ret,IntPoint(0,dk*nky*2,0));
13503 //ret2->insert_clip(pln,IntPoint(-nkx,(k-2)*nky*2,0));
13504// delete pln;
13505 for (int x=0; x<nkx-1; x++) {
13506 for (int y=-nky; y<nky; y++) {
13507 complex<float> val=ret->get_complex_at(x,y);
13508 if (x<=2&&abs(y)<=2&&Util::hypot_fast(x,y)<2.0) val=0.0; // We filter out the very low resolution
13509 ret2->set_value_at(x,y+nky+nky*dk*2,cbrt(std::real(val)));
13510// ret2->set_value_at(x,y+nky+nky*dk*2,std::real(val));
13511 }
13512 }
13513 }
13514// printf("%d\t%d\t%d\t%d\t%d\n",fp,nkx,nky,cimage->get_xsize(),cimage->get_ysize());
13515 delete ret;
13516 delete cimage;
13517 ret2->set_attr("is_bispec_fp",(int)fp);
13518 return ret2;
13519 }
13520 // In this mode we are making a translational invariant with information we can use for rotational alignment
13521 // It is not actually a bispectrum, but is based on a 2-point method instead
13522 else if (params.has_key("rfp")) {
13523// if (cimage->get_ysize()/2<nky*2) throw ImageDimensionException("Image size smaller than requested footprint size, this is invalid.");
13524 int rfp=(int)params.set_default("rfp",4);
13525 const int minr=4; // lowest frequency considered in the calculation (in pixels)
13526 int rsize=cimage->get_ysize()/4-minr-2;
13527 if (rsize>nky) rsize=nky;
13528// int nang=Util::calc_best_fft_size(int(M_PI*cimage->get_ysize())); // this gives basically 1 pixel precision at Nyquist. Accuracy is another question
13529 int nang=Util::calc_best_fft_size(int(M_PI*0.6*cimage->get_ysize())); // one pixel at 0.6 Nyquist. Refine alignment is probably necessary anyway
13530
13531 EMData *ret2=new EMData(nang,rsize*rfp*2,1);
13532 EMData *line=new EMData(minr*2+rsize*2,1,1); // one complex line
13533 line->set_complex(1);
13534 line->set_ri(1);
13535 line->set_fftpad(1); //correct?
13536
13537 for (int angi=0; angi<nang; angi++) { // this is x
13538 float ofs=M_PI/float(nang);
13539 float dx=cos(2.0*M_PI*angi/float(nang)+ofs);
13540 float dy=sin(2.0*M_PI*angi/float(nang)+ofs);
13541 line->to_zero();
13542
13543 // new invariant approach Steve 7/26/18
13544// for (int dr=0; dr<rfp; dr++) {
13545// for (int r=minr; r<rsize+minr; r++) {
13546// float jx=dx*r;
13547// float jy=dy*r;
13548// float kx=dx*(r+dr);
13549// float ky=dy*(r+dr);
13550//
13551// double pwr=double(r+dr)/double(r);
13552//
13553// complex<double> v1 = (complex<double>)cimage->get_complex_at_interp(jx,jy);
13554// complex<double> v2 = (complex<double>)cimage->get_complex_at_interp(kx,ky);
13555// line->set_complex_at(r,0,0,complex<float>(pow(v1,pwr)*std::conj(v2)));
13556// }
13557
13558
13559 // original bispectrum approach
13560 for (int dr=0; dr<rfp; dr++) {
13561 double lsq=0;
13562 for (int r=minr; r<rsize+minr; r++) {
13563 float jx=dx*r;
13564 float jy=dy*r;
13565 float kx=dx*(r+dr);
13566 float ky=dy*(r+dr);
13567
13568 // original bispectrum
13569 complex<double> v1 = (complex<double>)cimage->get_complex_at_interp(jx,jy);
13570 complex<double> v2 = (complex<double>)cimage->get_complex_at_interp(kx,ky);
13571 complex<double> v3 = (complex<double>)cimage->get_complex_at_interp(jx+kx,jy+ky);
13572 complex<double> bv = v1*v2*std::conj(v3);
13573 line->set_complex_at(r,0,0,complex<float>(bv));
13574 lsq+=std::norm(bv);
13575 }
13576 float rsig = 1.0/sqrt(lsq/double(rsize));
13577
13578 // Tried a bunch of different ways of representing the bispectral invariants, but the
13579 // pseudo real-space inverse seems to perform the best in testing on various targets
13580 // (decided this wasn't true in later testing)
13581
13582 // calling mult with 1D complex objects requires a lot of overhead for the RI/AP check, so we
13583 // combine it with output generation and compute rsig on the fly. Not exactly equivalent, but good enough.
13584// line->mult(1.0f/float(line->get_attr("sigma"))); // adjust intensities, but don't shift 0
13585
13586 // pseudo real-space
13587// EMData *lreal=line->do_ift();
13588// for (int i=0; i<rsize*2; i++) ret2->set_value_at(angi,i+(r2-minr)*rsize*2,lreal->get_value_at(i,0));
13589// for (int i=0; i<rsize*2; i++) ret2->set_value_at(angi,i+dr*rsize*2,lreal->get_value_at(i,0));
13590// delete lreal;
13591
13592 // stay in Fourier space
13593 // real/imaginary
13594// for (int i=0; i<rsize*2; i++) ret2->set_value_at(angi,i+(r2-minr)*rsize*2,line->get_value_at(i+minr*2,0));
13595 for (int i=0; i<rsize*2; i++) ret2->set_value_at(angi,i+dr*rsize*2,line->get_value_at(i+minr*2,0)*rsig);
13596
13597 // Amp/phase separation
13598// for (int i=0; i<rsize*2; i+=2) {
13599// complex<float> v(line->get_value_at(i,0),line->get_value_at(i+1,0));
13600// ret2->set_value_at(angi,i+ (r2-minr)*rsize*2,std::abs(v));
13601// ret2->set_value_at(angi,i+1+(r2-minr)*rsize*2,std::arg(v));
13602// }
13603
13604// // R/I with cube root
13605// for (int i=0; i<rsize*2; i+=2) {
13606// complex<float> v(line->get_value_at(i,0),line->get_value_at(i+1,0));
13607// v=pow(v,0.3333333f);
13608// ret2->set_value_at(angi,i+ (r2-minr)*rsize*2,std::real(v));
13609// ret2->set_value_at(angi,i+1+(r2-minr)*rsize*2,std::imag(v));
13610// }
13611
13612
13613
13614// printf("%d %d %d %d\n",rsize,angi,r2-minr,line->get_xsize());
13615 }
13616 }
13617 delete ret;
13618 delete line;
13619 delete cimage;
13620 ret2->set_attr("is_bispec_rfp",(int)rfp);
13621 return ret2;
13622 }
13623 // angular integrate mode, produces the simplest rotational invariant, ignoring rotational correlations
13624 else if (params.has_key("k")) {
13625 int k=(int)params.set_default("k",1);
13626// int jkx=k;
13627// int jky=0;
13628 int kx=k;
13629 int ky=0;
13630
13631 for (float ang=0; ang<360.0; ang+=360.0/(nky*M_PI) ) {
13632// for (float ang=0; ang<360.0; ang+=2 ) {// PRB changed this; change back
13633 EMData *cimage2=cimage->process("xform",Dict("alpha",ang));
13634
13635 for (int jy=-nky; jy<nky; jy++) {
13636 for (int jx=0; jx<nkx; jx++) {
13637// int kx=jkx-jx;
13638// int ky=jky-jy;
13639 int jkx=jx+kx;
13640 int jky=jy+ky;
13641
13642// if (abs(kx)>nkx || abs(ky)>nky) continue;
13643 complex<double> v1 = (complex<double>)cimage2->get_complex_at(jx,jy);
13644 complex<double> v2 = (complex<double>)cimage2->get_complex_at(kx,ky);
13645 complex<double> v3 = (complex<double>)cimage2->get_complex_at(jkx,jky);
13646 ret->add_complex_at(jx,jy,0,(complex<float>)(v1*v2*std::conj(v3)));
13647 }
13648 }
13649 delete cimage2;
13650 }
13651 }
13652 else if (params.has_key("prbk")) {
13653
13654 if (image->is_complex()) cimage = image->copy();
13655 else cimage = image->do_fft();
13656 cimage->process_inplace("xform.phaseorigin.tocorner");
13657
13658 // Decide how large the bispectrum will be
13659 //int nky=params.set_default("size",0)/2;
13660 int nky=params.set_default("size",0);
13661
13662 printf("nky = %d , size= %d \n",nky,image->get_xsize());
13663
13664 //int nkx=nky+1;
13665 int nkx=nky;
13666
13667 //cimage=new EMData((nkx*2+fp)*2-2,(nky*2+fp)*2,1);
13668 //cimage->set_complex(1);
13669 //cimage->set_ri(1);
13670 //cimage->set_fftpad(1);
13671
13672 //if (nky<4 || nky>=cimage->get_ysize()/2) nky=cimage->get_ysize()/8;
13673 //if (nkx<5 || nky>=cimage->get_xsize()/2) nkx=cimage->get_xsize()/8+1;
13674 //if (image->get_xsize()<nky*2) throw ImageDimensionException("Image size smaller than requested footprint size, this is invalid.");
13675
13676// EMData* ret=new EMData(nkx*2,nky*2,1);
13677 EMData* ret=new EMData(2*nkx,2*nky,1);
13678 ret->set_complex(0);
13679 //ret->set_fftpad(1);
13680 ret->to_zero();
13681
13682 int k=(int)params.set_default("prbk",1);
13683// int jkx=k;
13684// int jky=0;
13685 int kx=k;
13686 int ky=0;
13687
13688// for (float ang=0; ang<360.0; ang+=360.0/(nky*M_PI) ) { // Original
13689 for (float ang=0; ang<360.0; ang+=2 ) {
13690 EMData *cimage2=cimage->process("xform",Dict("alpha",ang));
13691
13692 for (int qy=-nky; qy<nky; qy++) {
13693 for (int qx=-nkx; qx<nkx; qx++) {
13694// int kx=jkx-jx;
13695// int ky=jky-jy;
13696 int kqx=qx+kx;
13697 int kqy=qy+ky;
13698 if (abs(kqx)>nkx || abs(kqy)>nky) continue;
13699
13700// if (abs(kx)>nkx || abs(ky)>nky) continue;
13701 complex<double> v1 = (complex<double>)cimage2->get_complex_at(kx,ky);
13702 complex<double> v2 = (complex<double>)cimage2->get_complex_at(qx,qy);
13703 complex<double> v3 = (complex<double>)cimage2->get_complex_at(kqx,kqy);
13704 complex<double> vTotal = (complex<double>)(v1*v2*std::conj(v3));
13705 double RealV = (real)(vTotal);
13706 int jxInd = (qx+2*nkx)%(2*nkx);
13707 int jyInd = (qy+2*nky)%(2*nky);
13708 //printf("RealV %g vTotal %g\n",RealV,vTotal);
13709 double OldVal = ret->get_value_at(jxInd,jyInd,0);
13710 ret->set_value_at(jxInd,jyInd,0,OldVal+RealV);
13711 //ret->set_value_at(jxInd,jyInd,0,1);
13712 }
13713 }
13714 delete cimage2;
13715 }
13716 return(ret);
13717 }
13718 else if (params.has_key("prbkv2")) {
13719
13720 if (image->is_complex()) cimage = image->copy();
13721 else cimage = image->do_fft();
13722 cimage->process_inplace("xform.phaseorigin.tocorner");
13723
13724 int k=(int)params.set_default("prbkv2",1);
13725
13726 // Decide how large the bispectrum will be
13727 int nky= 2*k+1;
13728 int nkx= 2*k+1;
13729
13730
13731 printf("nkx = %d, nky = %d , size= %d \n",nkx,nky,image->get_xsize());
13732
13733
13734 EMData* ret=new EMData(nkx,nky,1);
13735 ret->set_complex(0);
13736 ret->to_zero();
13737 int kx=k;
13738 int ky=0;
13739 int Nyquist = k;
13740// for (float ang=0; ang<360.0; ang+=360.0/(nky*M_PI) ) { // Original
13741 for (float ang=0; ang<360.0; ang+=2 ) {
13742 EMData *cimage2=cimage->process("xform",Dict("alpha",ang));
13743
13744 for (int qy=-k; qy<k+1; qy++) {
13745 if (qy*qy> (0.75)*k*k) continue;
13746 for (int qx=-k; qx<1; qx++) {
13747 if ((qx+k)*(qx+k) > (k*k - qy*qy)) continue;
13748 int kqx=qx+kx;
13749 int kqy=qy+ky;
13750 complex<double> v1 = (complex<double>)cimage2->get_complex_at(kx,ky);
13751 complex<double> v2 = (complex<double>)cimage2->get_complex_at(qx,qy);
13752 complex<double> v3 = (complex<double>)cimage2->get_complex_at(kqx,kqy);
13753 complex<double> vTotal = (complex<double>)(v1*v2*std::conj(v3));
13754 double RealV = (real)(vTotal);
13755 int qxInd = (qx+nkx)%(nkx);
13756 int qyInd = (qy+nky)%(nky);
13757 int NegQx = k/2 +qx ;
13758 int QyInd = qy+Nyquist;
13759 //printf("qy %d qyInd %d QyInd %d qx %d qxInd %d NegQx %d k %d\n",qy, qyInd, QyInd,qx,qxInd,NegQx, k);
13760 //printf("qx %d qy %d NegQx %d QyInd %d RealV %g vTotal %g\n",qx,qy, NegQx,QyInd, RealV,vTotal);
13761 //double OldVal = ret->get_value_at(qxInd,qyInd,0);
13762 double OldVal = ret->get_value_at(NegQx,QyInd,0);
13763 //ret->set_value_at(qxInd,qyInd,0,OldVal+RealV);
13764 ret->set_value_at(NegQx,QyInd,0,OldVal+RealV);
13765 //ret->set_value_at(qxInd,qyInd,0,1);
13766 }
13767 }
13768 delete cimage2;
13769 }
13770 return(ret);
13771 }
13772
13773 else if (params.has_key("prb3D")) {
13774
13775 int hp=(int)params.set_default("hp",0);
13776 int Nyquist=(int)params.set_default("prb3D",30);
13777
13778 if (image->is_complex()) cimage = image->copy();
13779 else cimage = image->do_fft();
13780
13781 cimage->process_inplace("xform.phaseorigin.tocorner");
13782
13783 //int k=(int)params.set_default("PRBkv3",1);
13784
13785 // Decide how large the bispectrum will be
13786
13787
13788 int csizex = cimage->get_xsize();
13789 int csizey = cimage->get_ysize();
13790
13791 printf("csizex = %d, csizey = %d \n", csizex , csizey );
13792
13793 Nyquist = (csizey+1)/2;
13794
13795 printf("csizex = %d, csizey = %d, Nyquist= %d \n", csizex , csizey, Nyquist);
13796
13797 int nkx= Nyquist/2+1;
13798 int nky= 2*Nyquist+2;
13799 int nkz= Nyquist+1;
13800
13801
13802 printf("nkx = %d, nky = %d , size= %d \n",nkx,nky,image->get_xsize());
13803
13804
13805
13806 EMData* ret=new EMData(nkx,nky,nkz);
13807 ret->set_complex(0);
13808 ret->to_zero();
13809
13810 //int Nyquist = 1* k;
13811// for (float ang=0; ang<360.0; ang+=360.0/(nky*M_PI) ) { // Original
13812 for (float ang=0; ang<360.0; ang+=2 ) {
13813 EMData *cimage2=cimage->process("xform",Dict("alpha",ang));
13814
13815 for (int k=0; k < Nyquist+1;k++){
13816 //printf(" k %d\n",k);
13817 int k2=k*k;
13818 int qyMax = sqrt(3*k2/4);
13819
13820 for (int qy=-qyMax; qy<qyMax+1; qy++) {
13821 //if (qy*qy> (0.75)*k*k) continue;
13822 int qxMin =k/2; qxMin *= -1;
13823 int qxMax = k - sqrt(k2-qy*qy)+0.999999; qxMax *= -1;
13824 for (int qx=qxMin; qx<qxMax+1; qx++) {
13825 //if ((qx+k)*(qx+k) > (k*k - qy*qy)) continue;
13826 //if (qx<(-k/2)) continue;
13827 int kqx=qx+k;
13828 int kqy=qy;
13829 complex<double> v1 = (complex<double>)cimage2->get_complex_at(k,0);
13830 complex<double> v2 = (complex<double>)cimage2->get_complex_at(qx,qy);
13831 complex<double> v3 = (complex<double>)cimage2->get_complex_at(kqx,kqy);
13832 complex<double> vTotal = (complex<double>)(v1*v2*std::conj(v3));
13833 double RealV = (real)(vTotal);
13834 int qxInd = (qx+nkx)%(nkx);
13835 int qyInd = (qy+nky)%(nky);
13836 int NegQx = k/2 +qx ;
13837 int QyInd = qy+Nyquist;
13838 //printf("qy %d qyInd %d QyInd %d qx %d qxInd %d NegQx %d k %d\n",qy, qyInd, QyInd,qx,qxInd,NegQx, k);
13839 //double OldVal = ret->get_value_at(qxInd,qyInd,0);
13840 double OldVal = ret->get_value_at(NegQx,QyInd,k);
13841 //ret->set_value_at(qxInd,qyInd,0,OldVal+RealV);
13842 ret->set_value_at(NegQx,QyInd,k,OldVal+RealV);
13843 //ret->set_value_at(qxInd,qyInd,0,1);
13844 }
13845 }
13846 }
13847 delete cimage2;
13848 }
13849 return(ret);
13850 }
13851 else if (params.has_key("jkx") && params.has_key("jky")) {
13852 int jkx=(int)params.set_default("jkx",0);
13853 int jky=(int)params.set_default("jky",0);
13854
13855 for (int jy=-nky; jy<nky; jy++) {
13856 for (int jx=0; jx<nkx; jx++) {
13857 int kx=jkx-jx;
13858 int ky=jky-jy;
13859
13860// if (abs(kx)>nkx || abs(ky)>nky) continue;
13861 complex<double> v1 = (complex<double>)cimage->get_complex_at(jx,jy);
13862 complex<double> v2 = (complex<double>)cimage->get_complex_at(kx,ky);
13863 complex<double> v3 = (complex<double>)cimage->get_complex_at(jkx,jky);
13864 ret->set_complex_at(jx,jy,(complex<float>)(v1*v2*std::conj(v3)));
13865// ret->set_complex_at(jx,jy,(complex<float>)(v1));
13866 }
13867 }
13868 }
13869 // normal slice mode
13870 else {
13871 int kx=(int)params.set_default("kx",0);
13872 int ky=(int)params.set_default("ky",0);
13873
13874// printf("KXY %d %d %d %d\n",kx,ky,nkx,nky);
13875
13876 for (int jy=-nky; jy<nky; jy++) {
13877 for (int jx=0; jx<nkx; jx++) {
13878// if (jx+kx>nkx || abs(jy+ky)>nky || jx+kx<0) continue;
13879 complex<double> v1 = (complex<double>)cimage->get_complex_at(jx,jy);
13880 complex<double> v2 = (complex<double>)cimage->get_complex_at(kx,ky);
13881 complex<double> v3 = (complex<double>)cimage->get_complex_at(jx+kx,jy+ky);
13882 ret->set_complex_at(jx,jy,(complex<float>)(v1*v2*std::conj(v3)));
13883// ret->set_complex_at(jx,jy,(complex<float>)(v1));
13884 }
13885 }
13886 }
13887
13888 for (int jy=-nky; jy<nky; jy++) {
13889 ret->set_complex_at(0,jy,ret->get_complex_at(0,jy)/sqrt(2.0f));
13890 ret->set_complex_at(nkx-1,jy,ret->get_complex_at(nkx-1,jy)/sqrt(2.0f));
13891 }
13892
13893 delete cimage;
13894
13895 return(ret);
13896}
EMAN::Dict
Dict is a dictionary to store <string, EMObject> pair.
Definition: emobject.h:385
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::EMData
EMData stores an image's data and defines core image processing routines.
Definition: emdata.h:82
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
sqrt
EMData * sqrt() const
return square root of current image
Definition: emdata_core.cpp:1084
real
EMData * real() const
return real part of a complex image as a real image format, if this image is a real image,...
get_ysize
int get_ysize() const
Get the image y-dimensional size.
Definition: emdata_metadata.h:553
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(), get_ysize(), EMAN::Dict::has_key(), EMAN::Util::hypot_fast(), ImageDimensionException, EMAN::Processor::params, real(), EMAN::Dict::set_default(), sqrt(), x, and y.

◆ process_inplace()

void EMAN::BispecSliceProcessor::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 744 of file processor.h.

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

References InvalidCallException.

Member Data Documentation

◆ NAME

const string BispecSliceProcessor::NAME = "math.bispectrum.slice"
static

Definition at line 776 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