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