Real space 3D reconstruction using per-voxel median. More...

#include <reconstructor.h>

Inheritance diagram for EMAN::RealMedianReconstructor:

[legend]

Collaboration diagram for EMAN::RealMedianReconstructor:

[legend]

Public Member Functions
	RealMedianReconstructor ()

virtual	~RealMedianReconstructor ()

virtual void	setup ()
	Initialize the reconstructor. More...

virtual EMData *	preprocess_slice (const EMData *const slice, const Transform &t=Transform())
	While you can just insert unprocessed slices, if you call preprocess_slice yourself, and insert the returned slice instead, repeatedly, it can save a fair bit of computation. More...

virtual int	insert_slice (const EMData *const slice, const Transform &euler, const float weight)
	Insert an image slice to the reconstructor. More...

virtual int	determine_slice_agreement (EMData *slice, const Transform &euler, const float weight=1.0, bool sub=true)
	Dummy function which always returns the same values. More...

virtual EMData *	finish (bool doift=true)
	Finish reconstruction and return the complete model. More...

virtual string	get_name () const
	Get the unique name of this class (especially for factory based instantiation access) More...

virtual string	get_desc () const
	Get a clear, concise description of this class. More...

virtual TypeDict	get_param_types () const

Public Member Functions inherited from EMAN::Reconstructor
	Reconstructor ()

virtual	~Reconstructor ()

virtual void	setup_seed (EMData *seed, float seed_weight)
	Initialize the reconstructor with a seed volume. More...

virtual void	setup_seedandweights (EMData seed, EMData weight)
	Initialize the reconstructor with a seed volume, as above. More...

int	insert_slice (const EMData *const slice, const Transform &euler)

virtual EMData *	projection (const Transform &euler, int ret_fourier=1)
	This will create a projection from the current reconstruction. More...

virtual void	clear ()
	set the volume and tmp_volume data to zero, for use in Monte Carlo reconstructors More...

void	print_params () const
	Print the current parameters to std::out. More...

EMObject &	operator[] (const string &key)

Public Member Functions inherited from EMAN::FactoryBase
	FactoryBase ()

virtual	~FactoryBase ()

Dict	get_params () const
	get a copy of the parameters of this class More...

void	set_params (const Dict &new_params)
	Set new parameters. More...

void	set_param (const string key, const EMObject val)

void	insert_params (const Dict &new_params)
	Insert parameters. More...

Dict	copy_relevant_params (const FactoryBase *const that) const

Public Member Functions inherited from EMAN::ReconstructorVolumeData
	ReconstructorVolumeData ()
	Only constructor All member variables are zeroed. More...

virtual	~ReconstructorVolumeData ()
	Destructor safely frees memory. More...

const EMData *	get_emdata ()
	Get the main image pointer, probably redundant (not used) More...

Static Public Member Functions
static Reconstructor *	NEW ()

Static Public Attributes
static const string	NAME = "real_median"

Private Member Functions
	RealMedianReconstructor (const RealMedianReconstructor &that)

RealMedianReconstructor &	operator= (const RealMedianReconstructor &)

Private Attributes
vector< EMData * >	slices

vector< Transform >	xforms

Additional Inherited Members
Protected Member Functions inherited from EMAN::ReconstructorVolumeData
void	free_memory ()
	Free allocated memorys The inherited class may have allocated image of tmp_data In either case you can safely call this function to delete either of those pointers, even if they bdb:refine_03::threed_00are NULL. More...

virtual void	normalize_threed (const bool sqrt_damp=false, const bool wiener=false)
	Normalize on the assumption that image is a Fourier volume and that tmp_data is a volume of weights corresponding in size to this Fourier volume. More...

virtual void	zero_memory ()
	Sends the pixels in tmp_data and image to zero Convenience only. More...

Protected Attributes inherited from EMAN::FactoryBase
Dict	params
	This is the dictionary the stores the parameters of the object. More...

Protected Attributes inherited from EMAN::ReconstructorVolumeData
EMData *	image
	Inheriting class allocates this, probably in setup(). More...

EMData *	tmp_data
	Inheriting class may allocate this, probably in setup() More...

int	nx

int	nx2

int	ny

int	ny2

int	nz

int	nz2

int	subnx

int	subny

int	subnz

int	subx0

int	suby0

int	subz0

Detailed Description

Real space 3D reconstruction using per-voxel median.

Real-space reconstruction similar to unfiltered back-projection, but computing the median rather than the mean value.

Definition at line 921 of file reconstructor.h.

Constructor & Destructor Documentation

◆ RealMedianReconstructor() [1/2]

EMAN::RealMedianReconstructor::RealMedianReconstructor ( )

inline

Definition at line 924 of file reconstructor.h.

924{ }

Referenced by NEW().

◆ ~RealMedianReconstructor()

virtual EMAN::RealMedianReconstructor::~RealMedianReconstructor ( )

inlinevirtual

Definition at line 926 of file reconstructor.h.

926{}

◆ RealMedianReconstructor() [2/2]

EMAN::RealMedianReconstructor::RealMedianReconstructor ( const RealMedianReconstructor & that )

private

Member Function Documentation

◆ determine_slice_agreement()

int RealMedianReconstructor::determine_slice_agreement	(	EMData *	slice,
		const Transform &	euler,
		const float	weight = `1.0`,
		bool	sub = `true`
	)

virtual

Dummy function which always returns the same values.

Parameters

input_slice	The EMData slice to be compared
euler	The orientation of the slice as a Transform object
weight	A weighting factor for this slice, generally the number of particles in a class-average. May be ignored by some reconstructors
sub	Flag indicating whether to subtract the slice from the volume before comparing. May be ignored by some reconstructors

Returns: 0 if OK. 1 if error.

Exceptions

NullPointerException	if the input EMData pointer is null
ImageFormatException	if the image is complex as opposed to real

Reimplemented from EMAN::Reconstructor.

Definition at line 2569 of file reconstructor.cpp.

{
        // Are these exceptions really necessary? (d.woolford)
        if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
 
        input_slice->set_attr("reconstruct_norm",1.0f);
        input_slice->set_attr("reconstruct_absqual",1.0f);
        input_slice->set_attr("reconstruct_weight",1.0f);
 
        return 0;
 
}

References NullPointerException.

◆ finish()

EMData * RealMedianReconstructor::finish ( bool doift = true )

virtual

Finish reconstruction and return the complete model.

Parameters

doift A flag indicating whether the returned object should be guaranteed to be in real-space (true) or should be left in whatever space the reconstructor generated

Returns: The result 3D model.

Reimplemented from EMAN::Reconstructor.

Definition at line 2582 of file reconstructor.cpp.

{
        
        int mode = params.set_default("mode",0);        // confusing, mode means mode of operation
 
        // We do the actual reconstruction here from all of the slices
        for (int z=0; z<nz; z++) {
                for (int y=0; y<ny; y++) {
                        for (int x=0; x<nx; x++) {
                                std::vector<float> vals;
                                std::vector<Transform>::iterator itt = xforms.begin();
                                for (std::vector<EMData *>::iterator it = slices.begin() ; it != slices.end(); ++it,++itt) {
                                        Vec3f imvec = itt->transform(x-nx/2.0f,y-ny/2.0f,z-nz/2.0f);
                                        
                                        // single interpolated value from 2-D
                                        float val=(*it)->sget_value_at_interp(imvec[0]+nx/2,imvec[1]+ny/2);
                                        if (val!=0) vals.push_back(val);
                                        
                                        // uses all 4 of the nearest values from the 2-D image
//                                      int xx=floor(imvec[0]+nx/2);
//                                      int yy=floor(imvec[1]+ny/2);
//                                      float val=(*it)->sget_value_at(xx,yy);
//                                      if (val!=0) vals.push_back(val);
//                                      val=(*it)->sget_value_at(xx+1,yy);
//                                      if (val!=0) vals.push_back(val);
//                                      val=(*it)->sget_value_at(xx+1,yy+1);
//                                      if (val!=0) vals.push_back(val);
//                                      val=(*it)->sget_value_at(xx,yy+1);
//                                      if (val!=0) vals.push_back(val);
                                }
                                if (vals.size()==0) continue;
                                
//                              // debugging, save all values for specific voxels to look at statistics
//                              if ((x==515&&y==461&&z==281) ||
//                              (x==516&&y==462&&z==299) ||
//                              (x==516&&y==461&&z==262) ||
//                              (x==212&&y==409&&z==304)) {
//                                      char fsp[80];
//                                      sprintf(fsp,"stat_%d_%d_%d.txt",x,y,z);
//                                      FILE *out=fopen(fsp,"w");
//                                      for (std::vector<float>::iterator v = vals.begin(); v!=vals.end(); ++v) fprintf(out,"%f\n",*v);
//                                      fclose(out);
//                              }
                                
                                switch(mode) {
                                // median, the namesake of the reconstructor
                                case 0:
                                        std::sort(vals.begin(),vals.end());
                                        //printf("%d %d %d    %d\n",x,y,z,vals.size());
                                        image->set_value_at(x,y,z,vals[vals.size()/2]);         // for even sizes, not quite right, and should include possibility of local average
                                        break;
 
                                // This is something like a mode. It progressively subdivides the numeric axis in halves
                                // until reaching a small number of values which are then averaged
                                case 1:
                                        {
                                        // We only use this method if we have at least 6 possible values for the voxel, otherwise we just average
                                        if (vals.size()<6) {
                                                float val=0;
                                                for (int i=0; i<vals.size(); i++) val+=vals[i];
                                                image->set_value_at(x,y,z,val/vals.size());
                                                break;
                                        }
                                        // Sort the list to simplify the remaining logic
                                        std::sort(vals.begin(),vals.end());
                                        int v0=0,v1=vals.size(),vs=v1/2;                // bracket values, v0 first, v1 last+1, vs is first value in second zone
                                        
                                        // we continue until we have a zone containing 4 or fewer values
                                        //int count=0;
                                        while (v1-v0>8) {
                                                float cen=(vals[v0]+vals[v1-1])/2.0f;           // value halfway between the extremes
                                                for (vs=v0; vs<v1; vs++) { if (vals[vs]>=cen) break; }          // find the first value in the upper zone
                                                // shift to the zone containing more values
                                                if (v1-vs>vs-v0) v0=vs;
                                                else v1=vs;
                                                //printf("%d %d %d %d\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n",v0,vs,v1,vals.size(),cen,vals[0],vals[v0],vals[vs-1],vals[vs],vals[v1],vals[vals.size()-1]);
                                                //count++;
                                                //if (count>20) exit(1);
                                        }
                                        float val=0;
                                        for (int i=v0; i<v1; i++) val+=vals[i];
                                        image->set_value_at(x,y,z,val/(v1-v0));
                                        break;
                                        }
                                // median with a bias towards zero (positive and negative), using 1/3 or 2/3 position in the sorted array
                                case 2:
                                {
                                        if (vals.size()<6) {
                                                float val=0;
                                                for (int i=0; i<vals.size(); i++) val+=vals[i];
                                                image->set_value_at(x,y,z,val/vals.size());
                                                break;
                                        }
                                        std::sort(vals.begin(),vals.end());
                                        //printf("%d %d %d    %d\n",x,y,z,vals.size());
                                        int p1=vals.size()/3-1;
                                        int p2=vals.size()*2/3;
                                        float v=(fabs(vals[p1])<fabs(vals[p2]))?(vals[p1]+vals[p1-1]+vals[p1+1])/3.0f:(vals[p2]+vals[p2-1]+vals[p2+1])/3.0f;
                                        image->set_value_at(x,y,z,v);           // for even sizes, not quite right, and should include possibility of local average
                                        break;
                                }
                                // bias the median towards the half of the points spanning the narrower range
                                case 3:
                                {
                                        if (vals.size()<6) {
                                                float val=0;
                                                for (int i=0; i<vals.size(); i++) val+=vals[i];
                                                image->set_value_at(x,y,z,val/vals.size());
                                                break;
                                        }
                                        std::sort(vals.begin(),vals.end());
                                        //printf("%d %d %d    %d\n",x,y,z,vals.size());
                                        int p1=vals.size()/3-1;
                                        int p2=vals.size()/2;
                                        int p3=vals.size()*2/3;
                                        float v=(vals[p2]-vals[0]<vals[vals.size()-1]-vals[p2])?(vals[p1]+vals[p1-1]+vals[p1+1])/3.0f:(vals[p3]+vals[p3-1]+vals[p3+1])/3.0f;
                                        image->set_value_at(x,y,z,v);           // for even sizes, not quite right, and should include possibility of local average
                                        break;
                                }
                                case 4:
                                {
                                        // too few values, return median
                                        if (vals.size()<10) {
                                                std::sort(vals.begin(),vals.end());
                                                image->set_value_at(x,y,z,vals[vals.size()/2]);         // for even sizes, not quite right, and should include possibility of local average
                                                break;
                                        }
                                        
                                        //compute stats
                                        int n=vals.size();
                                        float mean=0,sigma=0;
                                        for (int i=0; i<n; i++) { mean+=vals[i]; sigma+=vals[i]*vals[i]; }
                                        mean/=n;
                                        sigma=sqrt(sigma/n-mean);
                                        // low & high refer to angle not value
                                        float low=(vals[0]+vals[1]+vals[2]+vals[3])/4.0f;
                                        float high=(vals[n-1]+vals[n-2]+vals[n-3]+vals[n-4])/4.0f;
                                        
                                        // If both extreme tilts roughly agree, then we assume that's the right range to use and we average values near that
                                        // by computing the mean of all values within 1 sigma of all values
                                        if (fabs(low-high)<sigma) {
                                                mean=(low+high)*4.0f;
                                                int c=8;
                                                for (int i=4; i<n-4; i++) { 
                                                        if (fabs(vals[i]-mean/c)<sigma) { mean+=vals[i]; c++; } 
                                                }
                                                image->set_value_at(x,y,z,mean/c);
                                        }
                                        // This means the extremes disagree and we have to pick which value is best, we go with min(fabs)
                                        else {
                                                if (fabs(low)<fabs(high)) mean=low*4.0;
                                                else mean=high*4.0;
                                                int c=4;
                                                for (int i=4; i<n-4; i++) { 
                                                        if (fabs(vals[i]-mean/c)<sigma) { mean+=vals[i]; c++; } 
                                                }
                                                image->set_value_at(x,y,z,mean/c);
                                        }
                                        break;
                                }
                                }
                        }
                }
        } 
 
 
//      transform->set_scale(1.0);
//      transform->set_mirror(false);
//      transform->set_trans(0,0,0);
//      transform->invert();
 
//      tmp->transform(t.inverse());            // This was incorrect. inverse() was missing
//      image->add(*tmp);
 
//      if(transform) {delete transform; transform=0;}
 
        // apply symmetry if requested
        Symmetry3D* sym = Factory<Symmetry3D>::get((string)params["sym"]);
        vector<Transform> syms = sym->get_syms();
 
        for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
 
//              it->printme();
                Transform t=*it;
                EMData *tmpcopy = image->process("xform",Dict("transform",(EMObject)&t));
                image->add(*tmpcopy);
                delete tmpcopy;
        }
 
        image->mult(1.0f/(float)sym->get_nsym());
        delete sym;
//      if (image->get_xsize()==image->get_ysize() && image->get_ysize()==image->get_zsize()) {
//              image->process_inplace("mask.sharp",Dict("outer_radius",image->get_xsize()/2-1));
//      }
//      else printf("No masking %d %d %d\n",image->get_xsize(),image->get_ysize(),image->get_zsize());
        
        // erase out copy of the projection data
        if (!slices.empty()) {
                for (std::vector<EMData *>::iterator it = slices.begin() ; it != slices.end(); ++it) delete *it;
                slices.clear();
        }
        xforms.clear();
        
        EMData *ret = image;
        if (!doift) {
                ret=image->do_fft();
                ret->process_inplace("xform.phaseorigin.tocorner");
                delete image;
        }
        image = 0 ;
        return ret;
}

References EMAN::Factory< T >::get(), EMAN::Symmetry3D::get_nsym(), EMAN::Symmetry3D::get_syms(), EMAN::ReconstructorVolumeData::image, EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::nz, EMAN::FactoryBase::params, EMAN::Dict::set_default(), slices, sqrt(), v0, x, xforms, and y.

◆ get_desc()

virtual string EMAN::RealMedianReconstructor::get_desc ( ) const

inlinevirtual

Get a clear, concise description of this class.

Returns: a clear, concise description of this class

Implements EMAN::FactoryBase.

Definition at line 968 of file reconstructor.h.

                {
                        return "A back projection reconstructor with alternative statistics. mode:\n\
0 - median rather than mean per voxel (default)\n\
1 - similar to a mode. progressive subdivision of numeric axis to find dense value cluster\n\
2 - median with a bias towands zero (postive and negative)\n\
3 - median with a bias towards the 1/2 of the values spanning a narrower range\n\
4 - specifically for tomography, requires slices to be in order by tilt\n\
";
                }

◆ get_name()

virtual string EMAN::RealMedianReconstructor::get_name ( ) const

inlinevirtual

Get the unique name of this class (especially for factory based instantiation access)

Returns: the unique name of this class

Implements EMAN::FactoryBase.

Definition at line 963 of file reconstructor.h.

                {
                        return NAME;
                }

References NAME.

◆ get_param_types()

virtual TypeDict EMAN::RealMedianReconstructor::get_param_types ( ) const

inlinevirtual

Returns: a TypeDict defining and describing the feasible parameters of this class

Implements EMAN::FactoryBase.

Definition at line 984 of file reconstructor.h.

                {
                        TypeDict d;
                        d.put("mode", EMObject::INT, "Default is 0, median. 1 - compute estimated statistical mode.");
                        d.put("size", EMObject::INTARRAY, "Required. The dimensions of the real-space output volume, including any padding (must be handled by the calling application). Assumed that apix x/y/z identical.");
                        d.put("weight", EMObject::FLOAT, "Optional. A temporary value set prior to slice insertion, indicative of the inserted slice's weight. Default sis 1.");
                        d.put("sym", EMObject::STRING, "Optional. The symmetry to impose on the final reconstruction. Default is c1");
                        d.put("verbose", EMObject::INT, "Optional. Toggles writing diagnostic information to standard out. Default is 0.");
                        return d;
                }

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

◆ insert_slice()

int RealMedianReconstructor::insert_slice	(	const EMData *const	slice,
		const Transform &	euler,
		const float	weight
	)

virtual

Insert an image slice to the reconstructor.

To insert multiple image slices, call this function multiple times.

Parameters

slice	Image slice.
euler	Euler angle of this image slice.
weight	A weighting factor for this slice, generally the number of particles in a class-average. May be ignored by some reconstructors

Returns: 0 if OK. 1 if error.

Reimplemented from EMAN::Reconstructor.

Definition at line 2550 of file reconstructor.cpp.

{
        if (!input) {
                LOGERR("try to insert NULL slice");
                return 1;
        }
 
        if (input->get_xsize() != input->get_ysize() || input->get_xsize() != nx) {
                LOGERR("tried to insert image that was not correction dimensions");
                return 1;
        }
        
        slices.push_back(input->copy());
        xforms.push_back(t);
 
        return 0;
}

References LOGERR, EMAN::ReconstructorVolumeData::nx, slices, and xforms.

◆ NEW()

static Reconstructor * EMAN::RealMedianReconstructor::NEW ( )

inlinestatic

Definition at line 979 of file reconstructor.h.

                {
                        return new RealMedianReconstructor();
                }

References RealMedianReconstructor().

◆ operator=()

RealMedianReconstructor & EMAN::RealMedianReconstructor::operator= ( const RealMedianReconstructor & )

private

◆ preprocess_slice()

EMData * RealMedianReconstructor::preprocess_slice	(	const EMData *const	slice,
		const Transform &	t = `Transform()`
	)

virtual

While you can just insert unprocessed slices, if you call preprocess_slice yourself, and insert the returned slice instead, repeatedly, it can save a fair bit of computation.

The default operation just returns a copy of the image, as the preprocessing is reconstructor-specific.

Returns: the processed slice

Parameters

slice	the slice to be prepocessed
t	transform

Exceptions

InvalidValueException when the specified padding value is less than the size of the images

Reimplemented from EMAN::Reconstructor.

Definition at line 2524 of file reconstructor.cpp.

{
 
//      EMData* return_slice = slice->process("normalize.edgemean");
//      return_slice->process_inplace("filter.linearfourier");
 
        EMData* return_slice;
 
//      return_slice = slice->process("filter.linearfourier");
        return_slice = slice->copy();
        
//      Transform tmp(t);
//      tmp.set_rotation(Dict("type","eman")); // resets the rotation to 0 implicitly
//      Vec2f trans = tmp.get_trans_2d();
//      float scale = tmp.get_scale();
//      bool mirror = tmp.get_mirror();
//      if (trans[0] != 0 || trans[1] != 0 || scale != 1.0 ) {
//              return_slice->transform(tmp);
//      } 
//      if ( mirror == true ) {
//              return_slice->process_inplace("xform.flip",Dict("axis","x"));
//      }
 
        return return_slice;
}

◆ setup()

void RealMedianReconstructor::setup ( )

virtual

Initialize the reconstructor.

Implements EMAN::Reconstructor.

Definition at line 2508 of file reconstructor.cpp.

{
        vector<int> size=params["size"];
        nx = size[0];
        ny = size[1];
        nz = size[2];
        image = new EMData(nx,ny,nz);
        image->to_zero();
        
        if (!slices.empty()) {
                for (std::vector<EMData *>::iterator it = slices.begin() ; it != slices.end(); ++it) delete *it;
                slices.clear();
        }
        xforms.clear();
}

References EMAN::ReconstructorVolumeData::image, EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::nz, EMAN::FactoryBase::params, slices, and xforms.

Member Data Documentation

◆ NAME

const string RealMedianReconstructor::NAME = "real_median"

static

Definition at line 995 of file reconstructor.h.

Referenced by get_name().

◆ slices

vector<EMData *> EMAN::RealMedianReconstructor::slices

private

Definition at line 1003 of file reconstructor.h.

Referenced by finish(), insert_slice(), and setup().

◆ xforms

vector<Transform> EMAN::RealMedianReconstructor::xforms

private

Definition at line 1004 of file reconstructor.h.

Referenced by finish(), insert_slice(), and setup().

The documentation for this class was generated from the following files:

Public Member Functions

Static Public Member Functions

Static Public Attributes

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ RealMedianReconstructor() [1/2]

◆ ~RealMedianReconstructor()

◆ RealMedianReconstructor() [2/2]

Member Function Documentation

◆ determine_slice_agreement()

◆ finish()

◆ get_desc()

◆ get_name()

◆ get_param_types()

◆ insert_slice()

◆ NEW()

◆ operator=()

◆ preprocess_slice()

◆ setup()

Member Data Documentation

◆ NAME

◆ slices

◆ xforms