emdata.h

Go to the documentation of this file.

/*
 * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
 * Copyright (c) 2000-2006 Baylor College of Medicine
 *
 * This software is issued under a joint BSD/GNU license. You may use the
 * source code in this file under either license. However, note that the
 * complete EMAN2 and SPARX software packages have some GPL dependencies,
 * so you are responsible for compliance with the licenses of these packages
 * if you opt to use BSD licensing. The warranty disclaimer below holds
 * in either instance.
 *
 * This complete copyright notice must be included in any revised version of the
 * source code. Additional authorship citations may be added, but existing
 * author citations must be preserved.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 *
 * */
 
#ifndef eman__emdata_h__
#define eman__emdata_h__ 1
 
#include <cfloat>
#include <complex>
#include <fstream>
 
#include "sparx/fundamentals.h"
#include "emutil.h"
#include "util.h"
#include "sparx/emarray.h"
#include "geometry.h"
#include "transform.h"
#include "vec3.h"
#ifdef EMAN2_USING_CUDA
#include <cuda_runtime_api.h>
#include "cuda/cuda_util.h"
#endif // EMAN2_USING_CUDA
using std::string;
using std::vector;
using std::map;
//using std::complex;   //comment this out for conflict with ACML
using std::ostream;
 
#include <utility>
using std::pair;
 
namespace EMAN
{
        class ImageIO;
        class Ctf;
        class XYData;
        class Transform;
        class GLUtil;
        class EMBytes: public std::string {};
 
        typedef boost::multi_array_ref<float, 2> MArray2D;
        typedef boost::multi_array_ref<float, 3> MArray3D;
        typedef boost::multi_array_ref<std::complex<float>, 2> MCArray2D;
        typedef boost::multi_array_ref<std::complex<float>, 3> MCArray3D;
        typedef boost::multi_array<int, 2> MIArray2D;
        typedef boost::multi_array<int, 3> MIArray3D;
 
        class EMData
        {
                friend class GLUtil;
 
                #include "emdata_io.h"
 
                #include "emdata_metadata.h"
 
                #include "emdata_modular.h"
 
                #include "emdata_transform.h"
 
                #include "emdata_core.h"
 
                #include "sparx/emdata_sparx.h"
 
                #include "sphire/emdata_sphire.h"
 
#ifdef EMAN2_USING_CUDA
                #include "emdata_cuda.h"
#endif // EMAN2_USING_CUDA
 
        public:
                EMData();
                ~ EMData();
 
                explicit EMData(const string& filename, int image_index=0);
 
                EMData(int nx, int ny, int nz=1, bool is_real=true);
 
                EMData(float* data, const int nx, const int ny, const int nz, const Dict& attr_dict = Dict());
                
                EMData(float* data, float* cudadata, const int nx, const int ny, const int nz, const Dict& attr_dict = Dict());
                //I do not wish to use a default dict if none is provided. Copying Dicts arround is really expensive in CUDA.
 
                EMData(const EMData& that);
 
                EMData& operator=(const EMData& that);
 
 
                EMData *get_clip(const Region & area, const float fill = 0) const;
 
                void clip_inplace(const Region & area,const float& fill_value=0);
 
                EMData *get_top_half() const;
 
 
                EMData *get_rotated_clip(const Transform & xform, const IntSize &size, float scale=1.0);
 
                EMData* window_center(int l);
 
 
                float *setup4slice(bool redo = true);
 
 
                void scale(float scale_factor);
 
 
                void translate(float dx, float dy, float dz);
 
 
                void translate(const Vec3f &translation);
 
 
                void translate(int dx, int dy, int dz);
 
 
                void translate(const Vec3i &translation);
 
 
                void rotate(const Transform & t);
 
                float max_3D_pixel_error(const Transform &t1, const Transform &t2, float r);
 
                void rotate(float az, float alt, float phi);
 
 
//              void rotate_translate(const Transform3D & t);
 
                inline void transform(const Transform& t) {
                        ENTERFUNC;
                        process_inplace("xform",Dict("transform",(Transform*)(&t)));
                        //update(); no need, process_inplace did it
                        EXITFUNC;
                }
 
                inline void rotate_translate(const Transform & t) {
                        cout << "Deprecation warning. Please consider using EMData::transform() instead " << endl;
                        transform(t); }
 
                void rotate_translate(float az, float alt, float phi, float dx, float dy, float dz);
 
 
                void rotate_translate(float az, float alt, float phi, float dx, float dy,
                                                          float dz, float pdx, float pdy, float pdz);
 
 
                void rotate_x(int dx);
 
 
                inline void rotate_180() {
                        ENTERFUNC;
                        process_inplace("math.rotate.180",Dict());
                        EXITFUNC;
                }
 
 
                double dot_rotate_translate(EMData * with, float dx, float dy, float da,const bool mirror=false);
 
 
                EMData *little_big_dot(EMData * little_img, bool do_sigma = false);
 
 
                EMData *do_radon();
 
 
                EMData *calc_ccf(EMData * with = 0, fp_flag fpflag = CIRCULANT, bool center=false);
 
                EMData *calc_ccf_masked(EMData *with,EMData *withsq=0,EMData *mask=0);
                
                void zero_corner_circulant(const int radius = 0);
 
                EMData *calc_ccfx( EMData * const with, int y0 = 0, int y1 = -1, bool nosum = false, bool flip = false,bool usez=false);
 
 
                EMData *make_rotational_footprint(bool unwrap = true);
                EMData *make_rotational_footprint_e1(bool unwrap = true);
                EMData *make_rotational_footprint_cmc(bool unwrap = true);
 
                EMData *make_footprint(int type=0);
 
 
                EMData *calc_mutual_correlation(EMData * with, bool tocorner = false, EMData * filter = 0);
 
 
                EMData *unwrap(int r1 = -1, int r2 = -1, int xs = -1, int dx = 0,
                                                           int dy = 0, bool do360 = false, bool weight_radial=true) const;
 
                EMData * unwrap_largerR(int r1,int r2,int xs, float rmax_f);
 
                EMData *oneDfftPolar(int size, float rmax, float MAXR);
 
                
                void apply_radial_func(float x0, float dx, vector < float >array, bool interp = true);
 
 
                vector < float >calc_radial_dist(int n, float x0, float dx,int inten);
 
 
                vector < float >calc_radial_dist(int n, float x0, float dx, int nwedge, float offset,bool inten);
 
 
                void cconj();
 
 
                void add_incoherent(EMData * obj);
 
 
                vector <float> calc_hist(int hist_size = 128, float hist_min = 0, float hist_max = 0, const float& brt = 0.0f, const float& cont = 1.0f);
 
 
                vector<float> calc_az_dist(int n, float a0, float da, float rmin,
                                                                   float rmax);
 
#if 0
                void calc_rcf(EMData * with, vector < float >&sum_array);
#endif
                float calc_dist(EMData * second_img, int y_index = 0) const;
 
                EMData *calc_flcf(EMData * with);
 
                EMData *calc_fast_sigma_image( EMData* mask);
 
                EMData *convolute(EMData * with);
 
#if 0
                void create_ctf_map(CtfMapType type, XYData * sf = 0);
#endif
 
 
                void common_lines(EMData * image1, EMData * image2, int mode = 0,
                                                  int steps = 180, bool horizontal = false);
 
                void common_lines_real(EMData * image1, EMData * image2,
                                                           int steps = 180, bool horizontal = false);
 
                void cut_slice(const EMData * const map, const Transform& tr, bool interpolate = true);
 
                void uncut_slice(EMData * const map, const Transform& tr) const;
 
                EMData *extract_box(const Transform& cs, const Region& r);
                
                int getResolution() const {
                        int resolution = 0;
                        int num = 1;
                        while(num < get_xsize()) {
                                resolution++;
                                num = 1 << resolution;
                        }
 
                        return resolution;
                }
                
                void debug_print_parms()
                {
                        std::cout << "Printing EMData params" << std::endl;
                        for ( Dict::const_iterator it = attr_dict.begin(); it != attr_dict.end(); ++it )
                        {
                                std::cout << (it->first) << " " << (it->second).to_str() << std::endl;
                        }
                        std::cout << "Done printing EMData params" << std::endl;
                }
 
                void set_xyz_origin(float origin_x, float origin_y, float origin_z);
                
                EMData* compute_missingwedge(float wedgeangle, float start = 0.05, float stop = 0.5);
                
                static int totalalloc;
        private:
                enum EMDataFlags {
//                      EMDATA_COMPLEX = 1 << 1,
//                      EMDATA_RI = 1 << 2,            // real/imaginary or amp/phase
                        EMDATA_BUSY = 1 << 3,      // someone is modifying data, NO LONGER USED
                        EMDATA_HASCTFF = 1 << 4,   // has CTF info in the image file
                        EMDATA_NEEDUPD = 1 << 5,   // needs a real update
//                      EMDATA_COMPLEXX = 1 << 6,  // 1D fft's in X
                        EMDATA_FLIP = 1 << 7,      // is the image flipped
                        EMDATA_PAD = 1 << 8,       // is the image fft padded
                        EMDATA_FFTODD = 1 << 9,    // is the (real-space) nx odd
                        EMDATA_SHUFFLE = 1 << 10,  // fft been shuffled? (so O is centered) PRB
                        EMDATA_FH = 1 << 11,        // is the complex image a FH image
                        EMDATA_CPU_NEEDS_UPDATE = 1 << 12, // CUDA related: is the CPU version of the image out out data
                        EMDATA_GPU_NEEDS_UPDATE = 1 << 13, // CUDA related: is the GPU version of the image out out data
                        EMDATA_GPU_RO_NEEDS_UPDATE = 1 << 14 // // CUDA related: is the GPU RO version of the image out out data
                };
 
                void update_stat() const;
                void save_byteorder_to_dict(ImageIO * imageio);
 
        private:
                mutable Dict attr_dict;
                mutable float *rdata;
                float *supp;
 
                //Ctf *ctf;
 
                mutable int flags;
                // Incremented every time the image changes
                int changecount;
                int nx, ny, nz, nxy;
                size_t nxyz;
                int xoff, yoff, zoff;
 
                Vec3f all_translation;
//              Vec3f all_rotation;    /** rotation (az, alt, phi) from the original locaton*/
 
                string path;
                int pathnum;
 
                mutable EMData* rot_fp;
 
#ifdef FFT_CACHING
                mutable EMData *fftcache;
#endif
 
                // Clip inplace variables is a local class used from convenience in EMData::clip_inplace
                // Added by d.woolford
                class ClipInplaceVariables
                {
                        public:
                                ClipInplaceVariables(const int p_nx, const int p_ny, const int p_nz, const int n_nx, const int n_ny, const int n_nz,const int xtrans, const int ytrans, const int ztrans) :
                                        prv_nx(p_nx), prv_ny(p_ny), prv_nz(p_nz), new_nx(n_nx), new_ny(n_ny), new_nz(n_nz), xshift(xtrans), yshift(ytrans), zshift(ztrans),
                                 x_iter(prv_nx), y_iter(prv_ny), z_iter(prv_nz), new_z_top(0), new_z_bottom(0),  new_y_back(0), new_y_front(0),new_x_left(0), new_x_right(0),
                                prv_z_top(0), prv_z_bottom(0), prv_y_back(0), prv_y_front(0), prv_x_left(0), prv_x_right(0)
                        {
                                if ( xtrans > 0 ) x_iter -= xtrans;
                                if ( x_iter < 0 ) x_iter = 0;
                                if ( ytrans > 0 ) y_iter -= ytrans;
                                if ( y_iter < 0 ) y_iter = 0;
                                if ( ztrans > 0 ) z_iter -= ztrans;
                                if ( z_iter < 0 ) z_iter = 0;
 
                                // Get the depth in the new volume where slices are inserted
                                // if this value is zero it means that the last z-slice in the new
                                // volume contains image data
                                if ( (new_nz + ztrans) > prv_nz ) new_z_top = new_nz + ztrans - prv_nz;
                                if ( (new_ny + ytrans) > prv_ny ) new_y_back = new_ny + ytrans - prv_ny;
                                if ( (new_nx + xtrans) > prv_nx ) new_x_right = new_nx + xtrans - prv_nx;
 
                                if ( (new_nz + ztrans) < prv_nz )
                                {
                                        prv_z_top = prv_nz - new_nz - ztrans;
                                        z_iter -= prv_z_top;
                                }
                                if ( (new_ny + ytrans) < prv_ny )
                                {
                                        prv_y_back = prv_ny - new_ny - ytrans;
                                        y_iter -= prv_y_back;
                                }
                                if ( (new_nx + xtrans) < prv_nx )
                                {
                                        prv_x_right = prv_nx - new_nx - xtrans;
                                        x_iter -= prv_x_right;
                                }
 
                                if ( xtrans > 0 ) prv_x_left = xtrans;
                                if ( ytrans > 0 ) prv_y_front = ytrans;
                                if ( ztrans > 0 ) prv_z_bottom = ztrans;
 
                                if ( xtrans < 0 ) new_x_left = -xtrans;
                                if ( ytrans < 0 ) new_y_front = -ytrans;
                                if ( ztrans < 0 ) new_z_bottom = -ztrans;
 
                        }
                        ~ClipInplaceVariables() {}
 
                        int prv_nx, prv_ny, prv_nz, new_nx, new_ny, new_nz;
                        int xshift, yshift, zshift;
                        int x_iter, y_iter, z_iter;
                        int new_z_top, new_z_bottom, new_y_back, new_y_front, new_x_left, new_x_right;
                        int prv_z_top, prv_z_bottom,  prv_y_back, prv_y_front, prv_x_left, prv_x_right;
                };
 
 
 
                //              /**  Do the Fourier Harmonic Transform  PRB
                //               * Takes a real image, returns the FH
                //               * Sets the EMDATA_FH switch to indicate that it is an FH image
                //               * @exception ImageFormatException If the image is not a square real odd image.
                //               * @return the FH image.
                //               */
                //              EMData* do_FH();
 
                //              /**   Do the Inverse Fourier Harmonic Transform   PRB
                //               * Takes an FH image, returns a square  complex image with odd sides
                //               * @exception ImageFormatException If the image is not the FH of something
                //               * @return a square complex image with odd sides
                //               */
                //              EMData* do_FH2F();
 
 
 
                //              /** Caclulates normalization and phase residual for a slice in
                //               * an already existing volume. phase residual is calculated
                //               * using only the inner 1/2 of the fourier sphere. Both the
                //               * slice image and this image must be in complex image format.
                //               *
                //               * @param slice An slice image to be normalized.
                //               * @param orient Orientation of the slice.
                //               * @exception ImageFormatException If the images are not complex.
                //               * @exception ImageDimensionException If the image is 3D.
                //               * @return A float number pair (result, phase-residual).
                //               */
                //              FloatPoint normalize_slice(EMData * slice, const Transform3D & orient);
 
                //              /** Caclulates normalization and phase residual for a slice in
                //               * an already existing volume. phase residual is calculated
                //               * using only the inner 1/2 of the fourier sphere. Both the
                //               * slice image and this image must be in complex image format.
                //               *
                //               * @param slice An slice image to be normalized.
                //               * @param alt Orientation euler angle alt (in EMAN convention).
                //               * @param az  Orientation euler angle az  (in EMAN convention).
                //               * @param phi Orientation euler angle phi (in EMAN convention).
                //               * @exception ImageFormatException If the images are not complex.
                //               * @exception ImageDimensionException If the image is 3D.
                //               * @return A float number pair (result, phase-residual).
                //               */
                //              FloatPoint normalize_slice(EMData * slice, float az, float alt, float phi);
 
                //              /** Get the normalization and phase residual values
                //               * Used for normalizaton and error measurement when 2D slices are inserted into a 3D volume of Fourier pixels
                //               * Originally added for use by the FourierReconstructor object
                //               * @return the normalization const (pair.first) and the phase residual (pair.second)
                //               * @param slice -the slice to be inserted into the 3D volume
                //               * @param euler - the euler angle orientation of the slice
                //               * @exception ImageDimensionException If this image is not 3D.ImageFormatException
                //               * @exception ImageFormatException If this image is not complex
                //               * @exception ImageFormatException If the slice not complex
                //               */
                //              pair<float, float> get_normalization_and_phaseres( const EMData* const slice, const Transform3D& euler );
 
                //              /** cut a 2D slice out of a this 3D image and return it
                //               * An alternative to cut_slice
                //               * @param tr orientation of the slice as encapsulated in a Transform object
                //               * @exception ImageDimensionException If this image is not 3D.
                //               * @exception ImageFormatException If this image is complex
                //               * @author David Woolford (adapted from an original version by Steve Ludtke)
                //               * @date Feb 2009
                //               */
                //              EMData* get_cut_slice(const Transform& tr);
 
        };
 
 
        EMData * operator+(const EMData & em, float n);
        EMData * operator-(const EMData & em, float n);
        EMData * operator*(const EMData & em, float n);
        EMData * operator/(const EMData & em, float n);
 
        EMData * operator+(float n, const EMData & em);
        EMData * operator-(float n, const EMData & em);
        EMData * operator*(float n, const EMData & em);
        EMData * operator/(float n, const EMData & em);
 
        EMData * rsub(const EMData & em, float n);
        EMData * rdiv(const EMData & em, float n);
 
        EMData * operator+(const EMData & a, const EMData & b);
        EMData * operator-(const EMData & a, const EMData & b);
        EMData * operator*(const EMData & a, const EMData & b);
        EMData * operator/(const EMData & a, const EMData & b);
 
 
/*   Next  is Modified by PRB      Transform3D::EMAN,
        inline Transform3D EMData::get_transform() const
        {
                return Transform3D((float)attr_dict["euler_alt"],
                                   (float)attr_dict["euler_az"],
                                   (float)attr_dict["euler_phi"]);
        }
*/
 
 
}
 
 
#endif