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 
00456                 EMData *calc_ccfx( EMData * const with, int y0 = 0, int y1 = -1, bool nosum = false, bool flip = false,bool usez=false);
00457 
00458 
00467                 EMData *make_rotational_footprint(bool unwrap = true);
00468                 EMData *make_rotational_footprint_e1(bool unwrap = true);
00469                 EMData *make_rotational_footprint_cmc(bool unwrap = true);
00470 
00487                 EMData *make_footprint(int type=0);
00488 
00489 
00502                 EMData *calc_mutual_correlation(EMData * with, bool tocorner = false, EMData * filter = 0);
00503 
00504 
00520                 EMData *unwrap(int r1 = -1, int r2 = -1, int xs = -1, int dx = 0,
00521                                                            int dy = 0, bool do360 = false, bool weight_radial=true) const;
00522 
00523                 EMData * unwrap_largerR(int r1,int r2,int xs, float rmax_f);
00524 
00525                 EMData *oneDfftPolar(int size, float rmax, float MAXR);
00526 
00527                 
00535                 void apply_radial_func(float x0, float dx, vector < float >array, bool interp = true);
00536 
00537 
00549                 vector < float >calc_radial_dist(int n, float x0, float dx,bool inten);
00550 
00551 
00564                 vector < float >calc_radial_dist(int n, float x0, float dx, int nwedge, float offset,bool inten);
00565 
00566 
00570                 void cconj();
00571 
00572 
00580                 void add_incoherent(EMData * obj);
00581 
00582 
00594                 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);
00595 
00596 
00608                 vector<float> calc_az_dist(int n, float a0, float da, float rmin,
00609                                                                    float rmax);
00610 
00611 #if 0
00612                 void calc_rcf(EMData * with, vector < float >&sum_array);
00613 #endif
00614 
00627                 float calc_dist(EMData * second_img, int y_index = 0) const;
00628 
00645                 EMData *calc_flcf(EMData * with);
00646 
00665                 EMData *calc_fast_sigma_image( EMData* mask);
00666 
00672                 EMData *convolute(EMData * with);
00673 
00674 #if 0
00675                 void create_ctf_map(CtfMapType type, XYData * sf = 0);
00676 #endif
00677 
00678 
00700                 void common_lines(EMData * image1, EMData * image2, int mode = 0,
00701                                                   int steps = 180, bool horizontal = false);
00702 
00713                 void common_lines_real(EMData * image1, EMData * image2,
00714                                                            int steps = 180, bool horizontal = false);
00715 
00729                 void cut_slice(const EMData * const map, const Transform& tr, bool interpolate = true);
00730 
00745                 void uncut_slice(EMData * const map, const Transform& tr) const;
00746 
00755                 EMData *extract_box(const Transform& cs, const Region& r);
00756                 
00760                 int getResolution() const {
00761                         int resolution = 0;
00762                         int num = 1;
00763                         while(num < get_xsize()) {
00764                                 resolution++;
00765                                 num = 1 << resolution;
00766                         }
00767 
00768                         return resolution;
00769                 }
00770                 
00773                 void debug_print_parms()
00774                 {
00775                         std::cout << "Printing EMData params" << std::endl;
00776                         for ( Dict::const_iterator it = attr_dict.begin(); it != attr_dict.end(); ++it )
00777                         {
00778                                 std::cout << (it->first) << " " << (it->second).to_str() << std::endl;
00779                         }
00780                         std::cout << "Done printing EMData params" << std::endl;
00781                 }
00782 
00788                 void set_xyz_origin(float origin_x, float origin_y, float origin_z);
00789                 
00795                 EMData* compute_missingwedge(float wedgeangle, float start = 0.05, float stop = 0.5);
00796                 
00797                 static int totalalloc;
00798         private:
00805                 enum EMDataFlags {
00806 //                      EMDATA_COMPLEX = 1 << 1,
00807 //                      EMDATA_RI = 1 << 2,            // real/imaginary or amp/phase
00808                         EMDATA_BUSY = 1 << 3,      // someone is modifying data, NO LONGER USED
00809                         EMDATA_HASCTFF = 1 << 4,   // has CTF info in the image file
00810                         EMDATA_NEEDUPD = 1 << 5,   // needs a real update
00811 //                      EMDATA_COMPLEXX = 1 << 6,  // 1D fft's in X
00812                         EMDATA_FLIP = 1 << 7,      // is the image flipped
00813                         EMDATA_PAD = 1 << 8,       // is the image fft padded
00814                         EMDATA_FFTODD = 1 << 9,    // is the (real-space) nx odd
00815                         EMDATA_SHUFFLE = 1 << 10,  // fft been shuffled? (so O is centered) PRB
00816                         EMDATA_FH = 1 << 11,        // is the complex image a FH image
00817                         EMDATA_CPU_NEEDS_UPDATE = 1 << 12, // CUDA related: is the CPU version of the image out out data
00818                         EMDATA_GPU_NEEDS_UPDATE = 1 << 13, // CUDA related: is the GPU version of the image out out data
00819                         EMDATA_GPU_RO_NEEDS_UPDATE = 1 << 14 // // CUDA related: is the GPU RO version of the image out out data
00820                 };
00821 
00822                 void update_stat() const;
00823                 void save_byteorder_to_dict(ImageIO * imageio);
00824 
00825         private:
00827                 mutable Dict attr_dict;
00829                 mutable float *rdata;
00831                 float *supp;
00832 
00835                 //Ctf *ctf;
00836 
00838                 mutable int flags;
00839                 // Incremented every time the image changes
00840                 int changecount;
00842                 int nx, ny, nz, nxy;
00843                 size_t nxyz;
00845                 int xoff, yoff, zoff;
00846 
00848                 Vec3f all_translation;
00849 //              Vec3f all_rotation;    /** rotation (az, alt, phi) from the original locaton*/
00850 
00851                 string path;
00852                 int pathnum;
00853 
00855                 mutable EMData* rot_fp;
00856 
00857                 // Clip inplace variables is a local class used from convenience in EMData::clip_inplace
00858                 // Added by d.woolford
00859                 class ClipInplaceVariables
00860                 {
00861                         public:
00862                                 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) :
00863                                         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),
00864                                  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),
00865                                 prv_z_top(0), prv_z_bottom(0), prv_y_back(0), prv_y_front(0), prv_x_left(0), prv_x_right(0)
00866                         {
00867                                 if ( xtrans > 0 ) x_iter -= xtrans;
00868                                 if ( x_iter < 0 ) x_iter = 0;
00869                                 if ( ytrans > 0 ) y_iter -= ytrans;
00870                                 if ( y_iter < 0 ) y_iter = 0;
00871                                 if ( ztrans > 0 ) z_iter -= ztrans;
00872                                 if ( z_iter < 0 ) z_iter = 0;
00873 
00874                                 // Get the depth in the new volume where slices are inserted
00875                                 // if this value is zero it means that the last z-slice in the new
00876                                 // volume contains image data
00877                                 if ( (new_nz + ztrans) > prv_nz ) new_z_top = new_nz + ztrans - prv_nz;
00878                                 if ( (new_ny + ytrans) > prv_ny ) new_y_back = new_ny + ytrans - prv_ny;
00879                                 if ( (new_nx + xtrans) > prv_nx ) new_x_right = new_nx + xtrans - prv_nx;
00880 
00881                                 if ( (new_nz + ztrans) < prv_nz )
00882                                 {
00883                                         prv_z_top = prv_nz - new_nz - ztrans;
00884                                         z_iter -= prv_z_top;
00885                                 }
00886                                 if ( (new_ny + ytrans) < prv_ny )
00887                                 {
00888                                         prv_y_back = prv_ny - new_ny - ytrans;
00889                                         y_iter -= prv_y_back;
00890                                 }
00891                                 if ( (new_nx + xtrans) < prv_nx )
00892                                 {
00893                                         prv_x_right = prv_nx - new_nx - xtrans;
00894                                         x_iter -= prv_x_right;
00895                                 }
00896 
00897                                 if ( xtrans > 0 ) prv_x_left = xtrans;
00898                                 if ( ytrans > 0 ) prv_y_front = ytrans;
00899                                 if ( ztrans > 0 ) prv_z_bottom = ztrans;
00900 
00901                                 if ( xtrans < 0 ) new_x_left = -xtrans;
00902                                 if ( ytrans < 0 ) new_y_front = -ytrans;
00903                                 if ( ztrans < 0 ) new_z_bottom = -ztrans;
00904 
00905                         }
00906                         ~ClipInplaceVariables() {}
00907 
00908                         int prv_nx, prv_ny, prv_nz, new_nx, new_ny, new_nz;
00909                         int xshift, yshift, zshift;
00910                         int x_iter, y_iter, z_iter;
00911                         int new_z_top, new_z_bottom, new_y_back, new_y_front, new_x_left, new_x_right;
00912                         int prv_z_top, prv_z_bottom,  prv_y_back, prv_y_front, prv_x_left, prv_x_right;
00913                 };
00914 
00915 
00916 
00917                 //              /**  Do the Fourier Harmonic Transform  PRB
00918                 //               * Takes a real image, returns the FH
00919                 //               * Sets the EMDATA_FH switch to indicate that it is an FH image
00920                 //               * @exception ImageFormatException If the image is not a square real odd image.
00921                 //               * @return the FH image.
00922                 //               */
00923                 //              EMData* do_FH();
00924 
00925                 //              /**   Do the Inverse Fourier Harmonic Transform   PRB
00926                 //               * Takes an FH image, returns a square  complex image with odd sides
00927                 //               * @exception ImageFormatException If the image is not the FH of something
00928                 //               * @return a square complex image with odd sides
00929                 //               */
00930                 //              EMData* do_FH2F();
00931 
00932 
00933 
00934                 //              /** Caclulates normalization and phase residual for a slice in
00935                 //               * an already existing volume. phase residual is calculated
00936                 //               * using only the inner 1/2 of the fourier sphere. Both the
00937                 //               * slice image and this image must be in complex image format.
00938                 //               *
00939                 //               * @param slice An slice image to be normalized.
00940                 //               * @param orient Orientation of the slice.
00941                 //               * @exception ImageFormatException If the images are not complex.
00942                 //               * @exception ImageDimensionException If the image is 3D.
00943                 //               * @return A float number pair (result, phase-residual).
00944                 //               */
00945                 //              FloatPoint normalize_slice(EMData * slice, const Transform3D & orient);
00946 
00947                 //              /** Caclulates normalization and phase residual for a slice in
00948                 //               * an already existing volume. phase residual is calculated
00949                 //               * using only the inner 1/2 of the fourier sphere. Both the
00950                 //               * slice image and this image must be in complex image format.
00951                 //               *
00952                 //               * @param slice An slice image to be normalized.
00953                 //               * @param alt Orientation euler angle alt (in EMAN convention).
00954                 //               * @param az  Orientation euler angle az  (in EMAN convention).
00955                 //               * @param phi Orientation euler angle phi (in EMAN convention).
00956                 //               * @exception ImageFormatException If the images are not complex.
00957                 //               * @exception ImageDimensionException If the image is 3D.
00958                 //               * @return A float number pair (result, phase-residual).
00959                 //               */
00960                 //              FloatPoint normalize_slice(EMData * slice, float az, float alt, float phi);
00961 
00962                 //              /** Get the normalization and phase residual values
00963                 //               * Used for normalizaton and error measurement when 2D slices are inserted into a 3D volume of Fourier pixels
00964                 //               * Originally added for use by the FourierReconstructor object
00965                 //               * @return the normalization const (pair.first) and the phase residual (pair.second)
00966                 //               * @param slice -the slice to be inserted into the 3D volume
00967                 //               * @param euler - the euler angle orientation of the slice
00968                 //               * @exception ImageDimensionException If this image is not 3D.ImageFormatException
00969                 //               * @exception ImageFormatException If this image is not complex
00970                 //               * @exception ImageFormatException If the slice not complex
00971                 //               */
00972                 //              pair<float, float> get_normalization_and_phaseres( const EMData* const slice, const Transform3D& euler );
00973 
00974                 //              /** cut a 2D slice out of a this 3D image and return it
00975                 //               * An alternative to cut_slice
00976                 //               * @param tr orientation of the slice as encapsulated in a Transform object
00977                 //               * @exception ImageDimensionException If this image is not 3D.
00978                 //               * @exception ImageFormatException If this image is complex
00979                 //               * @author David Woolford (adapted from an original version by Steve Ludtke)
00980                 //               * @date Feb 2009
00981                 //               */
00982                 //              EMData* get_cut_slice(const Transform& tr);
00983 
00984         };
00985 
00986 
00987         EMData * operator+(const EMData & em, float n);
00988         EMData * operator-(const EMData & em, float n);
00989         EMData * operator*(const EMData & em, float n);
00990         EMData * operator/(const EMData & em, float n);
00991 
00992         EMData * operator+(float n, const EMData & em);
00993         EMData * operator-(float n, const EMData & em);
00994         EMData * operator*(float n, const EMData & em);
00995         EMData * operator/(float n, const EMData & em);
00996 
00997         EMData * rsub(const EMData & em, float n);
00998         EMData * rdiv(const EMData & em, float n);
00999 
01000         EMData * operator+(const EMData & a, const EMData & b);
01001         EMData * operator-(const EMData & a, const EMData & b);
01002         EMData * operator*(const EMData & a, const EMData & b);
01003         EMData * operator/(const EMData & a, const EMData & b);
01004 
01005 
01006 /*   Next  is Modified by PRB      Transform3D::EMAN,
01007         inline Transform3D EMData::get_transform() const
01008         {
01009                 return Transform3D((float)attr_dict["euler_alt"],
01010                                    (float)attr_dict["euler_az"],
01011                                    (float)attr_dict["euler_phi"]);
01012         }
01013 */
01014 
01015 
01016 }
01017 
01018 
01019 #endif
01020 
01021 /* vim: set ts=4 noet nospell: */

Generated on Tue Jun 11 12:40:21 2013 for EMAN2 by  doxygen 1.4.7