EMAN2
emdata.h
Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  *
00034  * */
00035 
00036 #ifndef eman__emdata_h__
00037 #define eman__emdata_h__ 1
00038 
00039 #ifdef _WIN32
00040         #pragma warning(disable:4819)
00041 #endif  //_WIN32
00042 
00043 #include <cfloat>
00044 #include <complex>
00045 #include <fstream>
00046 
00047 #include "sparx/fundamentals.h"
00048 #include "emutil.h"
00049 #include "util.h"
00050 #include "sparx/emarray.h"
00051 #include "geometry.h"
00052 #include "transform.h"
00053 #include "vec3.h"
00054 #ifdef EMAN2_USING_CUDA
00055 #include <cuda_runtime_api.h>
00056 #include "cuda/cuda_util.h"
00057 #endif // EMAN2_USING_CUDA
00058 using std::string;
00059 using std::vector;
00060 using std::map;
00061 //using std::complex;   //comment this out for conflict with ACML
00062 using std::ostream;
00063 
00064 #include <utility>
00065 using std::pair;
00066 
00067 namespace EMAN
00068 {
00069         class ImageIO;
00070         class Ctf;
00071         class XYData;
00072         class Transform;
00073         class GLUtil;
00074 
00075         typedef boost::multi_array_ref<float, 2> MArray2D;
00076         typedef boost::multi_array_ref<float, 3> MArray3D;
00077         typedef boost::multi_array_ref<std::complex<float>, 2> MCArray2D;
00078         typedef boost::multi_array_ref<std::complex<float>, 3> MCArray3D;
00079         typedef boost::multi_array<int, 2> MIArray2D;
00080         typedef boost::multi_array<int, 3> MIArray3D;
00081 
00088         class EMData
00089         {
00090                 friend class GLUtil;
00091 
00093                 #include "emdata_io.h"
00094 
00096                 #include "emdata_metadata.h"
00097 
00099                 #include "emdata_modular.h"
00100 
00102                 #include "emdata_transform.h"
00103 
00105                 #include "emdata_core.h"
00106 
00108                 #include "sparx/emdata_sparx.h"
00109 #ifdef EMAN2_USING_CUDA
00110 
00111                 #include "emdata_cuda.h"
00112 #endif // EMAN2_USING_CUDA
00113 
00114         public:
00115                 enum FFTPLACE { FFT_OUT_OF_PLACE, FFT_IN_PLACE };
00116                 enum WINDOWPLACE { WINDOW_OUT_OF_PLACE, WINDOW_IN_PLACE };
00117 
00119                 EMData();
00120                 ~ EMData();
00121 
00125                 explicit EMData(const string& filename, int image_index=0);
00126 
00134                 EMData(int nx, int ny, int nz=1, bool is_real=true);
00135 
00145                 EMData(float* data, const int nx, const int ny, const int nz, const Dict& attr_dict = Dict());
00146                 
00157                 EMData(float* data, float* cudadata, const int nx, const int ny, const int nz, const Dict& attr_dict = Dict());
00158                 //I do not wish to use a default dict if none is provided. Copying Dicts arround is really expensive in CUDA.
00159 
00165                 EMData(const EMData& that);
00166 
00172                 EMData& operator=(const EMData& that);
00173 
00174 
00182                 EMData *get_clip(const Region & area, const float fill = 0) const;
00183 
00190                 void clip_inplace(const Region & area,const float& fill_value=0);
00191 
00196                 EMData *get_top_half() const;
00197 
00198 
00207                 EMData *get_rotated_clip(const Transform & xform, const IntSize &size, float scale=1.0);
00208 
00219                 EMData* window_center(int l);
00220 
00221 
00236                 float *setup4slice(bool redo = true);
00237 
00238 
00242                 void scale(float scale_factor);
00243 
00244 
00250                 void translate(float dx, float dy, float dz);
00251 
00252 
00256                 void translate(const Vec3f &translation);
00257 
00258 
00265                 void translate(int dx, int dy, int dz);
00266 
00267 
00272                 void translate(const Vec3i &translation);
00273 
00274 
00279                 void rotate(const Transform & t);
00280 
00281                 float max_3D_pixel_error(const Transform &t1, const Transform &t2, float r);
00282 
00289                 void rotate(float az, float alt, float phi);
00290 
00291 
00296 //              void rotate_translate(const Transform3D & t);
00297 
00301                 inline void transform(const Transform& t) {
00302                         ENTERFUNC;
00303                         process_inplace("xform",Dict("transform",(Transform*)(&t)));
00304                         //update(); no need, process_inplace did it
00305                         EXITFUNC;
00306                 }
00307 
00312                 inline void rotate_translate(const Transform & t) {
00313                         cout << "Deprecation warning. Please consider using EMData::transform() instead " << endl;
00314                         transform(t); }
00315 
00325                 void rotate_translate(float az, float alt, float phi, float dx, float dy, float dz);
00326 
00327 
00340                 void rotate_translate(float az, float alt, float phi, float dx, float dy,
00341                                                           float dz, float pdx, float pdy, float pdz);
00342 
00343 
00349                 void rotate_x(int dx);
00350 
00351 
00356                 inline void rotate_180() {
00357                         ENTERFUNC;
00358                         process_inplace("math.rotate.180",Dict());
00359                         EXITFUNC;
00360                 }
00361 
00362 
00376                 double dot_rotate_translate(EMData * with, float dx, float dy, float da,const bool mirror=false);
00377 
00378 
00389                 EMData *little_big_dot(EMData * little_img, bool do_sigma = false);
00390 
00391 
00404                 EMData *do_radon();
00405 
00406 
00422                 EMData *calc_ccf(EMData * with = 0, fp_flag fpflag = CIRCULANT, bool center=false);
00423 
00436                 void zero_corner_circulant(const int radius = 0);
00437 
00457                 EMData *calc_ccfx( EMData * const with, int y0 = 0, int y1 = -1, bool nosum = false, bool flip = false,bool usez=false);
00458 
00459 
00468                 EMData *make_rotational_footprint(bool unwrap = true);
00469                 EMData *make_rotational_footprint_e1(bool unwrap = true);
00470                 EMData *make_rotational_footprint_cmc(bool unwrap = true);
00471 
00488                 EMData *make_footprint(int type=0);
00489 
00490 
00503                 EMData *calc_mutual_correlation(EMData * with, bool tocorner = false, EMData * filter = 0);
00504 
00505 
00521                 EMData *unwrap(int r1 = -1, int r2 = -1, int xs = -1, int dx = 0,
00522                                                            int dy = 0, bool do360 = false, bool weight_radial=true) const;
00523 
00524                 EMData * unwrap_largerR(int r1,int r2,int xs, float rmax_f);
00525 
00526                 EMData *oneDfftPolar(int size, float rmax, float MAXR);
00527 
00528                 
00536                 void apply_radial_func(float x0, float dx, vector < float >array, bool interp = true);
00537 
00538 
00550                 vector < float >calc_radial_dist(int n, float x0, float dx,bool inten);
00551 
00552 
00565                 vector < float >calc_radial_dist(int n, float x0, float dx, int nwedge, float offset,bool inten);
00566 
00567 
00571                 void cconj();
00572 
00573 
00581                 void add_incoherent(EMData * obj);
00582 
00583 
00595                 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);
00596 
00597 
00609                 vector<float> calc_az_dist(int n, float a0, float da, float rmin,
00610                                                                    float rmax);
00611 
00612 #if 0
00613                 void calc_rcf(EMData * with, vector < float >&sum_array);
00614 #endif
00615 
00628                 float calc_dist(EMData * second_img, int y_index = 0) const;
00629 
00646                 EMData *calc_flcf(EMData * with);
00647 
00666                 EMData *calc_fast_sigma_image( EMData* mask);
00667 
00673                 EMData *convolute(EMData * with);
00674 
00675 #if 0
00676                 void create_ctf_map(CtfMapType type, XYData * sf = 0);
00677 #endif
00678 
00679 
00701                 void common_lines(EMData * image1, EMData * image2, int mode = 0,
00702                                                   int steps = 180, bool horizontal = false);
00703 
00714                 void common_lines_real(EMData * image1, EMData * image2,
00715                                                            int steps = 180, bool horizontal = false);
00716 
00730                 void cut_slice(const EMData * const map, const Transform& tr, bool interpolate = true);
00731 
00746                 void uncut_slice(EMData * const map, const Transform& tr) const;
00747 
00756                 EMData *extract_box(const Transform& cs, const Region& r);
00757                 
00761                 int getResolution() const {
00762                         int resolution = 0;
00763                         int num = 1;
00764                         while(num < get_xsize()) {
00765                                 resolution++;
00766                                 num = 1 << resolution;
00767                         }
00768 
00769                         return resolution;
00770                 }
00771                 
00774                 void debug_print_parms()
00775                 {
00776                         std::cout << "Printing EMData params" << std::endl;
00777                         for ( Dict::const_iterator it = attr_dict.begin(); it != attr_dict.end(); ++it )
00778                         {
00779                                 std::cout << (it->first) << " " << (it->second).to_str() << std::endl;
00780                         }
00781                         std::cout << "Done printing EMData params" << std::endl;
00782                 }
00783 
00789                 void set_xyz_origin(float origin_x, float origin_y, float origin_z);
00790                 
00796                 EMData* compute_missingwedge(float wedgeangle, float start = 0.05, float stop = 0.5);
00797                 
00798                 static int totalalloc;
00799         private:
00806                 enum EMDataFlags {
00807 //                      EMDATA_COMPLEX = 1 << 1,
00808 //                      EMDATA_RI = 1 << 2,            // real/imaginary or amp/phase
00809                         EMDATA_BUSY = 1 << 3,      // someone is modifying data, NO LONGER USED
00810                         EMDATA_HASCTFF = 1 << 4,   // has CTF info in the image file
00811                         EMDATA_NEEDUPD = 1 << 5,   // needs a real update
00812 //                      EMDATA_COMPLEXX = 1 << 6,  // 1D fft's in X
00813                         EMDATA_FLIP = 1 << 7,      // is the image flipped
00814                         EMDATA_PAD = 1 << 8,       // is the image fft padded
00815                         EMDATA_FFTODD = 1 << 9,    // is the (real-space) nx odd
00816                         EMDATA_SHUFFLE = 1 << 10,  // fft been shuffled? (so O is centered) PRB
00817                         EMDATA_FH = 1 << 11,        // is the complex image a FH image
00818                         EMDATA_CPU_NEEDS_UPDATE = 1 << 12, // CUDA related: is the CPU version of the image out out data
00819                         EMDATA_GPU_NEEDS_UPDATE = 1 << 13, // CUDA related: is the GPU version of the image out out data
00820                         EMDATA_GPU_RO_NEEDS_UPDATE = 1 << 14 // // CUDA related: is the GPU RO version of the image out out data
00821                 };
00822 
00823                 void update_stat() const;
00824                 void save_byteorder_to_dict(ImageIO * imageio);
00825 
00826         private:
00828                 mutable Dict attr_dict;
00830                 mutable float *rdata;
00832                 float *supp;
00833 
00836                 //Ctf *ctf;
00837 
00839                 mutable int flags;
00840                 // Incremented every time the image changes
00841                 int changecount;
00843                 int nx, ny, nz, nxy;
00844                 size_t nxyz;
00846                 int xoff, yoff, zoff;
00847 
00849                 Vec3f all_translation;
00850 //              Vec3f all_rotation;    /** rotation (az, alt, phi) from the original locaton*/
00851 
00852                 string path;
00853                 int pathnum;
00854 
00856                 mutable EMData* rot_fp;
00857 
00858 #ifdef FFT_CACHING
00859                 mutable EMData *fftcache;
00860 #endif
00861 
00862                 // Clip inplace variables is a local class used from convenience in EMData::clip_inplace
00863                 // Added by d.woolford
00864                 class ClipInplaceVariables
00865                 {
00866                         public:
00867                                 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) :
00868                                         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),
00869                                  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),
00870                                 prv_z_top(0), prv_z_bottom(0), prv_y_back(0), prv_y_front(0), prv_x_left(0), prv_x_right(0)
00871                         {
00872                                 if ( xtrans > 0 ) x_iter -= xtrans;
00873                                 if ( x_iter < 0 ) x_iter = 0;
00874                                 if ( ytrans > 0 ) y_iter -= ytrans;
00875                                 if ( y_iter < 0 ) y_iter = 0;
00876                                 if ( ztrans > 0 ) z_iter -= ztrans;
00877                                 if ( z_iter < 0 ) z_iter = 0;
00878 
00879                                 // Get the depth in the new volume where slices are inserted
00880                                 // if this value is zero it means that the last z-slice in the new
00881                                 // volume contains image data
00882                                 if ( (new_nz + ztrans) > prv_nz ) new_z_top = new_nz + ztrans - prv_nz;
00883                                 if ( (new_ny + ytrans) > prv_ny ) new_y_back = new_ny + ytrans - prv_ny;
00884                                 if ( (new_nx + xtrans) > prv_nx ) new_x_right = new_nx + xtrans - prv_nx;
00885 
00886                                 if ( (new_nz + ztrans) < prv_nz )
00887                                 {
00888                                         prv_z_top = prv_nz - new_nz - ztrans;
00889                                         z_iter -= prv_z_top;
00890                                 }
00891                                 if ( (new_ny + ytrans) < prv_ny )
00892                                 {
00893                                         prv_y_back = prv_ny - new_ny - ytrans;
00894                                         y_iter -= prv_y_back;
00895                                 }
00896                                 if ( (new_nx + xtrans) < prv_nx )
00897                                 {
00898                                         prv_x_right = prv_nx - new_nx - xtrans;
00899                                         x_iter -= prv_x_right;
00900                                 }
00901 
00902                                 if ( xtrans > 0 ) prv_x_left = xtrans;
00903                                 if ( ytrans > 0 ) prv_y_front = ytrans;
00904                                 if ( ztrans > 0 ) prv_z_bottom = ztrans;
00905 
00906                                 if ( xtrans < 0 ) new_x_left = -xtrans;
00907                                 if ( ytrans < 0 ) new_y_front = -ytrans;
00908                                 if ( ztrans < 0 ) new_z_bottom = -ztrans;
00909 
00910                         }
00911                         ~ClipInplaceVariables() {}
00912 
00913                         int prv_nx, prv_ny, prv_nz, new_nx, new_ny, new_nz;
00914                         int xshift, yshift, zshift;
00915                         int x_iter, y_iter, z_iter;
00916                         int new_z_top, new_z_bottom, new_y_back, new_y_front, new_x_left, new_x_right;
00917                         int prv_z_top, prv_z_bottom,  prv_y_back, prv_y_front, prv_x_left, prv_x_right;
00918                 };
00919 
00920 
00921 
00922                 //              /**  Do the Fourier Harmonic Transform  PRB
00923                 //               * Takes a real image, returns the FH
00924                 //               * Sets the EMDATA_FH switch to indicate that it is an FH image
00925                 //               * @exception ImageFormatException If the image is not a square real odd image.
00926                 //               * @return the FH image.
00927                 //               */
00928                 //              EMData* do_FH();
00929 
00930                 //              /**   Do the Inverse Fourier Harmonic Transform   PRB
00931                 //               * Takes an FH image, returns a square  complex image with odd sides
00932                 //               * @exception ImageFormatException If the image is not the FH of something
00933                 //               * @return a square complex image with odd sides
00934                 //               */
00935                 //              EMData* do_FH2F();
00936 
00937 
00938 
00939                 //              /** Caclulates normalization and phase residual for a slice in
00940                 //               * an already existing volume. phase residual is calculated
00941                 //               * using only the inner 1/2 of the fourier sphere. Both the
00942                 //               * slice image and this image must be in complex image format.
00943                 //               *
00944                 //               * @param slice An slice image to be normalized.
00945                 //               * @param orient Orientation of the slice.
00946                 //               * @exception ImageFormatException If the images are not complex.
00947                 //               * @exception ImageDimensionException If the image is 3D.
00948                 //               * @return A float number pair (result, phase-residual).
00949                 //               */
00950                 //              FloatPoint normalize_slice(EMData * slice, const Transform3D & orient);
00951 
00952                 //              /** Caclulates normalization and phase residual for a slice in
00953                 //               * an already existing volume. phase residual is calculated
00954                 //               * using only the inner 1/2 of the fourier sphere. Both the
00955                 //               * slice image and this image must be in complex image format.
00956                 //               *
00957                 //               * @param slice An slice image to be normalized.
00958                 //               * @param alt Orientation euler angle alt (in EMAN convention).
00959                 //               * @param az  Orientation euler angle az  (in EMAN convention).
00960                 //               * @param phi Orientation euler angle phi (in EMAN convention).
00961                 //               * @exception ImageFormatException If the images are not complex.
00962                 //               * @exception ImageDimensionException If the image is 3D.
00963                 //               * @return A float number pair (result, phase-residual).
00964                 //               */
00965                 //              FloatPoint normalize_slice(EMData * slice, float az, float alt, float phi);
00966 
00967                 //              /** Get the normalization and phase residual values
00968                 //               * Used for normalizaton and error measurement when 2D slices are inserted into a 3D volume of Fourier pixels
00969                 //               * Originally added for use by the FourierReconstructor object
00970                 //               * @return the normalization const (pair.first) and the phase residual (pair.second)
00971                 //               * @param slice -the slice to be inserted into the 3D volume
00972                 //               * @param euler - the euler angle orientation of the slice
00973                 //               * @exception ImageDimensionException If this image is not 3D.ImageFormatException
00974                 //               * @exception ImageFormatException If this image is not complex
00975                 //               * @exception ImageFormatException If the slice not complex
00976                 //               */
00977                 //              pair<float, float> get_normalization_and_phaseres( const EMData* const slice, const Transform3D& euler );
00978 
00979                 //              /** cut a 2D slice out of a this 3D image and return it
00980                 //               * An alternative to cut_slice
00981                 //               * @param tr orientation of the slice as encapsulated in a Transform object
00982                 //               * @exception ImageDimensionException If this image is not 3D.
00983                 //               * @exception ImageFormatException If this image is complex
00984                 //               * @author David Woolford (adapted from an original version by Steve Ludtke)
00985                 //               * @date Feb 2009
00986                 //               */
00987                 //              EMData* get_cut_slice(const Transform& tr);
00988 
00989         };
00990 
00991 
00992         EMData * operator+(const EMData & em, float n);
00993         EMData * operator-(const EMData & em, float n);
00994         EMData * operator*(const EMData & em, float n);
00995         EMData * operator/(const EMData & em, float n);
00996 
00997         EMData * operator+(float n, const EMData & em);
00998         EMData * operator-(float n, const EMData & em);
00999         EMData * operator*(float n, const EMData & em);
01000         EMData * operator/(float n, const EMData & em);
01001 
01002         EMData * rsub(const EMData & em, float n);
01003         EMData * rdiv(const EMData & em, float n);
01004 
01005         EMData * operator+(const EMData & a, const EMData & b);
01006         EMData * operator-(const EMData & a, const EMData & b);
01007         EMData * operator*(const EMData & a, const EMData & b);
01008         EMData * operator/(const EMData & a, const EMData & b);
01009 
01010 
01011 /*   Next  is Modified by PRB      Transform3D::EMAN,
01012         inline Transform3D EMData::get_transform() const
01013         {
01014                 return Transform3D((float)attr_dict["euler_alt"],
01015                                    (float)attr_dict["euler_az"],
01016                                    (float)attr_dict["euler_phi"]);
01017         }
01018 */
01019 
01020 
01021 }
01022 
01023 
01024 #endif
01025 
01026 /* vim: set ts=4 noet nospell: */