emdata_core.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 emdata__core_h__
#define emdata__core_h__
 
public:
EMData *copy() const;
 
 
EMData *copy_head() const;
 
 
void add(float f,int keepzero=0);
 
 
void add(const EMData & image);
 
void addsquare(const EMData & image);
 
 
void sub(float f);
 
 
void sub(const EMData & image);
 
void subsquare(const EMData & image);
 
 
void mult(int n)
{
        mult((float)n);
}
 
 
void mult(float f);
 
 
void mult(const EMData & image, bool prevent_complex_multiplication=false);
 
void mult_ri(const EMData & image);
 
 
void mult_complex_efficient(const EMData & em, const int radius);
 
void div(float f);
 
 
void div(const EMData & image);
 
void update_min(const EMData & image);
 
 
void to_zero();
 
 
void to_one();
 
void to_value(const float& value);
 
float dot(EMData * with);
 
 
EMData *get_row(int row_index) const;
 
 
void set_row(const EMData * data, int row_index);
 
 
EMData *get_col(int col_index) const;
 
 
void set_col(const EMData * data, int col_index);
 
 
inline float get_value_at(int x, int y, int z) const
{
        return get_data()[(size_t)x + (size_t)y * (size_t)nx + (size_t)z * (size_t)nxy];
}
 
inline float get_value_at_index(size_t i) const
{
        return *(rdata + i);
}
 
inline float get_value_at(int x, int y) const
{
        return get_data()[x + y * nx];
}
 
 
inline float get_value_at(size_t i) const
{
        return get_data()[i];
}
 
std::complex<float> get_complex_at(const int &x,const int &y) const;
 
std::complex<float> get_complex_at(const int &x,const int &y,const int &z) const;
 
size_t get_complex_index(const int &x,const int &y,const int &z) const;
 
size_t get_complex_index(int x,int y,int z,const int &subx0,const int &suby0,const int &subz0,const int &fullnx,const int &fullny,const int &fullnz) const;
 
inline size_t get_complex_index_fast(const int &x,const int &y,const int &z) const {
//      if (abs(x)>=nx/2 || abs(y)>ny/2 || abs(z)>nz/2) return nxyz;
        if (x<0) {
                return (size_t)x*-2+(y<=0?-y:ny-y)*(size_t)nx+(z<=0?-z:nz-z)*(size_t)nxy;
        }
        return x*2+(y<0?ny+y:y)*(size_t)nx+(z<0?nz+z:z)*(size_t)nxy;
}
 
inline void set_complex_at_idx(const int &x,const int &y,const int &z,const std::complex<float> &val) {
        size_t idx=x*2+y*(size_t)nx+z*(size_t)nxy;
        rdata[idx]=(float)val.real();
        rdata[idx+1]=(float)val.imag();
}
 
void set_complex_at(const int &x,const int &y,const std::complex<float> &val);
 
void set_complex_at(const int &x,const int &y,const int &z,const std::complex<float> &val);
 
size_t add_complex_at(const int &x,const int &y,const int &z,const std::complex<float> &val);
 
inline size_t add_complex_at_fast(const int &x,const int &y,const int &z,const std::complex<float> &val) {
//if (x>=nx/2 || y>ny/2 || z>nz/2 || x<=-nx/2 || y<-ny/2 || z<-nz/2) return nxyz;
if (abs(x)>=nx/2 || abs(y)>ny/2 || abs(z)>nz/2) return nxyz;
 
//if (x==0 && abs(y)==16 && abs(z)==1) printf("## %d %d %d\n",x,y,z);
size_t idx;
// for x=0, we need to insert the value in 2 places
if (x==0) {
        if (y==0 && z==0) {
                rdata[0]+=(float)val.real();
                rdata[1]=0;
                return 0;
        }
        // complex conjugate in x=0 plane
        size_t idx=(y<=0?-y:ny-y)*(size_t)nx+(z<=0?-z:nz-z)*(size_t)nxy;
//      if (idx==16*34+1*34*32) printf("a %d %d %d\t%1.5f+%1.5fi\t%1.5f+%1.5fi\n",x,y,z,val.real(),val.imag(),rdata[idx],rdata[idx+1]);
        rdata[idx]+=(float)val.real();
        rdata[idx+1]+=(float)-val.imag();
}
if (abs(x)==nx/2-1) {
        if (y==0 && z==0) {
                rdata[nx-2]+=(float)val.real();
                rdata[nx-1]=0;
                return nx-2;
        }
        // complex conjugate in x=0 plane
        size_t idx=nx-2+(y<=0?-y:ny-y)*(size_t)nx+(z<=0?-z:nz-z)*(size_t)nxy;
//      if (idx==16*34+1*34*32) printf("b %d %d %d\t%1.5f+%1.5fi\t%1.5f+%1.5fi\n",x,y,z,val.real(),val.imag(),rdata[idx],rdata[idx+1]);
        rdata[idx]+=(float)val.real();
        rdata[idx+1]+=(float)-val.imag();
}
if (x<0) {
        idx=-x*2+(y<=0?-y:ny-y)*(size_t)nx+(z<=0?-z:nz-z)*(size_t)nxy;
//      if (idx==16*34+1*34*32) printf("c %d %d %d\t%1.5f+%1.5fi\t%1.5f+%1.5fi\n",x,y,z,val.real(),val.imag(),rdata[idx],rdata[idx+1]);
        rdata[idx]+=(float)val.real();
        rdata[idx+1]+=-(float)val.imag();
        return idx;
}
 
idx=x*2+(y<0?ny+y:y)*(size_t)nx+(z<0?nz+z:z)*(size_t)nxy;
//if (idx==16*34+1*34*32) printf("d %d %d %d\t%1.5f+%1.5fi\t%1.5f+%1.5fi\n",x,y,z,val.real(),val.imag(),rdata[idx],rdata[idx+1]);
rdata[idx]+=(float)val.real();
rdata[idx+1]+=(float)val.imag();
 
return idx;
}
 
size_t add_complex_at(int x,int y,int z,const int &subx0,const int &suby0,const int &subz0,const int &fullnx,const int &fullny,const int &fullnz,const std::complex<float> &val);
 
float get_value_at_wrap(int x, int y, int z) const;
float& get_value_at_wrap(int x, int y, int z);
 
float get_value_at_wrap(int x, int y) const;
float& get_value_at_wrap(int x, int y);
 
float get_value_at_wrap(int x) const;
float& get_value_at_wrap(int x);
 
inline float sget_value_at(Vec3i v) { return sget_value_at(v[0],v[1],v[2]); }
 
float sget_value_at(int x, int y, int z) const;
 
 
float sget_value_at(int x, int y) const;
 
 
float sget_value_at(size_t i) const;
 
 
std::complex<float> get_complex_at_interp(float x, float y) const;
 
std::complex<float> get_complex_at_ginterp(float x, float y) const;
 
std::complex<float> get_complex_at_3ginterp(float x, float y) const;
 
 
float sget_value_at_interp(float x, float y) const;
 
 
float sget_value_at_interp(float x, float y, float z) const;
 
inline void set_value_at(Vec3i loc,float val) { set_value_at(loc[0],loc[1],loc[2],val); }
 
inline void set_value_at(int x, int y, int z, float v)
{
        if( x>=nx || x<0 )
        {
                throw OutofRangeException(0, nx-1, x, "x dimension index");
        }
        else if( y>=ny || y<0 )
        {
                throw OutofRangeException(0, ny-1, y, "y dimension index");
        }
        else if( z>=nz || z<0 )
        {
                throw OutofRangeException(0, nz-1, z, "z dimension index");
        }
        else
        {
                get_data()[(size_t)x + (size_t)y * (size_t)nx + (size_t)z * (size_t)nxy] = v;
                update();
        }
}
 
 
inline void mult_value_at_fast(int x, int y, int z, float v)
{
        get_data()[(size_t)x + (size_t)y * (size_t)nx + (size_t)z * (size_t)nxy] *= v;
        update();
}
 
inline void set_value_at_fast(int x, int y, int z, float v)
{
        get_data()[(size_t)x + (size_t)y * (size_t)nx + (size_t)z * (size_t)nxy] = v;
        update();
}
 
inline void set_value_at_index(size_t i, float v)
{
        *(rdata + i) = v;
}
 
inline void set_value_at(int x, int y, float v)
{
        if( x>=nx || x<0 )
        {
                throw OutofRangeException(0, nx-1, x, "x dimension index");
        }
        else if( y>=ny || y<0 )
        {
                throw OutofRangeException(0, ny-1, y, "y dimension index");
        }
        else
        {
                get_data()[x + y * nx] = v;
                update();
                
        }
}
 
 
inline void set_value_at_fast(int x, int y, float v)
{
        get_data()[x + y * nx] = v;
        update();
}
 
 
inline void set_value_at(int x, float v)
{
        if( x>=nx || x<0 )
        {
                throw OutofRangeException(0, nx-1, x, "x dimension index");
        }
        else
        {
                get_data()[x] = v;
 
                update();
        }
}
 
inline void set_value_at_fast(int x, float v)
{
        get_data()[x] = v;
 
        update();
}
 
 
void free_memory();
 
void free_rdata();
 
EMData & operator+=(float n);
EMData & operator-=(float n);
EMData & operator*=(float n);
EMData & operator/=(float n);
 
EMData & operator+=(const EMData & em);
EMData & operator-=(const EMData & em);
EMData & operator*=(const EMData & em);
EMData & operator/=(const EMData & em);
 
bool operator==(const EMData& that) const;
bool equal(const EMData& that) const;
 
inline float& operator()(const int ix, const int iy, const int iz) const {
        ptrdiff_t pos = (size_t)(ix-xoff) + ((iy-yoff) + (size_t)(iz-zoff)*ny)*nx;
#ifdef BOUNDS_CHECKING
        if (pos < 0 || pos >= (size_t)nx*ny*nz) {
                throw OutofRangeException(0, (size_t)nx*ny*nz-1, pos, "EMData");
        }
#endif // BOUNDS_CHECKING
        return *(get_data() + pos);
}
 
inline float& operator()(const int ix, const int iy) const {
        ptrdiff_t pos = (ix - xoff) + (iy-yoff)*nx;
#ifdef BOUNDS_CHECKING
        if (pos < 0 || pos >= (size_t)nx*ny*nz)
        {
                throw OutofRangeException(0, (size_t)nx*ny*nz-1, pos, "EMData");
        }
#endif // BOUNDS_CHECKING
        return *(get_data() + pos);
}
 
 
inline float& operator()(const size_t ix) const {
        ptrdiff_t pos = ix - xoff;
#ifdef BOUNDS_CHECKING
        if (pos < 0 || pos >= (size_t)nx*ny*nz)
                throw OutofRangeException(0, (size_t)nx*ny*nz-1, pos, "EMData");
#endif // BOUNDS_CHECKING
        return *(get_data() + pos);
}
 
 
void set_array_offsets(const int xoff_=0, const int yoff_=0,
                               const int zoff_=0) {
        xoff=xoff_; yoff=yoff_; zoff=zoff_;
}
 
 
void set_array_offsets(vector<int> offsets) {
        set_array_offsets(offsets[0],offsets[1],offsets[2]);
}
 
 
vector<int> get_array_offsets() {
        vector<int> offsets;
        offsets.push_back(xoff);
        offsets.push_back(yoff);
        offsets.push_back(zoff);
        return offsets;
}
 
 
std::complex<float>& cmplx(const int ix, const int iy, const int iz) {
        ptrdiff_t pos = 2*(ix-xoff)+((iy-yoff)+(iz-zoff)*ny)*(size_t)nx;
#ifdef BOUNDS_CHECKING
        if (pos < 0 || pos >= (size_t)nx*ny*nz)
                throw OutofRangeException(0, (size_t)nx*ny*nz-1, pos, "EMData");
#endif // BOUNDS_CHECKING
        float* begin = get_data() + pos;
        return *(reinterpret_cast<std::complex<float>* >(begin));
}
 
 
std::complex<float>& cmplx(const int ix, const int iy) {
        ptrdiff_t pos = 2*(ix-xoff)+(iy-yoff)*nx;
#ifdef BOUNDS_CHECKING
        if (pos < 0 || pos >= (size_t)nx*ny*nz)
                throw OutofRangeException(0, (size_t)nx*ny*nz-1, pos, "EMData");
#endif // BOUNDS_CHECKING
        float* begin = get_data() + pos;
        return *(reinterpret_cast<std::complex<float>* >(begin));
}
 
 
std::complex<float>& cmplx(const int ix) {
        ptrdiff_t pos = 2*(ix-xoff);
#ifdef BOUNDS_CHECKING
        if (pos < 0 || pos >= (size_t)nx*ny*nz)
                throw OutofRangeException(0, (size_t)nx*ny*nz-1, pos, "EMData");
#endif // BOUNDS_CHECKING
        float* begin = get_data() + pos;
        return *(reinterpret_cast<std::complex<float>* >(begin));
}
 
 
EMData * power(int n) const;
 
EMData * sqrt() const;
 
 
EMData * log() const;
 
 
EMData * log10() const;
 
 
EMData * real() const;
 
 
EMData * imag() const;
 
EMData * absi() const;
 
 
EMData * amplitude() const;
 
 
EMData * phase() const;
 
 
EMData * real2complex(float img = 0.0f) const;
 
#endif  //emdata__core_h__