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 #ifdef EMAN2_USING_CUDA
00054 #include <cuda_runtime_api.h>
00055 #include "cuda/cuda_util.h"
00056 #endif // EMAN2_USING_CUDA
00057 using std::string;
00058 using std::vector;
00059 using std::map;
00060 //using std::complex;   //comment this out for conflict with ACML
00061 using std::ostream;
00062 
00063 #include <utility>
00064 using std::pair;
00065 
00066 namespace EMAN
00067 {
00068         class ImageIO;
00069         class Ctf;
00070         class XYData;
00071         class Transform;
00072         class GLUtil;
00073 
00074         typedef boost::multi_array_ref<float, 2> MArray2D;
00075         typedef boost::multi_array_ref<float, 3> MArray3D;
00076         typedef boost::multi_array_ref<std::complex<float>, 2> MCArray2D;
00077         typedef boost::multi_array_ref<std::complex<float>, 3> MCArray3D;
00078         typedef boost::multi_array<int, 2> MIArray2D;
00079         typedef boost::multi_array<int, 3> MIArray3D;
00080 
00087         class EMData
00088         {
00089                 friend class GLUtil;
00090 
00092                 #include "emdata_io.h"
00093 
00095                 #include "emdata_metadata.h"
00096 
00098                 #include "emdata_modular.h"
00099 
00101                 #include "emdata_transform.h"
00102 
00104                 #include "emdata_core.h"
00105 
00107                 #include "sparx/emdata_sparx.h"
00108 #ifdef EMAN2_USING_CUDA
00109 
00110                 #include "emdata_cuda.h"
00111 #endif // EMAN2_USING_CUDA
00112 
00113         public:
00114                 enum FFTPLACE { FFT_OUT_OF_PLACE, FFT_IN_PLACE };
00115                 enum WINDOWPLACE { WINDOW_OUT_OF_PLACE, WINDOW_IN_PLACE };
00116 
00118                 EMData();
00119                 ~ EMData();
00120 
00124                 explicit EMData(const string& filename, int image_index=0);
00125 
00133                 EMData(int nx, int ny, int nz=1, bool is_real=true);
00134 
00144                 EMData(float* data, const int nx, const int ny, const int nz, const Dict& attr_dict = Dict());
00145                 
00156                 EMData(float* data, float* cudadata, const int nx, const int ny, const int nz, const Dict& attr_dict = Dict());
00157                 //I do not wish to use a default dict if none is provided. Copying Dicts arround is really expensive in CUDA.
00158 
00164                 EMData(const EMData& that);
00165 
00171                 EMData& operator=(const EMData& that);
00172 
00173 
00181                 EMData *get_clip(const Region & area, const float fill = 0) const;
00182 
00189                 void clip_inplace(const Region & area,const float& fill_value=0);
00190 
00195                 EMData *get_top_half() const;
00196 
00197 
00206                 EMData *get_rotated_clip(const Transform & xform, const IntSize &size, float scale=1.0);
00207 
00218                 EMData* window_center(int l);
00219 
00220 
00235                 float *setup4slice(bool redo = true);
00236 
00237 
00241                 void scale(float scale_factor);
00242 
00243 
00249                 void translate(float dx, float dy, float dz);
00250 
00251 
00255                 void translate(const Vec3f &translation);
00256 
00257 
00264                 void translate(int dx, int dy, int dz);
00265 
00266 
00271                 void translate(const Vec3i &translation);
00272 
00273 
00278                 void rotate(const Transform & t);
00279 
00280                 float max_3D_pixel_error(const Transform &t1, const Transform &t2, float r);
00281 
00288                 void rotate(float az, float alt, float phi);
00289 
00290 
00295 //              void rotate_translate(const Transform3D & t);
00296 
00300                 inline void transform(const Transform& t) {
00301                         ENTERFUNC;
00302                         process_inplace("xform",Dict("transform",(Transform*)(&t)));
00303                         //update(); no need, process_inplace did it
00304                         EXITFUNC;
00305                 }
00306 
00311                 inline void rotate_translate(const Transform & t) {
00312                         cout << "Deprecation warning. Please consider using EMData::transform() instead " << endl;
00313                         transform(t); }
00314 
00324                 void rotate_translate(float az, float alt, float phi, float dx, float dy, float dz);
00325 
00326 
00339                 void rotate_translate(float az, float alt, float phi, float dx, float dy,
00340                                                           float dz, float pdx, float pdy, float pdz);
00341 
00342 
00348                 void rotate_x(int dx);
00349 
00350 
00355                 inline void rotate_180() {
00356                         ENTERFUNC;
00357                         process_inplace("math.rotate.180",Dict());
00358                         EXITFUNC;
00359                 }
00360 
00361 
00375                 double dot_rotate_translate(EMData * with, float dx, float dy, float da,const bool mirror=false);
00376 
00377 
00388                 EMData *little_big_dot(EMData * little_img, bool do_sigma = false);
00389 
00390 
00403                 EMData *do_radon();
00404 
00405 
00421                 EMData *calc_ccf(EMData * with = 0, fp_flag fpflag = CIRCULANT, bool center=false);
00422 
00435                 void zero_corner_circulant(const int radius = 0);
00436 
00455                 EMData *calc_ccfx( EMData * const with, int y0 = 0, int y1 = -1, bool nosum = false, bool flip = false);
00456 
00457 
00466                 EMData *make_rotational_footprint(bool unwrap = true);
00467                 EMData *make_rotational_footprint_e1(bool unwrap = true);
00468                 EMData *make_rotational_footprint_cmc(bool unwrap = true);
00469 
00486                 EMData *make_footprint(int type=0);
00487 
00488 
00501                 EMData *calc_mutual_correlation(EMData * with, bool tocorner = false, EMData * filter = 0);
00502 
00503 
00519                 EMData *unwrap(int r1 = -1, int r2 = -1, int xs = -1, int dx = 0,
00520                                                            int dy = 0, bool do360 = false, bool weight_radial=true) const;
00521 
00522                 EMData * unwrap_largerR(int r1,int r2,int xs, float rmax_f);
00523 
00524                 EMData *oneDfftPolar(int size, float rmax, float MAXR);
00525 
00526                 
00534                 void apply_radial_func(float x0, float dx, vector < float >array, bool interp = true);
00535 
00536 
00548                 vector < float >calc_radial_dist(int n, float x0, float dx,bool inten);
00549 
00550 
00562                 vector < float >calc_radial_dist(int n, float x0, float dx, int nwedge, bool inten);
00563 
00564 
00568                 void cconj();
00569 
00570 
00578                 void add_incoherent(EMData * obj);
00579 
00580 
00592                 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);
00593 
00594 
00606                 vector<float> calc_az_dist(int n, float a0, float da, float rmin,
00607                                                                    float rmax);
00608 
00609 #if 0
00610                 void calc_rcf(EMData * with, vector < float >&sum_array);
00611 #endif
00612 
00625                 float calc_dist(EMData * second_img, int y_index = 0) const;
00626 
00643                 EMData *calc_flcf(EMData * with);
00644 
00663                 EMData *calc_fast_sigma_image( EMData* mask);
00664 
00670                 EMData *convolute(EMData * with);
00671 
00672 #if 0
00673                 void create_ctf_map(CtfMapType type, XYData * sf = 0);
00674 #endif
00675 
00676 
00698                 void common_lines(EMData * image1, EMData * image2, int mode = 0,
00699                                                   int steps = 180, bool horizontal = false);
00700 
00711                 void common_lines_real(EMData * image1, EMData * image2,
00712                                                            int steps = 180, bool horizontal = false);
00713 
00727                 void cut_slice(const EMData * const map, const Transform& tr, bool interpolate = true);
00728 
00743                 void uncut_slice(EMData * const map, const Transform& tr) const;
00744 
00748                 int getResolution() const {
00749                         int resolution = 0;
00750                         int num = 1;
00751                         while(num < get_xsize()) {
00752                                 resolution++;
00753                                 num = 1 << resolution;
00754                         }
00755 
00756                         return resolution;
00757                 }
00758 
00761                 void debug_print_parms()
00762                 {
00763                         std::cout << "Printing EMData params" << std::endl;
00764                         for ( Dict::const_iterator it = attr_dict.begin(); it != attr_dict.end(); ++it )
00765                         {
00766                                 std::cout << (it->first) << " " << (it->second).to_str() << std::endl;
00767                         }
00768                         std::cout << "Done printing EMData params" << std::endl;
00769                 }
00770 
00776                 void set_xyz_origin(float origin_x, float origin_y, float origin_z);
00777                 
00778                 static int totalalloc;
00779         private:
00786                 enum EMDataFlags {
00787 //                      EMDATA_COMPLEX = 1 << 1,
00788 //                      EMDATA_RI = 1 << 2,            // real/imaginary or amp/phase
00789                         EMDATA_BUSY = 1 << 3,      // someone is modifying data, NO LONGER USED
00790                         EMDATA_HASCTFF = 1 << 4,   // has CTF info in the image file
00791                         EMDATA_NEEDUPD = 1 << 5,   // needs a real update
00792 //                      EMDATA_COMPLEXX = 1 << 6,  // 1D fft's in X
00793                         EMDATA_FLIP = 1 << 7,      // is the image flipped
00794                         EMDATA_PAD = 1 << 8,       // is the image fft padded
00795                         EMDATA_FFTODD = 1 << 9,    // is the (real-space) nx odd
00796                         EMDATA_SHUFFLE = 1 << 10,  // fft been shuffled? (so O is centered) PRB
00797                         EMDATA_FH = 1 << 11,        // is the complex image a FH image
00798                         EMDATA_CPU_NEEDS_UPDATE = 1 << 12, // CUDA related: is the CPU version of the image out out data
00799                         EMDATA_GPU_NEEDS_UPDATE = 1 << 13, // CUDA related: is the GPU version of the image out out data
00800                         EMDATA_GPU_RO_NEEDS_UPDATE = 1 << 14 // // CUDA related: is the GPU RO version of the image out out data
00801                 };
00802 
00803                 void update_stat() const;
00804                 void save_byteorder_to_dict(ImageIO * imageio);
00805 
00806         private:
00808                 mutable Dict attr_dict;
00810                 mutable float *rdata;
00812                 float *supp;
00813 
00816                 //Ctf *ctf;
00817 
00819                 mutable int flags;
00820                 // Incremented every time the image changes
00821                 int changecount;
00823                 int nx, ny, nz, nxy;
00824                 size_t nxyz;
00826                 int xoff, yoff, zoff;
00827 
00829                 Vec3f all_translation;
00830 //              Vec3f all_rotation;    /** rotation (az, alt, phi) from the original locaton*/
00831 
00832                 string path;
00833                 int pathnum;
00834 
00836                 mutable EMData* rot_fp;
00837 
00838                 // Clip inplace variables is a local class used from convenience in EMData::clip_inplace
00839                 // Added by d.woolford
00840                 class ClipInplaceVariables
00841                 {
00842                         public:
00843                                 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) :
00844                                         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),
00845                                  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),
00846                                 prv_z_top(0), prv_z_bottom(0), prv_y_back(0), prv_y_front(0), prv_x_left(0), prv_x_right(0)
00847                         {
00848                                 if ( xtrans > 0 ) x_iter -= xtrans;
00849                                 if ( x_iter < 0 ) x_iter = 0;
00850                                 if ( ytrans > 0 ) y_iter -= ytrans;
00851                                 if ( y_iter < 0 ) y_iter = 0;
00852                                 if ( ztrans > 0 ) z_iter -= ztrans;
00853                                 if ( z_iter < 0 ) z_iter = 0;
00854 
00855                                 // Get the depth in the new volume where slices are inserted
00856                                 // if this value is zero it means that the last z-slice in the new
00857                                 // volume contains image data
00858                                 if ( (new_nz + ztrans) > prv_nz ) new_z_top = new_nz + ztrans - prv_nz;
00859                                 if ( (new_ny + ytrans) > prv_ny ) new_y_back = new_ny + ytrans - prv_ny;
00860                                 if ( (new_nx + xtrans) > prv_nx ) new_x_right = new_nx + xtrans - prv_nx;
00861 
00862                                 if ( (new_nz + ztrans) < prv_nz )
00863                                 {
00864                                         prv_z_top = prv_nz - new_nz - ztrans;
00865                                         z_iter -= prv_z_top;
00866                                 }
00867                                 if ( (new_ny + ytrans) < prv_ny )
00868                                 {
00869                                         prv_y_back = prv_ny - new_ny - ytrans;
00870                                         y_iter -= prv_y_back;
00871                                 }
00872                                 if ( (new_nx + xtrans) < prv_nx )
00873                                 {
00874                                         prv_x_right = prv_nx - new_nx - xtrans;
00875                                         x_iter -= prv_x_right;
00876                                 }
00877 
00878                                 if ( xtrans > 0 ) prv_x_left = xtrans;
00879                                 if ( ytrans > 0 ) prv_y_front = ytrans;
00880                                 if ( ztrans > 0 ) prv_z_bottom = ztrans;
00881 
00882                                 if ( xtrans < 0 ) new_x_left = -xtrans;
00883                                 if ( ytrans < 0 ) new_y_front = -ytrans;
00884                                 if ( ztrans < 0 ) new_z_bottom = -ztrans;
00885 
00886                         }
00887                         ~ClipInplaceVariables() {}
00888 
00889                         int prv_nx, prv_ny, prv_nz, new_nx, new_ny, new_nz;
00890                         int xshift, yshift, zshift;
00891                         int x_iter, y_iter, z_iter;
00892                         int new_z_top, new_z_bottom, new_y_back, new_y_front, new_x_left, new_x_right;
00893                         int prv_z_top, prv_z_bottom,  prv_y_back, prv_y_front, prv_x_left, prv_x_right;
00894                 };
00895 
00896 
00897 
00898                 //              /**  Do the Fourier Harmonic Transform  PRB
00899                 //               * Takes a real image, returns the FH
00900                 //               * Sets the EMDATA_FH switch to indicate that it is an FH image
00901                 //               * @exception ImageFormatException If the image is not a square real odd image.
00902                 //               * @return the FH image.
00903                 //               */
00904                 //              EMData* do_FH();
00905 
00906                 //              /**   Do the Inverse Fourier Harmonic Transform   PRB
00907                 //               * Takes an FH image, returns a square  complex image with odd sides
00908                 //               * @exception ImageFormatException If the image is not the FH of something
00909                 //               * @return a square complex image with odd sides
00910                 //               */
00911                 //              EMData* do_FH2F();
00912 
00913 
00914 
00915                 //              /** Caclulates normalization and phase residual for a slice in
00916                 //               * an already existing volume. phase residual is calculated
00917                 //               * using only the inner 1/2 of the fourier sphere. Both the
00918                 //               * slice image and this image must be in complex image format.
00919                 //               *
00920                 //               * @param slice An slice image to be normalized.
00921                 //               * @param orient Orientation of the slice.
00922                 //               * @exception ImageFormatException If the images are not complex.
00923                 //               * @exception ImageDimensionException If the image is 3D.
00924                 //               * @return A float number pair (result, phase-residual).
00925                 //               */
00926                 //              FloatPoint normalize_slice(EMData * slice, const Transform3D & orient);
00927 
00928                 //              /** Caclulates normalization and phase residual for a slice in
00929                 //               * an already existing volume. phase residual is calculated
00930                 //               * using only the inner 1/2 of the fourier sphere. Both the
00931                 //               * slice image and this image must be in complex image format.
00932                 //               *
00933                 //               * @param slice An slice image to be normalized.
00934                 //               * @param alt Orientation euler angle alt (in EMAN convention).
00935                 //               * @param az  Orientation euler angle az  (in EMAN convention).
00936                 //               * @param phi Orientation euler angle phi (in EMAN convention).
00937                 //               * @exception ImageFormatException If the images are not complex.
00938                 //               * @exception ImageDimensionException If the image is 3D.
00939                 //               * @return A float number pair (result, phase-residual).
00940                 //               */
00941                 //              FloatPoint normalize_slice(EMData * slice, float az, float alt, float phi);
00942 
00943                 //              /** Get the normalization and phase residual values
00944                 //               * Used for normalizaton and error measurement when 2D slices are inserted into a 3D volume of Fourier pixels
00945                 //               * Originally added for use by the FourierReconstructor object
00946                 //               * @return the normalization const (pair.first) and the phase residual (pair.second)
00947                 //               * @param slice -the slice to be inserted into the 3D volume
00948                 //               * @param euler - the euler angle orientation of the slice
00949                 //               * @exception ImageDimensionException If this image is not 3D.ImageFormatException
00950                 //               * @exception ImageFormatException If this image is not complex
00951                 //               * @exception ImageFormatException If the slice not complex
00952                 //               */
00953                 //              pair<float, float> get_normalization_and_phaseres( const EMData* const slice, const Transform3D& euler );
00954 
00955                 //              /** cut a 2D slice out of a this 3D image and return it
00956                 //               * An alternative to cut_slice
00957                 //               * @param tr orientation of the slice as encapsulated in a Transform object
00958                 //               * @exception ImageDimensionException If this image is not 3D.
00959                 //               * @exception ImageFormatException If this image is complex
00960                 //               * @author David Woolford (adapted from an original version by Steve Ludtke)
00961                 //               * @date Feb 2009
00962                 //               */
00963                 //              EMData* get_cut_slice(const Transform& tr);
00964 
00965         };
00966 
00967 
00968         EMData * operator+(const EMData & em, float n);
00969         EMData * operator-(const EMData & em, float n);
00970         EMData * operator*(const EMData & em, float n);
00971         EMData * operator/(const EMData & em, float n);
00972 
00973         EMData * operator+(float n, const EMData & em);
00974         EMData * operator-(float n, const EMData & em);
00975         EMData * operator*(float n, const EMData & em);
00976         EMData * operator/(float n, const EMData & em);
00977 
00978         EMData * rsub(const EMData & em, float n);
00979         EMData * rdiv(const EMData & em, float n);
00980 
00981         EMData * operator+(const EMData & a, const EMData & b);
00982         EMData * operator-(const EMData & a, const EMData & b);
00983         EMData * operator*(const EMData & a, const EMData & b);
00984         EMData * operator/(const EMData & a, const EMData & b);
00985 
00986 
00987 /*   Next  is Modified by PRB      Transform3D::EMAN,
00988         inline Transform3D EMData::get_transform() const
00989         {
00990                 return Transform3D((float)attr_dict["euler_alt"],
00991                                    (float)attr_dict["euler_az"],
00992                                    (float)attr_dict["euler_phi"]);
00993         }
00994 */
00995 
00996 
00997 }
00998 
00999 
01000 #endif
01001 
01002 /* vim: set ts=4 noet nospell: */

Generated on Tue Jul 12 13:45:47 2011 for EMAN2 by  doxygen 1.4.7