Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

emdata_metadata.cpp

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 #include "emdata.h"
00037 #include "ctf.h"
00038 #include "portable_fileio.h"
00039 #include "imageio.h"
00040 
00041 #include <cstring>
00042 #include <sstream>
00043 using std::stringstream;
00044 
00045 #include <iomanip>
00046 using std::setprecision;
00047 using namespace EMAN;
00048 
00049 #ifdef EMAN2_USING_CUDA
00050 #include "cuda/cuda_processor.h"
00051 #endif // EMAN2_USING_CUDA
00052 
00053 EMData* EMData::get_fft_amplitude2D()
00054 {
00055         ENTERFUNC;
00056 
00057 //      int ndim = get_ndim();
00058         if (!is_complex()) {
00059                 LOGERR("complex image expected. Input image is real image.");
00060                 throw ImageFormatException("complex image expected. Input image is a real image.");
00061         }
00062         if (nz>1) {
00063                 LOGERR("2D image expected. Input image is 3D");
00064                 throw ImageFormatException("2D odd square complex image"
00065                         " expected Input image is 3D.");
00066         }
00067 
00068         int nx2 = nx/2;
00069 
00070         EMData *dat = copy_head();
00071 
00072         dat->set_size(nx2, ny, nz);
00073         dat->to_zero();
00074 
00075         float temp=0;
00076 
00077         for (int j = 0; j < ny; j++) {
00078                 for (int i = 0; i < nx2; i++) {
00079                         temp = (*this)(2*i,j)*(*this)(2*i,j);
00080                         temp += (*this)(2*i+1,j)*(*this)(2*i+1,j);
00081                         (*dat)(i,j) = std::sqrt(temp);
00082                 }
00083         }
00084 
00085         dat->update();
00086         dat->set_complex(false);
00087         dat->set_ri(false);
00088 
00089         EXITFUNC;
00090         return dat;
00091 }
00092 
00093 
00094 EMData* EMData::get_fft_amplitude()
00095 {
00096         ENTERFUNC;
00097 
00098         if (!is_complex()) {
00099                 LOGERR("complex image expected. Input image is real image.");
00100                 throw ImageFormatException("complex image expected. Input image is a real image.");
00101         }
00102 
00103         ri2ap();
00104 
00105         int nx2 = nx - 2;
00106         EMData *dat = copy_head();
00107         dat->set_size(nx2, ny, nz);
00108         dat->to_zero();
00109 
00110         float *d = dat->get_data();
00111         float *data = get_data();
00112         int ndim = get_ndim();
00113 
00114         size_t idx1, idx2, idx3;
00115         if (ndim == 3) {
00116                 for (int k = 1; k < nz; ++k) {
00117                         for (int j = 1; j < ny; ++j) {
00118                                 for (int i = 0; i < nx2/2; ++i) {
00119                                         idx1 = (size_t)k*nx2*ny+j*nx2+nx2/2+i;
00120                                         idx2 = (size_t)k*nx*ny+j*nx+2*i;
00121                                         idx3 = (size_t)(nz-k)*nx2*ny+(ny-j)*nx2+nx2/2-i;
00122                                         d[idx1] = data[idx2];
00123                                         d[idx3] = data[idx2];
00124                                 }
00125                         }
00126                 }
00127         }
00128         else if (ndim == 2) {
00129                 for (int j = 1; j < ny; ++j) {
00130                         for (int i = 0; i < nx2/2; ++i) {
00131                                 d[j*nx2+nx2/2+i] = data[j*nx+2*i];
00132                                 d[(ny-j)*nx2+nx2/2-i] = data[j*nx+2*i];
00133                         }
00134                 }
00135         }
00136 
00137         dat->update();
00138         dat->set_complex(false);
00139         if(dat->get_ysize()==1 && dat->get_zsize()==1) {
00140                 dat->set_complex_x(false);
00141         }
00142         dat->set_ri(false);
00143 
00144         EXITFUNC;
00145         return dat;
00146 }
00147 
00148 
00149 EMData* EMData::get_fft_phase()
00150 {
00151         ENTERFUNC;
00152 
00153         if (!is_complex()) {
00154                 LOGERR("complex image expected. Input image is real image.");
00155                 throw ImageFormatException("complex image expected. Input image is a real image.");
00156         }
00157 
00158         ri2ap();
00159 
00160         int nx2 = nx - 2;
00161         EMData *dat = copy_head();
00162         dat->set_size(nx2, ny, nz);
00163         dat->to_zero();
00164 
00165         float *d = dat->get_data();
00166         float * data = get_data();
00167 
00168         int ndim = get_ndim();
00169         size_t idx1, idx2, idx3;
00170         if (ndim == 3) {
00171                 for (int k = 1; k < nz; ++k) {
00172                         for (int j = 1; j < ny; ++j) {
00173                                 for (int i = 0; i < nx2/2; ++i) {
00174                                         idx1 = (size_t)k*nx2*ny+j*nx2+nx2/2+i;
00175                                         idx2 = (size_t)k*nx*ny+j*nx+2*i+1;
00176                                         idx3 = (size_t)(nz-k)*nx2*ny+(ny-j)*nx2+nx2/2-i;
00177                                         d[idx1] = data[idx2];
00178                                         d[idx3] = -data[idx2];
00179                                 }
00180                         }
00181                 }
00182         }
00183         else if (ndim == 2) {
00184                 for (int j = 1; j < ny; ++j) {
00185                         for (int i = 0; i < nx2/2; ++i) {
00186                                 d[j*nx2+nx2/2+i] = data[j*nx+2*i+1];
00187                                 d[(ny-j)*nx2+nx2/2-i] = -data[j*nx+2*i+1];
00188                         }
00189                 }
00190         }
00191 
00192         dat->update();
00193         dat->set_complex(false);
00194         if(dat->get_ysize()==1 && dat->get_zsize()==1) {
00195                 dat->set_complex_x(false);
00196         }
00197         dat->set_ri(false);
00198 
00199         EXITFUNC;
00200         return dat;
00201 }
00202 
00203 #include <sys/stat.h>
00204 
00205 void EMData::write_data(string fsp,size_t loc,const Region* area,const int file_nx, const int file_ny, const int file_nz) {
00206 
00207         if (area) {
00208                 struct stat fileinfo;
00209                 if ( stat(fsp.c_str(),&fileinfo) != 0 ) throw UnexpectedBehaviorException("To write an image using a region the file must already exist and be the correct dimensions");
00210         }
00211 
00212 
00213         FILE *f = 0;
00214         f=fopen(fsp.c_str(), "rb+");
00215         if (!f) f=fopen(fsp.c_str(), "wb");
00216         if (!f) throw FileAccessException(fsp);
00217         portable_fseek(f,loc,SEEK_SET);
00218         if (!area) {
00219                 if (fwrite(get_data(),nx*ny,nz*4,f)!=(size_t)(nz*4)) throw FileAccessException(fsp);
00220         } else {
00221                 int fnx = nx;
00222                 if (file_nx != 0) fnx = file_nx;
00223                 int fny = ny;
00224                 if (file_ny != 0) fny = file_ny;
00225                 int fnz = nz;
00226                 if (file_nz != 0) fnz = file_nz;
00227 
00228                 EMUtil::process_region_io(get_data(), f, ImageIO::READ_WRITE,
00229                                                                   0, 4,fnx,fny,fnz,area);
00230         }
00231         fclose(f);
00232 //      printf("WROTE %s %ld %ld\n",fsp.c_str(),loc,nx*ny*nz);
00233 }
00234 
00235 void EMData::read_data(string fsp,size_t loc,const Region* area, const int file_nx, const int file_ny, const int file_nz) {
00236         FILE *f = 0;
00237         f=fopen(fsp.c_str(), "rb");
00238         if (!f) throw FileAccessException(fsp);
00239         int fnx = nx;
00240         if (file_nx != 0) fnx = file_nx;
00241         int fny = ny;
00242         if (file_ny != 0) fny = file_ny;
00243         int fnz = nz;
00244         if (file_nz != 0) fnz = file_nz;
00245 
00246         portable_fseek(f,loc,SEEK_SET);
00247         EMUtil::process_region_io(get_data(), f, ImageIO::READ_ONLY,
00248                                                           0, 4,fnx,fny,fnz,area);
00249 //      portable_fseek(f,loc,SEEK_SET);
00250 //      if (fread(get_data(),nx*ny,nz*4,f)!=(size_t)(nz*4)) throw FileAccessException(fsp);
00251         fclose(f);
00252 }
00253 
00254 float EMData::calc_center_density()
00255 {
00256         ENTERFUNC;
00257 
00258         float center = get_attr("mean");
00259         float sigma = get_attr("sigma");
00260         float ds = sigma / 2;
00261         size_t size = (size_t)nx * ny * nz;
00262         float *d = get_data();
00263         float sigma1 = sigma / 20;
00264         float sigma2 = sigma / 1000;
00265 
00266         while (ds > sigma1) {
00267                 double sum = 0;
00268                 int norm = 0;
00269 
00270                 for (size_t i = 0; i < size; ++i) {
00271                         if (fabs(d[i] - center) < ds) {
00272                                 sum += d[i];
00273                                 norm++;
00274                         }
00275                 }
00276                 if (!norm) {
00277                         break;
00278                 }
00279                 float mean = (float)(sum / norm);
00280                 if (fabs(mean - center) < sigma2) {
00281                         ds *= 0.5f;
00282                 }
00283                 center = mean;
00284         }
00285         EXITFUNC;
00286 
00287         return center;
00288 }
00289 
00290 
00291 float EMData::calc_sigma_diff()
00292 {
00293         ENTERFUNC;
00294 
00295         float *d = get_data();
00296         float mean = get_attr("mean");
00297         float sigma = get_attr("sigma");
00298 
00299         double sum_up = 0;
00300         double sum_down = 0;
00301         int nup = 0;
00302         int ndown = 0;
00303 
00304         size_t size = (size_t)nx * ny * nz;
00305 
00306         for (size_t i = 0; i < size; ++i) {
00307                 if (d[i] > mean) {
00308                         sum_up += Util::square(d[i] - mean);
00309                         nup++;
00310                 }
00311                 else {
00312                         sum_down += Util::square(mean - d[i]);
00313                         ndown++;
00314                 }
00315         }
00316 
00317         float sigup = std::sqrt((float)sum_up / nup);
00318         float sigdown = std::sqrt((float)sum_down / ndown);
00319         float sig_diff = fabs(sigup - sigdown) / sigma;
00320 
00321 
00322         EXITFUNC;
00323         return sig_diff;
00324 
00325 }
00326 
00327 
00328 IntPoint EMData::calc_min_location() const
00329 {
00330         ENTERFUNC;
00331 
00332         int di = 1;
00333         if (is_complex() && !is_ri()) {
00334                 di = 2;
00335         }
00336 
00337         float min = FLT_MAX;
00338         int min_x = 0;
00339         int min_y = 0;
00340         int min_z = 0;
00341         int nxy = nx * ny;
00342         float * data = get_data();
00343 
00344         for (int j = 0; j < nz; ++j) {
00345                 size_t cur_z = (size_t)j * nxy;
00346 
00347                 for (int k = 0; k < ny; ++k) {
00348                         size_t cur_y = k * nx + cur_z;
00349 
00350                         for (int l = 0; l < nx; l += di) {
00351                                 float t = data[l + cur_y];
00352                                 if (t < min) {
00353                                         min_x = l;
00354                                         min_y = k;
00355                                         min_z = j;
00356                                         min = t;
00357                                 }
00358                         }
00359                 }
00360         }
00361 
00362         return IntPoint(min_x, min_y, min_z);
00363 }
00364 
00365 
00366 IntPoint EMData::calc_max_location() const
00367 {
00368         ENTERFUNC;
00369 
00370         int di = 1;
00371         if (is_complex() && !is_ri()) {
00372                 di = 2;
00373         }
00374 
00375         float max = -FLT_MAX;
00376         int max_x = 0;
00377         int max_y = 0;
00378         int max_z = 0;
00379         int nxy = nx * ny;
00380         float * data = get_data();
00381 
00382         for (int j = 0; j < nz; ++j) {
00383                 size_t cur_z = (size_t)j * nxy;
00384 
00385                 for (int k = 0; k < ny; ++k) {
00386                         size_t cur_y = k * nx + cur_z;
00387 
00388                         for (int l = 0; l < nx; l += di) {
00389                                 float t = data[l + cur_y];
00390                                 if (t > max) {
00391                                         max_x = l;
00392                                         max_y = k;
00393                                         max_z = j;
00394                                         max = t;
00395                                 }
00396                         }
00397                 }
00398         }
00399 
00400         EXITFUNC;
00401         return IntPoint(max_x, max_y, max_z);
00402 }
00403 
00404 
00405 IntPoint EMData::calc_max_location_wrap(const int maxdx, const int maxdy, const int maxdz)
00406 {
00407         int maxshiftx = maxdx, maxshifty = maxdy, maxshiftz = maxdz;
00408         if (maxdx == -1) maxshiftx = get_xsize()/4;
00409         if (maxdy == -1) maxshifty = get_ysize()/4;
00410         if (maxdz == -1) maxshiftz = get_zsize()/4;
00411 
00412         float max_value = -FLT_MAX;
00413 
00414         IntPoint peak(0,0,0);
00415         
00416 #ifdef EMAN2_USING_CUDA //CUDA
00417         if(EMData::usecuda == 1 && cudarwdata){
00418         
00419                 CudaPeakInfo* soln = calc_max_location_wrap_cuda(cudarwdata, nx, ny, nz, maxdx, maxdy, maxdz);
00420                 
00421                 peak[0] = soln->px;
00422                 peak[1] = soln->py;
00423                 peak[2] = soln->pz;
00424                 free(soln);
00425                 
00426 //              cout << "x " << peak[0] << " y " << peak[1] << " z " << peak[2] << endl;
00427                 return peak;
00428         }
00429 #endif
00430         for (int k = -maxshiftz; k <= maxshiftz; k++) {
00431                 for (int j = -maxshifty; j <= maxshifty; j++) {
00432                         for (int i = -maxshiftx; i <= maxshiftx; i++) {
00433 
00434                                 float value = get_value_at_wrap(i,j,k);
00435 
00436                                 if (value > max_value) {
00437                                         max_value = value;
00438                                         peak[0] = i;
00439                                         peak[1] = j;
00440                                         peak[2] = k;
00441                                 }
00442                         }
00443                 }
00444         }
00445 
00446         return peak;
00447 }
00448 
00449 vector<float> EMData::calc_max_location_wrap_intp(const int maxdx, const int maxdy, const int maxdz)
00450 {
00451         int maxshiftx = maxdx, maxshifty = maxdy, maxshiftz = maxdz;
00452         if (maxdx == -1) maxshiftx = get_xsize()/4;
00453         if (maxdy == -1) maxshifty = get_ysize()/4;
00454         if (maxdz == -1) maxshiftz = get_zsize()/4;
00455 
00456         float max_value = -FLT_MAX;
00457 
00458         IntPoint peak(0,0,0);
00459 
00460 // NOT yet working......
00477         for (int k = -maxshiftz; k <= maxshiftz; k++) {
00478                 for (int j = -maxshifty; j <= maxshifty; j++) {
00479                         for (int i = -maxshiftx; i <= maxshiftx; i++) {
00480 
00481                                 float value = get_value_at_wrap(i,j,k);
00482 
00483                                 if (value > max_value) {
00484                                         max_value = value;
00485                                         peak[0] = i;
00486                                         peak[1] = j;
00487                                         peak[2] = k;
00488                                 }
00489                         }
00490                 }
00491         }
00492         
00493         // compute the center of mass
00494         float cmx = 0.0; float cmy = 0.0f; float cmz = 0.0f;
00495         float sval = 0.0f;
00496         for (float x = float(peak[0])-2.0f; x <= float(peak[0])+2.0f; x++) {
00497                 for (float y = float(peak[1])-2.0f; y <= float(peak[1])+2.0f; y++) {
00498                         for (float z = float(peak[2])-2.0f; z <= float(peak[2])+2.0f; z++) {
00499                                 //Compute center of mass
00500                                 float val = get_value_at_wrap(x,y,z);
00501                                 cmx += x*val;
00502                                 cmy += y*val;
00503                                 cmz += z*val;
00504                                 sval += val;
00505                         }
00506                 }
00507         }
00508         cmx /= sval;
00509         cmy /= sval;
00510         cmz /= sval;
00511 
00512         vector<float> mydata;
00513         mydata.push_back(cmx);
00514         mydata.push_back(cmy);
00515         mydata.push_back(cmz);
00516         mydata.push_back(max_value);
00517 
00564         return mydata;
00565 }
00566 
00567 FloatPoint EMData::calc_center_of_mass(float threshold)
00568 {
00569         float *data = get_data();
00570 
00571         //float sigma = get_attr("sigma");
00572         //float mean = get_attr("mean");
00573         float m = 0.0;
00574 
00575         FloatPoint com(0,0,0);
00576         for (int i = 0; i < nx; ++i) {
00577                 for (int j = 0; j < ny; ++j) {
00578                         int j2 = nx * j;
00579                         for (int k = 0; k < nz; ++k) {
00580                                 size_t l = i + j2 + (size_t)k * nxy;
00581                                 if (data[l] >= threshold) {             // threshold out noise (and negative density)
00582                                         com[0] += i * data[l];
00583                                         com[1] += j * data[l];
00584                                         com[2] += k * data[l];
00585                                         m += data[l];
00586                                 }
00587                         }
00588                 }
00589         }
00590 
00591         com[0] /= m;
00592         com[1] /= m;
00593         com[2] /= m;
00594 
00595         return com;
00596 }
00597 
00598 
00599 size_t EMData::calc_min_index() const
00600 {
00601         IntPoint min_location = calc_min_location();
00602         size_t i = min_location[0] + min_location[1] * nx + (size_t)min_location[2] * nx * ny;
00603         return i;
00604 }
00605 
00606 
00607 size_t EMData::calc_max_index() const
00608 {
00609         IntPoint max_location = calc_max_location();
00610         size_t i = max_location[0] + max_location[1] * nx + (size_t)max_location[2] * nx * ny;
00611         return i;
00612 }
00613 
00614 
00615 vector<Pixel> EMData::calc_highest_locations(float threshold) const 
00616 {
00617         ENTERFUNC;
00618 
00619         vector<Pixel> result;
00620 
00621         int di = 1;
00622         if (is_complex() && !is_ri()) {
00623                 di = 2;
00624         }
00625 
00626         int nxy = nx * ny;
00627         float * data = get_data();
00628 
00629         for (int j = 0; j < nz; ++j) {
00630                 size_t cur_z = (size_t)j * nxy;
00631 
00632                 for (int k = 0; k < ny; ++k) {
00633                         size_t cur_y = k * nx + cur_z;
00634 
00635                         for (int l = 0; l < nx; l += di) {
00636                                 float v =data[l + cur_y];
00637                                 if (v > threshold) {
00638                                         result.push_back(Pixel(l, k, j, v));
00639                                 }
00640                         }
00641                 }
00642         }
00643 
00644         std::sort(result.begin(), result.end());
00645 
00646         EXITFUNC;
00647         return result;
00648 }
00649 
00650 vector<Pixel> EMData::calc_n_highest_locations(int n)
00651 {
00652         ENTERFUNC;
00653 
00654         vector<Pixel> result;
00655 
00656         int di = 1;
00657         if (is_complex() && !is_ri()) {
00658                 di = 2;
00659         }
00660 
00661         // initialize with n elements
00662         float * data = get_data();
00663         for ( int i=0; i<n; i++) result.push_back(Pixel(0,0,0,data[0]));
00664 
00665         int nxy = nx * ny;
00666 
00667         for (int j = 0; j < nz; ++j) {
00668                 size_t cur_z = (size_t)j * nxy;
00669 
00670                 for (int k = 0; k < ny; ++k) {
00671                         size_t cur_y = k * nx + cur_z;
00672 
00673                         for (int l = 0; l < nx; l += di) {
00674                                 float v =data[l + cur_y];
00675                                 if (v<result[n-1].value) continue;
00676                                 for (vector<Pixel>::iterator i=result.begin(); i<result.end(); i++) {
00677                                         if (v>(*i).value) { result.insert(i,Pixel(l, k, j, v)); result.pop_back(); break; }
00678                                 }
00679                         }
00680                 }
00681         }
00682 
00683         EXITFUNC;
00684         return result;
00685 }
00686 
00687 vector<Pixel> EMData::find_pixels_with_value(float val) 
00688 {
00689         ENTERFUNC;
00690         
00691         if ( is_complex() ) throw ImageFormatException("Error - find_pixels_with_value real only");
00692 
00693         vector<Pixel> result;
00694 
00695         for (int k = 0; k < nz; k++) {
00696                 for (int j = 0; j < ny; j++) {
00697                         for (int i = 0; i < nx; i++) {
00698                                 if (get_value_at(i,j,k)==val) result.push_back(Pixel(i,j,k,val));
00699                         }
00700                 }
00701         }
00702 
00703         EXITFUNC;
00704         return result;
00705 }
00706 
00707 float EMData::get_edge_mean() const
00708 {
00709         ENTERFUNC;
00710 #ifdef EMAN2_USING_CUDA
00711         if(EMData::usecuda == 1 && cudarwdata){
00712                 
00713                         return get_edgemean_cuda(cudarwdata, nx, ny, nz);
00714                         
00715         }
00716 #endif
00717         int di = 0;
00718         double edge_sum = 0;
00719         float edge_mean = 0;
00720         size_t nxy = nx * ny;
00721         float * data = get_data();
00722         if (nz == 1) {
00723                 for (int i = 0, j = (ny - 1) * nx; i < nx; ++i, ++j) {
00724                         edge_sum += data[i] + data[j];
00725                 }
00726                 for (size_t i = 0, j = nx - 1; i < nxy; i += nx, j += nx) {
00727                         edge_sum += data[i] + data[j];
00728                 }
00729                 edge_mean = (float)edge_sum / (nx * 2 + ny * 2);
00730         }
00731         else {
00732                 if (nx == ny && nx == nz * 2 - 1) {
00733                         for (size_t j = (nxy * (nz - 1)); j < nxy * nz; ++j, ++di) {
00734                                 edge_sum += data[j];
00735                         }
00736                 }
00737                 else {
00738                         for (size_t i = 0, j = (nxy * (nz - 1)); i < nxy; ++i, ++j, ++di) {
00739                                 edge_sum += data[i] + data[j];
00740                         }
00741                 }
00742 
00743                 int nxy2 = nx * (ny - 1);
00744                 for (int k = 1; k < nz - 1; ++k) {
00745                         size_t k2 = k * nxy;
00746                         size_t k3 = k2 + nxy2;
00747                         for (int i = 0; i < nx; ++i, ++di) {
00748                                 edge_sum += data[i + k2] + data[i + k3];
00749                         }
00750                 }
00751                 for (int k = 1; k < nz - 1; ++k) {
00752                         size_t k2 = k * nxy;
00753                         size_t k3 = nx - 1 + k2;
00754                         for (int i = 1; i < ny - 1; ++i, ++di) {
00755                                 edge_sum += data[i * nx + k2] + data[i * nx + k3];
00756                         }
00757                 }
00758 
00759                 edge_mean = (float)edge_sum / (di * 2);
00760         }
00761         EXITFUNC;
00762 
00763         return edge_mean;
00764 }
00765 
00766 
00767 float EMData::get_circle_mean()
00768 {
00769         ENTERFUNC;
00770 
00771         static bool busy = false;
00772         static EMData *mask = 0;
00773 
00774         while (busy);
00775         busy = true;
00776 
00777         if (!mask || !EMUtil::is_same_size(this, mask)) {
00778                 if (!mask) {
00779                         mask = new EMData();
00780                 }
00781                 mask->set_size(nx, ny, nz);
00782                 mask->to_one();
00783 
00784                 float radius = (float)(ny / 2 - 2);
00785                 mask->process_inplace("mask.sharp", Dict("inner_radius", radius - 1,
00786                                                                            "outer_radius", radius + 1));
00787 
00788         }
00789         double n = 0,s=0;
00790         float *d = mask->get_data();
00791         float * data = get_data();
00792         size_t size = (size_t)nx*ny*nz;
00793         for (size_t i = 0; i < size; ++i) {
00794                 if (d[i]) { n+=1.0; s+=data[i]; }
00795         }
00796 
00797 
00798         float result = (float)(s/n);
00799         busy = false;
00800 
00801         EXITFUNC;
00802         return result;
00803 }
00804 
00805 
00806 void EMData::set_ctf(Ctf * new_ctf)
00807 {
00808         ENTERFUNC;
00809 
00810         vector<float> vctf = new_ctf->to_vector();
00811         attr_dict["ctf"] = vctf;
00812 
00813         EXITFUNC;
00814 }
00815 
00816 Ctf * EMData::get_ctf() const
00817 {
00818         if(attr_dict.has_key("ctf")) {
00819                 EMAN1Ctf * ctf = new EMAN1Ctf();
00820                 ctf->from_vector(attr_dict["ctf"]);
00821 
00822                 return dynamic_cast<Ctf *>(ctf);
00823         }
00824         else {
00825                 return 0;
00826         }
00827 }
00828 
00829 #include <iostream>
00830 using std::cout;
00831 using std::endl;
00832 
00833 void EMData::set_size(int x, int y, int z, bool noalloc)
00834 {
00835         ENTERFUNC;
00836 
00837         if (x <= 0) {
00838                 throw InvalidValueException(x, "x size <= 0");
00839         }
00840         else if (y <= 0) {
00841                 throw InvalidValueException(y, "y size <= 0");
00842         }
00843         else if (z <= 0) {
00844                 throw InvalidValueException(z, "z size <= 0");
00845         }
00846 
00847 #ifdef MEMDEBUG2
00848         printf("EMDATA sz  %4d    %p (%d,%d,%d)\n",EMData::totalalloc,this,x,y,z);
00849 #endif
00850 
00851 
00852         int old_nx = nx;
00853 
00854         size_t size = (size_t)x*y*z*sizeof(float);
00855         
00856         if (noalloc) {
00857                 nx = x;
00858                 ny = y;
00859                 nz = z;
00860                 nxy = nx*ny;
00861                 nxyz = (size_t)nx*ny*nz;
00862                 return;
00863         }
00864         
00865         if (rdata != 0) {
00866                 rdata = (float*)EMUtil::em_realloc(rdata,size);
00867         } else {
00868                 // Just pass on this for a while....see what happens
00869                 rdata = (float*)EMUtil::em_malloc(size);
00870         }
00871 //      rdata = static_cast < float *>(realloc(rdata, size));
00872         if ( rdata == 0 )
00873         {
00874                 stringstream ss;
00875                 string gigs;
00876                 ss << (float) size/1000000000.0;
00877                 ss >> gigs;
00878                 string message = "Cannot allocate " + gigs + " GB - not enough memory.";
00879                 throw BadAllocException(message);
00880         }
00881 
00882         nx = x;
00883         ny = y;
00884         nz = z;
00885         nxy = nx*ny;
00886         nxyz = (size_t)nx*ny*nz;
00887 
00888 // once the object is resized, the CUDA need to be updated
00889 #ifdef EMAN2_USING_CUDA
00890         if(cudarwdata) {
00891                 //cout << "rw free on set size" << endl;
00892                 rw_free();
00893                 rw_alloc();
00894         }
00895         if(cudarodata) {
00896                 ro_free();
00897                 ro_alloc();
00898         }
00899 #endif // EMAN2_USING_CUDA
00900 
00901         if (old_nx == 0) {
00902                 EMUtil::em_memset(get_data(),0,size);
00903         }
00904 
00905         if (supp) {
00906                 EMUtil::em_free(supp);
00907                 supp = 0;
00908         }
00909 
00910         update();
00911         EXITFUNC;
00912 }
00913 
00914 #ifdef EMAN2_USING_CUDA
00915 
00916 void EMData::set_size_cuda(int x, int y, int z)
00917 {
00918         ENTERFUNC;
00919 
00920         if (x <= 0) {
00921                 throw InvalidValueException(x, "x size <= 0");
00922         }
00923         else if (y <= 0) {
00924                 throw InvalidValueException(y, "y size <= 0");
00925         }
00926         else if (z <= 0) {
00927                 throw InvalidValueException(z, "z size <= 0");
00928         }
00929 
00930         nx = x;
00931         ny = y;
00932         nz = z;
00933 
00934         nxy = nx*ny;
00935         nxyz = (size_t)nx*ny*nz;
00936 
00937 //      gpu_update();
00938 
00939         EXITFUNC;
00940 }
00941 
00942 #endif // EMAN2_USING_CUDA
00943 
00944 MArray2D EMData::get_2dview() const
00945 {
00946         const int ndims = 2;
00947         if (get_ndim() != ndims) {
00948                 throw ImageDimensionException("2D only");
00949         }
00950         boost::array<std::size_t,ndims> dims = {{nx, ny}};
00951         MArray2D marray(get_data(), dims, boost::fortran_storage_order());
00952         return marray;
00953 }
00954 
00955 
00956 MArray3D EMData::get_3dview() const
00957 {
00958         const int ndims = 3;
00959         boost::array<std::size_t,ndims> dims = {{nx, ny, nz}};
00960         MArray3D marray(get_data(), dims, boost::fortran_storage_order());
00961         return marray;
00962 }
00963 
00964 
00965 MCArray2D EMData::get_2dcview() const
00966 {
00967         const int ndims = 2;
00968         if (get_ndim() != ndims) {
00969                 throw ImageDimensionException("2D only");
00970         }
00971         boost::array<std::size_t,ndims> dims = {{nx/2, ny}};
00972         std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
00973         MCArray2D marray(cdata, dims, boost::fortran_storage_order());
00974         return marray;
00975 }
00976 
00977 
00978 MCArray3D EMData::get_3dcview() const
00979 {
00980         const int ndims = 3;
00981         boost::array<std::size_t,ndims> dims = {{nx/2, ny, nz}};
00982         std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
00983         MCArray3D marray(cdata, dims, boost::fortran_storage_order());
00984         return marray;
00985 }
00986 
00987 
00988 MCArray3D* EMData::get_3dcviewptr() const
00989 {
00990         const int ndims = 3;
00991         boost::array<std::size_t,ndims> dims = {{nx/2, ny, nz}};
00992         std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
00993         MCArray3D* marray = new MCArray3D(cdata, dims,
00994                                                                           boost::fortran_storage_order());
00995         return marray;
00996 }
00997 
00998 
00999 MArray2D EMData::get_2dview(int x0, int y0) const
01000 {
01001         const int ndims = 2;
01002         if (get_ndim() != ndims) {
01003                 throw ImageDimensionException("2D only");
01004         }
01005         boost::array<std::size_t,ndims> dims = {{nx, ny}};
01006         MArray2D marray(get_data(), dims, boost::fortran_storage_order());
01007         boost::array<std::size_t,ndims> bases={{x0, y0}};
01008         marray.reindex(bases);
01009         return marray;
01010 }
01011 
01012 
01013 MArray3D EMData::get_3dview(int x0, int y0, int z0) const
01014 {
01015         const int ndims = 3;
01016         boost::array<std::size_t,ndims> dims = {{nx, ny, nz}};
01017         MArray3D marray(get_data(), dims, boost::fortran_storage_order());
01018         boost::array<std::size_t,ndims> bases={{x0, y0, z0}};
01019         marray.reindex(bases);
01020         return marray;
01021 }
01022 
01023 
01024 MCArray2D EMData::get_2dcview(int x0, int y0) const
01025 {
01026         const int ndims = 2;
01027         if (get_ndim() != ndims) {
01028                 throw ImageDimensionException("2D only");
01029         }
01030         boost::array<std::size_t,ndims> dims = {{nx/2, ny}};
01031         std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
01032         MCArray2D marray(cdata, dims, boost::fortran_storage_order());
01033         boost::array<std::size_t,ndims> bases={{x0, y0}};
01034         marray.reindex(bases);
01035         return marray;
01036 }
01037 
01038 
01039 MCArray3D EMData::get_3dcview(int x0, int y0, int z0) const
01040 {
01041         const int ndims = 3;
01042         boost::array<std::size_t,ndims> dims = {{nx/2, ny, nz}};
01043         std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
01044         MCArray3D marray(cdata, dims, boost::fortran_storage_order());
01045         boost::array<std::size_t,ndims> bases={{x0, y0, z0}};
01046         marray.reindex(bases);
01047         return marray;
01048 }
01049 
01050 int greaterthan( const void* p1, const void* p2 )
01051 {
01052         float*  v1 = (float*) p1;
01053         float*  v2 = (float*) p2;
01054 
01055         if ( *v1 < *v2 )  return 0;
01056         else return 1;
01057 }
01058 
01059 
01060 EMObject EMData::get_attr(const string & key) const
01061 {
01062         ENTERFUNC;
01063         
01064         if ((flags & EMDATA_NEEDUPD) && (key != "is_fftpad") && (key != "xform.align2d")){update_stat();} //this gives a spped up of 7.3% according to e2speedtest
01065                 
01066         size_t size = (size_t)nx * ny * nz;
01067         if (key == "kurtosis") {
01068                 float mean = attr_dict["mean"];
01069                 float sigma = attr_dict["sigma"];
01070 
01071                 float *data = get_data();
01072                 double kurtosis_sum = 0;
01073 
01074                 for (size_t k = 0; k < size; ++k) {
01075                         float t = (data[k] - mean) / sigma;
01076                         float tt = t * t;
01077                         kurtosis_sum += tt * tt;
01078                 }
01079 
01080                 float kurtosis = (float)(kurtosis_sum / size - 3.0);
01081                 return kurtosis;
01082         }
01083         else if (key == "skewness") {
01084                 float mean = attr_dict["mean"];
01085                 float sigma = attr_dict["sigma"];
01086 
01087                 float *data = get_data();
01088                 double skewness_sum = 0;
01089                 for (size_t k = 0; k < size; ++k) {
01090                         float t = (data[k] - mean) / sigma;
01091                         skewness_sum +=  t * t * t;
01092                 }
01093                 float skewness = (float)(skewness_sum / size);
01094                 return skewness;
01095         }
01096         else if (key == "median")
01097         {
01098                 if ( is_complex() ) throw ImageFormatException("Error - can not calculate the median of a complex image");
01099                 size_t n = size;
01100                 float* tmp = new float[n];
01101                 float* d = get_data();
01102                 if (tmp == 0 ) throw BadAllocException("Error - could not create deep copy of image data");
01103 //              for(size_t i=0; i < n; ++i) tmp[i] = d[i]; // should just be a memcpy
01104                 std::copy(d, d+n, tmp);
01105                 qsort(tmp, n, sizeof(float), &greaterthan);
01106                 float median;
01107                 if (n%2==1) median = tmp[n/2];
01108                 else median = (tmp[n/2-1]+tmp[n/2])/2.0f;
01109                 delete [] tmp;
01110                 return median;
01111         }
01112         else if (key == "nonzero_median")
01113         {
01114                 if ( is_complex() ) throw ImageFormatException("Error - can not calculate the median of a complex image");
01115                 vector<float> tmp;
01116                 size_t n = size;
01117                 float* d = get_data();
01118                 for( size_t i = 0; i < n; ++i ) {
01119                         if ( d[i] != 0 ) tmp.push_back(d[i]);
01120                 }
01121                 sort(tmp.begin(), tmp.end());
01122                 unsigned int vsize = tmp.size();
01123                 float median;
01124                 if (vsize%2==1) median = tmp[vsize/2];
01125                 else median = (tmp[vsize/2-1]+tmp[vsize/2])/2.0f;
01126                 return median;
01127         }
01128         else if (key == "changecount") return EMObject(changecount);
01129         else if (key == "nx") {
01130                 return nx;
01131         }
01132         else if (key == "ny") {
01133                 return ny;
01134         }
01135         else if (key == "nz") {
01136                 return nz;
01137         }
01138 
01139         if(attr_dict.has_key(key)) {
01140                 return attr_dict[key];
01141         }
01142         else {
01143                 throw NotExistingObjectException(key, "The requested key does not exist");
01144         }
01145 
01146         EXITFUNC;
01147 }
01148 
01149 EMObject EMData::get_attr_default(const string & key, const EMObject & em_obj) const
01150 {
01151         ENTERFUNC;
01152 
01153         if(attr_dict.has_key(key)) {
01154                 return get_attr(key);
01155         }
01156         else {
01157                 return em_obj;
01158         }
01159 
01160         EXITFUNC;
01161 }
01162 
01163 Dict EMData::get_attr_dict() const
01164 {
01165         if(rdata) {
01166                 update_stat();
01167         }
01168 
01169         Dict tmp=Dict(attr_dict);
01170         tmp["nx"]=nx;
01171         tmp["ny"]=ny;
01172         tmp["nz"]=nz;
01173         tmp["changecount"]=changecount;
01174 
01175         return tmp;
01176 }
01177 
01178 void EMData::set_attr_dict(const Dict & new_dict)
01179 {
01180         /*set nx, ny nz may resize the image*/
01181         // This wasn't supposed to 'clip' the image, but just redefine the size --steve
01182         if( new_dict.has_key("nx") || new_dict.has_key("ny") || new_dict.has_key("nz") ) {
01183                 LOGWARN("Warning: Ignored setting dimension size by modifying attribute!!!");
01184                 const_cast<Dict&>(new_dict).erase("nx");
01185                 const_cast<Dict&>(new_dict).erase("ny");
01186                 const_cast<Dict&>(new_dict).erase("nz");
01187         }
01188 
01189         vector<string> new_keys = new_dict.keys();
01190         vector<string>::const_iterator it;
01191         for(it = new_keys.begin(); it!=new_keys.end(); ++it) {
01192                 this->set_attr(*it, new_dict[*it]);
01193         }
01194 }
01195 
01196 void EMData::set_attr_dict_explicit(const Dict & new_dict)
01197 {
01198         attr_dict = new_dict;
01199 }
01200 
01201 void EMData::del_attr(const string & attr_name)
01202 {
01203         attr_dict.erase(attr_name);
01204 }
01205 
01206 void EMData::del_attr_dict(const vector<string> & del_keys)
01207 {
01208         vector<string>::const_iterator it;
01209         for(it=del_keys.begin(); it!=del_keys.end(); ++it) {
01210                 this->del_attr(*it);
01211         }
01212 }
01213 
01214 void EMData::set_attr(const string & key, EMObject val)
01215 {
01216         /* Ignore dimension attribute. */
01217         if(key == "nx" || key == "ny" || key == "nz")
01218         {
01219                 printf("Ignore setting dimension attribute %s. Use set_size if you need resize this EMData object.", key.c_str());
01220                 return;
01221         }
01222 
01223         if(rdata) {     //skip following for header only image
01224                 /* Ignore 'read only' attribute. */
01225                 if(key == "sigma" ||
01226                         key == "sigma_nonzero" ||
01227                         key == "square_sum" ||
01228                         key == "maximum" ||
01229                         key == "minimum" ||
01230                         key == "mean" ||
01231                         key == "mean_nonzero" )
01232                 {
01233                         LOGWARN("Ignore setting read only attribute %s", key.c_str());
01234                         return;
01235                 }
01236         }
01237 
01238         attr_dict[key] = val;
01239 }
01240 
01241 void EMData::set_attr_python(const string & key, EMObject val)
01242 {
01243         /* Ignore dimension attribute. */
01244         if(key == "nx" || key == "ny" || key == "nz")
01245         {
01246                 printf("Ignore setting dimension attribute %s. Use set_size if you need resize this EMData object.", key.c_str());
01247                 return;
01248         }
01249 
01250         /* Ignore 'read only' attribute. */
01251         if(key == "sigma" ||
01252                   key == "sigma_nonzero" ||
01253                   key == "square_sum" ||
01254                   key == "maximum" ||
01255                   key == "minimum" ||
01256                   key == "mean" ||
01257                   key == "mean_nonzero" )
01258         {
01259                 LOGWARN("Ignore setting read only attribute %s", key.c_str());
01260                 return;
01261         }
01262 
01263         EMObject::ObjectType argtype = val.get_type();
01264         if (argtype == EMObject::EMDATA) {
01265                 EMData* e = (EMData*) val;
01266                 e = e->copy();
01267                 EMObject v(e);
01268                 attr_dict[key] = v;
01269         }
01270         else if (argtype == EMObject::TRANSFORM) {
01271                 Transform* t = new Transform(*((Transform*) val));
01272                 EMObject v(t);
01273                 attr_dict[key] = v;
01274                 delete t; t=0;
01275         } else {
01276                 attr_dict[key] = val;
01277         }
01278 
01279 }
01280 
01281 void EMData::scale_pixel(float scale) const
01282 {
01283         attr_dict["apix_x"] = ((float) attr_dict["apix_x"]) * scale;
01284         attr_dict["apix_y"] = ((float) attr_dict["apix_y"]) * scale;
01285         attr_dict["apix_z"] = ((float) attr_dict["apix_z"]) * scale;
01286         if (attr_dict.has_key("ctf")) {
01287                 Ctf *ctf=(Ctf *)attr_dict["ctf"];
01288                 ctf->apix*=scale;
01289                 attr_dict["ctf"]=ctf;
01290                 if(ctf) {delete ctf; ctf=0;}
01291         }
01292 }
01293 
01294 
01295 //vector<float> EMData::get_data_pickle() const
01296 std::string EMData::get_data_pickle() const
01297 {
01298 //      vector<float> vf;
01299 //      vf.resize(nx*ny*nz);
01300 //      std::copy(rdata, rdata+nx*ny*nz, vf.begin());
01301 
01302         std::string vf((const char *)get_data(),nx*ny*nz*sizeof(float));
01303 
01304         return vf;
01305 }
01306 
01307 //void EMData::set_data_pickle(const vector<float>& vf)
01308 void EMData::set_data_pickle(std::string vf)
01309 {
01310 //      if (rdata) printf("rdata exists\n");
01311 //      rdata = (float *)malloc(nx*ny*nz*sizeof(float));
01312 //      std::copy(vf.begin(), vf.end(), rdata);
01313         EMUtil::em_memcpy(get_data(),vf.data(),(size_t)nx*ny*nz*sizeof(float));
01314 
01315 }
01316 
01317 int EMData::get_supp_pickle() const
01318 {
01319         return 0;
01320 }
01321 
01322 void EMData::set_supp_pickle(int)
01323 {
01324         this->supp = 0;
01325 }
01326 
01327 float EMData::get_amplitude_thres(float thres)
01328 {
01329 
01330         if (thres < 0 || thres > 1){
01331                 LOGERR("threshold bust be between 0 and 1.");
01332                 throw InvalidValueException(thres, "thres: 0 <= thres <= 1");
01333         }
01334                 
01335         EMData * amps = get_fft_amplitude();
01336         vector<float> ampvector = amps->get_data_as_vector();
01337         // yes I realize this may be slow if the map is big, but then again this function is only suited for tomo alignments, which if you have a big map will be VERY slow anyways!
01338         sort (ampvector.begin(), ampvector.end()); 
01339         int thresidx = int(thres * ampvector.size());
01340         float thresamp =  ampvector[thresidx];
01341 
01342         return thresamp;
01343 }
01344 
01345 vector<Vec3i > find_region(EMData* image,const vector<Vec3i >& coords, const float value, vector<Vec3i >& region)
01346 {
01347         static vector<Vec3i> two_six_connected;
01348         if (two_six_connected.size() == 0) {
01349                 for(int i = -1; i <= 1; ++i) {
01350                         for(int j = -1; j <= 1; ++j) {
01351                                 for(int  k = -1; k <= 1; ++k) {
01352                                         if ( j != 0 || i != 0 || k != 0) {
01353                                                 two_six_connected.push_back(Vec3i(i,j,k));
01354                                         }
01355                                 }
01356                         }
01357                 }
01358         }
01359 
01360         vector<Vec3i> ret;
01361         for(vector<Vec3i>::const_iterator it = two_six_connected.begin(); it != two_six_connected.end(); ++it ) {
01362                 for(vector<Vec3i>::const_iterator it2 = coords.begin(); it2 != coords.end(); ++it2 ) {
01363                         if  (image->get_value_at((*it2)[0],(*it2)[1],(*it2)[2]) != value) throw;
01364                         Vec3i c = (*it)+(*it2);
01365 
01366                         if ( c[0] < 0 || c[0] >= image->get_xsize()) continue;
01367                         if ( c[1] < 0 || c[1] >= image->get_ysize()) continue;
01368                         if ( c[2] < 0 || c[2] >= image->get_zsize()) continue;
01369 
01370                         if( image->get_value_at(c[0],c[1],c[2]) == value ) {
01371                                 if (find(ret.begin(),ret.end(),c) == ret.end()) {
01372                                         if (find(region.begin(),region.end(),c) == region.end()) {
01373                                                 region.push_back(c);
01374                                                 ret.push_back(c);
01375                                         }
01376                                 }
01377                         }
01378                 }
01379         }
01380         return ret;
01381 }
01382 
01383 vector<Vec3i> EMData::mask_contig_region(const float& value, const Vec3i& seed) {
01384         Vec3i coord(seed[0],seed[1],seed[2]);
01385         vector<Vec3i> region;
01386         region.push_back(coord);
01387         vector<Vec3i> find_region_input = region;
01388         while (true) {
01389                 vector<Vec3i> v = find_region(this,find_region_input, value, region);
01390                 if (v.size() == 0 ) break;
01391                 else find_region_input = v;
01392         }
01393         return region;
01394 }

Generated on Tue Jun 11 13:40:37 2013 for EMAN2 by  doxygen 1.3.9.1