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

#ifdef _WIN32
        #pragma warning(disable:4819)
#endif  //_WIN32

#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"
#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 Transform3D;
        class Transform;
        class GLUtil;

        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"
#ifdef EMAN2_USING_CUDA

                #include "emdata_cuda.h"
#endif // EMAN2_USING_CUDA

        static int totalalloc;
        public:
                enum FFTPLACE { FFT_OUT_OF_PLACE, FFT_IN_PLACE };
                enum WINDOWPLACE { WINDOW_OUT_OF_PLACE, WINDOW_IN_PLACE };

                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(const EMData& that);

                EMData& operator=(const EMData& that);

//              /**  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);


                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;

//              /** 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 );


                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 Transform3D & t);


                void rotate(float az, float alt, float phi);


                void rotate_translate(const Transform3D & t);

                inline void transform(const Transform& t) {
                        ENTERFUNC;
                        process_inplace("math.transform",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, fp_flag fpflag = CIRCULANT, bool center=false);

                void zero_corner_circulant(const int radius = 0);

                EMData *calc_ccfx( EMData * const with, int y0 = 0, int y1 = -1, bool nosum = 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;


                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,bool inten);


                vector < float >calc_radial_dist(int n, float x0, float dx, int nwedge, 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);

//              EMData* get_cut_slice(const Transform& tr);
#ifdef EMAN2_USING_CUDA

                EMData* cut_slice_cuda(const Transform& tr);
#endif // EMAN2_USING_CUDA

                void uncut_slice(EMData * const map, const Transform& tr) const;

                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);
        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;
                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;

                // 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;
                };
        };


        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

/* vim: set ts=4 noet nospell: */