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(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 FloatPoint EMData::calc_center_of_mass(float threshold)
00450 {
00451         float *data = get_data();
00452 
00453         //float sigma = get_attr("sigma");
00454         //float mean = get_attr("mean");
00455         float m = 0.0;
00456 
00457         FloatPoint com(0,0,0);
00458         for (int i = 0; i < nx; ++i) {
00459                 for (int j = 0; j < ny; ++j) {
00460                         int j2 = nx * j;
00461                         for (int k = 0; k < nz; ++k) {
00462                                 size_t l = i + j2 + (size_t)k * nxy;
00463                                 if (data[l] >= threshold) {             // threshold out noise (and negative density)
00464                                         com[0] += i * data[l];
00465                                         com[1] += j * data[l];
00466                                         com[2] += k * data[l];
00467                                         m += data[l];
00468                                 }
00469                         }
00470                 }
00471         }
00472 
00473         com[0] /= m;
00474         com[1] /= m;
00475         com[2] /= m;
00476 
00477         return com;
00478 }
00479 
00480 
00481 size_t EMData::calc_min_index() const
00482 {
00483         IntPoint min_location = calc_min_location();
00484         size_t i = min_location[0] + min_location[1] * nx + (size_t)min_location[2] * nx * ny;
00485         return i;
00486 }
00487 
00488 
00489 size_t EMData::calc_max_index() const
00490 {
00491         IntPoint max_location = calc_max_location();
00492         size_t i = max_location[0] + max_location[1] * nx + (size_t)max_location[2] * nx * ny;
00493         return i;
00494 }
00495 
00496 
00497 vector<Pixel> EMData::calc_highest_locations(float threshold) const 
00498 {
00499         ENTERFUNC;
00500 
00501         vector<Pixel> result;
00502 
00503         int di = 1;
00504         if (is_complex() && !is_ri()) {
00505                 di = 2;
00506         }
00507 
00508         int nxy = nx * ny;
00509         float * data = get_data();
00510 
00511         for (int j = 0; j < nz; ++j) {
00512                 size_t cur_z = (size_t)j * nxy;
00513 
00514                 for (int k = 0; k < ny; ++k) {
00515                         size_t cur_y = k * nx + cur_z;
00516 
00517                         for (int l = 0; l < nx; l += di) {
00518                                 float v =data[l + cur_y];
00519                                 if (v > threshold) {
00520                                         result.push_back(Pixel(l, k, j, v));
00521                                 }
00522                         }
00523                 }
00524         }
00525 
00526         std::sort(result.begin(), result.end());
00527 
00528         EXITFUNC;
00529         return result;
00530 }
00531 
00532 vector<Pixel> EMData::calc_n_highest_locations(int n)
00533 {
00534         ENTERFUNC;
00535 
00536         vector<Pixel> result;
00537 
00538         int di = 1;
00539         if (is_complex() && !is_ri()) {
00540                 di = 2;
00541         }
00542 
00543         // initialize with n elements
00544         float * data = get_data();
00545         for ( int i=0; i<n; i++) result.push_back(Pixel(0,0,0,data[0]));
00546 
00547         int nxy = nx * ny;
00548 
00549         for (int j = 0; j < nz; ++j) {
00550                 size_t cur_z = (size_t)j * nxy;
00551 
00552                 for (int k = 0; k < ny; ++k) {
00553                         size_t cur_y = k * nx + cur_z;
00554 
00555                         for (int l = 0; l < nx; l += di) {
00556                                 float v =data[l + cur_y];
00557                                 if (v<result[n-1].value) continue;
00558                                 for (vector<Pixel>::iterator i=result.begin(); i<result.end(); i++) {
00559                                         if (v>(*i).value) { result.insert(i,Pixel(l, k, j, v)); result.pop_back(); break; }
00560                                 }
00561                         }
00562                 }
00563         }
00564 
00565         EXITFUNC;
00566         return result;
00567 }
00568 
00569 vector<Pixel> EMData::find_pixels_with_value(float val) 
00570 {
00571         ENTERFUNC;
00572         
00573         if ( is_complex() ) throw ImageFormatException("Error - find_pixels_with_value real only");
00574 
00575         vector<Pixel> result;
00576 
00577         for (int k = 0; k < nz; k++) {
00578                 for (int j = 0; j < ny; j++) {
00579                         for (int i = 0; i < nx; i++) {
00580                                 if (get_value_at(i,j,k)==val) result.push_back(Pixel(i,j,k,val));
00581                         }
00582                 }
00583         }
00584 
00585         EXITFUNC;
00586         return result;
00587 }
00588 
00589 float EMData::get_edge_mean() const
00590 {
00591         ENTERFUNC;
00592 #ifdef EMAN2_USING_CUDA
00593         if(cudarwdata){
00594                 
00595                         return get_edgemean_cuda(cudarwdata, nx, ny, nz);
00596                         
00597         }
00598 #endif
00599         int di = 0;
00600         double edge_sum = 0;
00601         float edge_mean = 0;
00602         size_t nxy = nx * ny;
00603         float * data = get_data();
00604         if (nz == 1) {
00605                 for (int i = 0, j = (ny - 1) * nx; i < nx; ++i, ++j) {
00606                         edge_sum += data[i] + data[j];
00607                 }
00608                 for (size_t i = 0, j = nx - 1; i < nxy; i += nx, j += nx) {
00609                         edge_sum += data[i] + data[j];
00610                 }
00611                 edge_mean = (float)edge_sum / (nx * 2 + ny * 2);
00612         }
00613         else {
00614                 if (nx == ny && nx == nz * 2 - 1) {
00615                         for (size_t j = (nxy * (nz - 1)); j < nxy * nz; ++j, ++di) {
00616                                 edge_sum += data[j];
00617                         }
00618                 }
00619                 else {
00620                         for (size_t i = 0, j = (nxy * (nz - 1)); i < nxy; ++i, ++j, ++di) {
00621                                 edge_sum += data[i] + data[j];
00622                         }
00623                 }
00624 
00625                 int nxy2 = nx * (ny - 1);
00626                 for (int k = 1; k < nz - 1; ++k) {
00627                         size_t k2 = k * nxy;
00628                         size_t k3 = k2 + nxy2;
00629                         for (int i = 0; i < nx; ++i, ++di) {
00630                                 edge_sum += data[i + k2] + data[i + k3];
00631                         }
00632                 }
00633                 for (int k = 1; k < nz - 1; ++k) {
00634                         size_t k2 = k * nxy;
00635                         size_t k3 = nx - 1 + k2;
00636                         for (int i = 1; i < ny - 1; ++i, ++di) {
00637                                 edge_sum += data[i * nx + k2] + data[i * nx + k3];
00638                         }
00639                 }
00640 
00641                 edge_mean = (float)edge_sum / (di * 2);
00642         }
00643         EXITFUNC;
00644 
00645         return edge_mean;
00646 }
00647 
00648 
00649 float EMData::get_circle_mean()
00650 {
00651         ENTERFUNC;
00652 
00653         static bool busy = false;
00654         static EMData *mask = 0;
00655 
00656         while (busy);
00657         busy = true;
00658 
00659         if (!mask || !EMUtil::is_same_size(this, mask)) {
00660                 if (!mask) {
00661                         mask = new EMData();
00662                 }
00663                 mask->set_size(nx, ny, nz);
00664                 mask->to_one();
00665 
00666                 float radius = (float)(ny / 2 - 2);
00667                 mask->process_inplace("mask.sharp", Dict("inner_radius", radius - 1,
00668                                                                            "outer_radius", radius + 1));
00669 
00670         }
00671         double n = 0,s=0;
00672         float *d = mask->get_data();
00673         float * data = get_data();
00674         size_t size = (size_t)nx*ny*nz;
00675         for (size_t i = 0; i < size; ++i) {
00676                 if (d[i]) { n+=1.0; s+=data[i]; }
00677         }
00678 
00679 
00680         float result = (float)(s/n);
00681         busy = false;
00682 
00683         EXITFUNC;
00684         return result;
00685 }
00686 
00687 
00688 void EMData::set_ctf(Ctf * new_ctf)
00689 {
00690         ENTERFUNC;
00691 
00692         vector<float> vctf = new_ctf->to_vector();
00693         attr_dict["ctf"] = vctf;
00694 
00695         EXITFUNC;
00696 }
00697 
00698 Ctf * EMData::get_ctf() const
00699 {
00700         if(attr_dict.has_key("ctf")) {
00701                 EMAN1Ctf * ctf = new EMAN1Ctf();
00702                 ctf->from_vector(attr_dict["ctf"]);
00703 
00704                 return dynamic_cast<Ctf *>(ctf);
00705         }
00706         else {
00707                 return 0;
00708         }
00709 }
00710 
00711 #include <iostream>
00712 using std::cout;
00713 using std::endl;
00714 
00715 void EMData::set_size(int x, int y, int z)
00716 {
00717         ENTERFUNC;
00718 
00719         if (x <= 0) {
00720                 throw InvalidValueException(x, "x size <= 0");
00721         }
00722         else if (y <= 0) {
00723                 throw InvalidValueException(y, "y size <= 0");
00724         }
00725         else if (z <= 0) {
00726                 throw InvalidValueException(z, "z size <= 0");
00727         }
00728 
00729 #ifdef MEMDEBUG2
00730         printf("EMDATA sz  %4d    %p (%d,%d,%d)\n",EMData::totalalloc,this,x,y,z);
00731 #endif
00732 
00733 
00734         int old_nx = nx;
00735 
00736         size_t size = x*y*z*sizeof(float);
00737 
00738         if (rdata != 0) {
00739                 rdata = (float*)EMUtil::em_realloc(rdata,size);
00740         } else {
00741                 // Just pass on this for a while....see what happens
00742                 rdata = (float*)EMUtil::em_malloc(size);
00743         }
00744 //      rdata = static_cast < float *>(realloc(rdata, size));
00745         if ( rdata == 0 )
00746         {
00747                 stringstream ss;
00748                 string gigs;
00749                 ss << (float) size/1000000000.0;
00750                 ss >> gigs;
00751                 string message = "Cannot allocate " + gigs + " GB - not enough memory.";
00752                 throw BadAllocException(message);
00753         }
00754 
00755         nx = x;
00756         ny = y;
00757         nz = z;
00758         nxy = nx*ny;
00759         nxyz = (size_t)nx*ny*nz;
00760 
00761 // once the object is resized, the CUDA need to be updated
00762 #ifdef EMAN2_USING_CUDA
00763         if(cudarwdata) {
00764                 rw_free();
00765                 rw_alloc();
00766         }
00767         if(cudarodata) {
00768                 ro_free();
00769                 ro_alloc();
00770         }
00771 #endif // EMAN2_USING_CUDA
00772 
00773         if (old_nx == 0) {
00774                 EMUtil::em_memset(get_data(),0,size);
00775         }
00776 
00777         if (supp) {
00778                 EMUtil::em_free(supp);
00779                 supp = 0;
00780         }
00781 
00782         update();
00783         EXITFUNC;
00784 }
00785 
00786 #ifdef EMAN2_USING_CUDA
00787 
00788 void EMData::set_size_cuda(int x, int y, int z)
00789 {
00790         ENTERFUNC;
00791 
00792         if (x <= 0) {
00793                 throw InvalidValueException(x, "x size <= 0");
00794         }
00795         else if (y <= 0) {
00796                 throw InvalidValueException(y, "y size <= 0");
00797         }
00798         else if (z <= 0) {
00799                 throw InvalidValueException(z, "z size <= 0");
00800         }
00801 
00802         nx = x;
00803         ny = y;
00804         nz = z;
00805 
00806         nxy = nx*ny;
00807         nxyz = (size_t)nx*ny*nz;
00808 
00809 //      gpu_update();
00810 
00811         EXITFUNC;
00812 }
00813 
00814 #endif // EMAN2_USING_CUDA
00815 
00816 MArray2D EMData::get_2dview() const
00817 {
00818         const int ndims = 2;
00819         if (get_ndim() != ndims) {
00820                 throw ImageDimensionException("2D only");
00821         }
00822         boost::array<std::size_t,ndims> dims = {{nx, ny}};
00823         MArray2D marray(get_data(), dims, boost::fortran_storage_order());
00824         return marray;
00825 }
00826 
00827 
00828 MArray3D EMData::get_3dview() const
00829 {
00830         const int ndims = 3;
00831         boost::array<std::size_t,ndims> dims = {{nx, ny, nz}};
00832         MArray3D marray(get_data(), dims, boost::fortran_storage_order());
00833         return marray;
00834 }
00835 
00836 
00837 MCArray2D EMData::get_2dcview() const
00838 {
00839         const int ndims = 2;
00840         if (get_ndim() != ndims) {
00841                 throw ImageDimensionException("2D only");
00842         }
00843         boost::array<std::size_t,ndims> dims = {{nx/2, ny}};
00844         std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
00845         MCArray2D marray(cdata, dims, boost::fortran_storage_order());
00846         return marray;
00847 }
00848 
00849 
00850 MCArray3D EMData::get_3dcview() const
00851 {
00852         const int ndims = 3;
00853         boost::array<std::size_t,ndims> dims = {{nx/2, ny, nz}};
00854         std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
00855         MCArray3D marray(cdata, dims, boost::fortran_storage_order());
00856         return marray;
00857 }
00858 
00859 
00860 MCArray3D* EMData::get_3dcviewptr() const
00861 {
00862         const int ndims = 3;
00863         boost::array<std::size_t,ndims> dims = {{nx/2, ny, nz}};
00864         std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
00865         MCArray3D* marray = new MCArray3D(cdata, dims,
00866                                                                           boost::fortran_storage_order());
00867         return marray;
00868 }
00869 
00870 
00871 MArray2D EMData::get_2dview(int x0, int y0) const
00872 {
00873         const int ndims = 2;
00874         if (get_ndim() != ndims) {
00875                 throw ImageDimensionException("2D only");
00876         }
00877         boost::array<std::size_t,ndims> dims = {{nx, ny}};
00878         MArray2D marray(get_data(), dims, boost::fortran_storage_order());
00879         boost::array<std::size_t,ndims> bases={{x0, y0}};
00880         marray.reindex(bases);
00881         return marray;
00882 }
00883 
00884 
00885 MArray3D EMData::get_3dview(int x0, int y0, int z0) const
00886 {
00887         const int ndims = 3;
00888         boost::array<std::size_t,ndims> dims = {{nx, ny, nz}};
00889         MArray3D marray(get_data(), dims, boost::fortran_storage_order());
00890         boost::array<std::size_t,ndims> bases={{x0, y0, z0}};
00891         marray.reindex(bases);
00892         return marray;
00893 }
00894 
00895 
00896 MCArray2D EMData::get_2dcview(int x0, int y0) const
00897 {
00898         const int ndims = 2;
00899         if (get_ndim() != ndims) {
00900                 throw ImageDimensionException("2D only");
00901         }
00902         boost::array<std::size_t,ndims> dims = {{nx/2, ny}};
00903         std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
00904         MCArray2D marray(cdata, dims, boost::fortran_storage_order());
00905         boost::array<std::size_t,ndims> bases={{x0, y0}};
00906         marray.reindex(bases);
00907         return marray;
00908 }
00909 
00910 
00911 MCArray3D EMData::get_3dcview(int x0, int y0, int z0) const
00912 {
00913         const int ndims = 3;
00914         boost::array<std::size_t,ndims> dims = {{nx/2, ny, nz}};
00915         std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
00916         MCArray3D marray(cdata, dims, boost::fortran_storage_order());
00917         boost::array<std::size_t,ndims> bases={{x0, y0, z0}};
00918         marray.reindex(bases);
00919         return marray;
00920 }
00921 
00922 int greaterthan( const void* p1, const void* p2 )
00923 {
00924         float*  v1 = (float*) p1;
00925         float*  v2 = (float*) p2;
00926 
00927         if ( *v1 < *v2 )  return 0;
00928         else return 1;
00929 }
00930 
00931 
00932 EMObject EMData::get_attr(const string & key) const
00933 {
00934         ENTERFUNC;
00935         
00936         if ((flags & EMDATA_NEEDUPD) && (key != "is_fftpad") && (key != "xform.align2d")){update_stat();} //this gives a spped up of 7.3% according to e2speedtest
00937                 
00938         size_t size = (size_t)nx * ny * nz;
00939         if (key == "kurtosis") {
00940                 float mean = attr_dict["mean"];
00941                 float sigma = attr_dict["sigma"];
00942 
00943                 float *data = get_data();
00944                 double kurtosis_sum = 0;
00945 
00946                 for (size_t k = 0; k < size; ++k) {
00947                         float t = (data[k] - mean) / sigma;
00948                         float tt = t * t;
00949                         kurtosis_sum += tt * tt;
00950                 }
00951 
00952                 float kurtosis = (float)(kurtosis_sum / size - 3.0);
00953                 return kurtosis;
00954         }
00955         else if (key == "skewness") {
00956                 float mean = attr_dict["mean"];
00957                 float sigma = attr_dict["sigma"];
00958 
00959                 float *data = get_data();
00960                 double skewness_sum = 0;
00961                 for (size_t k = 0; k < size; ++k) {
00962                         float t = (data[k] - mean) / sigma;
00963                         skewness_sum +=  t * t * t;
00964                 }
00965                 float skewness = (float)(skewness_sum / size);
00966                 return skewness;
00967         }
00968         else if (key == "median")
00969         {
00970                 if ( is_complex() ) throw ImageFormatException("Error - can not calculate the median of a complex image");
00971                 size_t n = size;
00972                 float* tmp = new float[n];
00973                 float* d = get_data();
00974                 if (tmp == 0 ) throw BadAllocException("Error - could not create deep copy of image data");
00975 //              for(size_t i=0; i < n; ++i) tmp[i] = d[i]; // should just be a memcpy
00976                 std::copy(d, d+n, tmp);
00977                 qsort(tmp, n, sizeof(float), &greaterthan);
00978                 float median;
00979                 if (n%2==1) median = tmp[n/2];
00980                 else median = (tmp[n/2-1]+tmp[n/2])/2.0f;
00981                 delete [] tmp;
00982                 return median;
00983         }
00984         else if (key == "nonzero_median")
00985         {
00986                 if ( is_complex() ) throw ImageFormatException("Error - can not calculate the median of a complex image");
00987                 vector<float> tmp;
00988                 size_t n = size;
00989                 float* d = get_data();
00990                 for( size_t i = 0; i < n; ++i ) {
00991                         if ( d[i] != 0 ) tmp.push_back(d[i]);
00992                 }
00993                 sort(tmp.begin(), tmp.end());
00994                 unsigned int vsize = tmp.size();
00995                 float median;
00996                 if (vsize%2==1) median = tmp[vsize/2];
00997                 else median = (tmp[vsize/2-1]+tmp[vsize/2])/2.0f;
00998                 return median;
00999         }
01000         else if (key == "changecount") return EMObject(changecount);
01001         else if (key == "nx") return nx;
01002         else if (key == "ny") return ny;
01003         else if (key == "nz") return nz;
01004 
01005         if(attr_dict.has_key(key)) {
01006                 return attr_dict[key];
01007         }
01008         else {
01009                 throw NotExistingObjectException(key, "The requested key does not exist");
01010         }
01011 
01012         EXITFUNC;
01013 }
01014 
01015 EMObject EMData::get_attr_default(const string & key, const EMObject & em_obj) const
01016 {
01017         ENTERFUNC;
01018 
01019         if(attr_dict.has_key(key)) {
01020                 return get_attr(key);
01021         }
01022         else {
01023                 return em_obj;
01024         }
01025 
01026         EXITFUNC;
01027 }
01028 
01029 Dict EMData::get_attr_dict() const
01030 {
01031         update_stat();
01032         
01033         Dict tmp=Dict(attr_dict);
01034         tmp["nx"]=nx;
01035         tmp["ny"]=ny;
01036         tmp["nz"]=nz;
01037         tmp["changecount"]=changecount;
01038 
01039         return tmp;
01040 }
01041 
01042 void EMData::set_attr_dict(const Dict & new_dict)
01043 {
01044         /*set nx, ny nz may resize the image*/
01045         // This wasn't supposed to 'clip' the image, but just redefine the size --steve
01046         if( ( new_dict.has_key("nx") && nx!=(int)new_dict["nx"] )
01047                 || ( new_dict.has_key("ny") && ny!=(int)new_dict["ny"] )
01048                 || ( new_dict.has_key("nz") && nz!=(int)new_dict["nz"] ) ) {
01049 
01050                 int newx, newy, newz;
01051                 newx = new_dict.has_key("nx") ? (int)new_dict["nx"] : nx;
01052                 newy = new_dict.has_key("ny") ? (int)new_dict["ny"] : ny;
01053                 newz = new_dict.has_key("nz") ? (int)new_dict["nz"] : nz;
01054 
01055                 set_size(newx,newy,newz);
01056 
01057 //              EMData * new_image = get_clip(Region((nx-newx)/2, (ny-newy)/2, (nz=newz)/2, newx, newy, newz));
01058 //              if(new_image) {
01059 //                      this->operator=(*new_image);
01060 //                      delete new_image;
01061 //                      new_image = 0;
01062 //              }
01063         }
01064 
01065         vector<string> new_keys = new_dict.keys();
01066         vector<string>::const_iterator it;
01067         for(it = new_keys.begin(); it!=new_keys.end(); ++it) {
01068                 this->set_attr(*it, new_dict[*it]);
01069         }
01070 }
01071 
01072 void EMData::set_attr_dict_explicit(const Dict & new_dict)
01073 {
01074         attr_dict = new_dict;
01075 }
01076 
01077 void EMData::del_attr(const string & attr_name)
01078 {
01079         attr_dict.erase(attr_name);
01080 }
01081 
01082 void EMData::del_attr_dict(const vector<string> & del_keys)
01083 {
01084         vector<string>::const_iterator it;
01085         for(it=del_keys.begin(); it!=del_keys.end(); ++it) {
01086                 this->del_attr(*it);
01087         }
01088 }
01089 
01090 void EMData::set_attr(const string & key, EMObject val)
01091 {
01092         if( key == "nx" && nx != (int)val) { set_size((int)val,ny,nz); return; }
01093         if( key == "ny" && ny != (int)val) { set_size(nx,(int)val,nz); return; }
01094         if( key == "nz" && nz != (int)val) { set_size(nx,ny,(int)val); return; }
01095 
01096         /* Ignore 'read only' attribute. */
01097         if(key == "sigma" ||
01098                 key == "sigma_nonzero" ||
01099                 key == "square_sum" ||
01100                 key == "maximum" ||
01101                 key == "minimum" ||
01102                 key == "mean" ||
01103                 key == "mean_nonzero" )
01104         {
01105                 LOGWARN("Ignore setting read only attribute %s", key.c_str());
01106                 return;
01107         }
01108 
01109         attr_dict[key] = val;
01110 
01111 
01112 
01113 }
01114 
01115 void EMData::set_attr_python(const string & key, EMObject val)
01116 {
01117         if( key == "nx" && nx != (int)val) { set_size((int)val,ny,nz); return; }
01118         if( key == "ny" && ny != (int)val) { set_size(nx,(int)val,nz); return; }
01119         if( key == "nz" && nz != (int)val) { set_size(nx,ny,(int)val); return; }
01120 
01121         /* Ignore 'read only' attribute. */
01122         if(key == "sigma" ||
01123                   key == "sigma_nonzero" ||
01124                   key == "square_sum" ||
01125                   key == "maximum" ||
01126                   key == "minimum" ||
01127                   key == "mean" ||
01128                   key == "mean_nonzero" )
01129         {
01130                 LOGWARN("Ignore setting read only attribute %s", key.c_str());
01131                 return;
01132         }
01133 
01134         EMObject::ObjectType argtype = val.get_type();
01135         if (argtype == EMObject::EMDATA) {
01136                 EMData* e = (EMData*) val;
01137                 e = e->copy();
01138                 EMObject v(e);
01139                 attr_dict[key] = v;
01140         }
01141         else if (argtype == EMObject::TRANSFORM) {
01142                 Transform* t = new Transform(*((Transform*) val));
01143                 EMObject v(t);
01144                 attr_dict[key] = v;
01145                 delete t; t=0;
01146         } else {
01147                 attr_dict[key] = val;
01148         }
01149 
01150 }
01151 
01152 void EMData::scale_pixel(float scale) const
01153 {
01154         attr_dict["apix_x"] = ((float) attr_dict["apix_x"]) * scale;
01155         attr_dict["apix_y"] = ((float) attr_dict["apix_y"]) * scale;
01156         attr_dict["apix_z"] = ((float) attr_dict["apix_z"]) * scale;
01157         if (attr_dict.has_key("ctf")) {
01158                 Ctf *ctf=(Ctf *)attr_dict["ctf"];
01159                 ctf->apix*=scale;
01160                 attr_dict["ctf"]=ctf;
01161                 if(ctf) {delete ctf; ctf=0;}
01162         }
01163 }
01164 
01165 
01166 //vector<float> EMData::get_data_pickle() const
01167 std::string EMData::get_data_pickle() const
01168 {
01169 //      vector<float> vf;
01170 //      vf.resize(nx*ny*nz);
01171 //      std::copy(rdata, rdata+nx*ny*nz, vf.begin());
01172 
01173         std::string vf((const char *)get_data(),nx*ny*nz*sizeof(float));
01174 
01175         return vf;
01176 }
01177 
01178 //void EMData::set_data_pickle(const vector<float>& vf)
01179 void EMData::set_data_pickle(std::string vf)
01180 {
01181 //      if (rdata) printf("rdata exists\n");
01182 //      rdata = (float *)malloc(nx*ny*nz*sizeof(float));
01183 //      std::copy(vf.begin(), vf.end(), rdata);
01184         EMUtil::em_memcpy(get_data(),vf.data(),nx*ny*nz*sizeof(float));
01185 
01186 }
01187 
01188 int EMData::get_supp_pickle() const
01189 {
01190         return 0;
01191 }
01192 
01193 void EMData::set_supp_pickle(int)
01194 {
01195         this->supp = 0;
01196 }
01197 
01198 float EMData::get_amplitude_thres(float thres)
01199 {
01200 
01201         if (thres < 0 || thres > 1){
01202                 LOGERR("threshold bust be between 0 and 1.");
01203                 throw InvalidValueException(thres, "thres: 0 <= thres <= 1");
01204         }
01205                 
01206         EMData * amps = get_fft_amplitude();
01207         vector<float> ampvector = amps->get_data_as_vector();
01208         // 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!
01209         sort (ampvector.begin(), ampvector.end()); 
01210         int thresidx = int(thres * ampvector.size());
01211         float thresamp =  ampvector[thresidx];
01212 
01213         return thresamp;
01214 }
01215 
01216 vector<Vec3i > find_region(EMData* image,const vector<Vec3i >& coords, const float value, vector<Vec3i >& region)
01217 {
01218         static vector<Vec3i> two_six_connected;
01219         if (two_six_connected.size() == 0) {
01220                 for(int i = -1; i <= 1; ++i) {
01221                         for(int j = -1; j <= 1; ++j) {
01222                                 for(int  k = -1; k <= 1; ++k) {
01223                                         if ( j != 0 || i != 0 || k != 0) {
01224                                                 two_six_connected.push_back(Vec3i(i,j,k));
01225                                         }
01226                                 }
01227                         }
01228                 }
01229         }
01230 
01231         vector<Vec3i> ret;
01232         for(vector<Vec3i>::const_iterator it = two_six_connected.begin(); it != two_six_connected.end(); ++it ) {
01233                 for(vector<Vec3i>::const_iterator it2 = coords.begin(); it2 != coords.end(); ++it2 ) {
01234                         if  (image->get_value_at((*it2)[0],(*it2)[1],(*it2)[2]) != value) throw;
01235                         Vec3i c = (*it)+(*it2);
01236 
01237                         if ( c[0] < 0 || c[0] >= image->get_xsize()) continue;
01238                         if ( c[1] < 0 || c[1] >= image->get_ysize()) continue;
01239                         if ( c[2] < 0 || c[2] >= image->get_zsize()) continue;
01240 
01241                         if( image->get_value_at(c[0],c[1],c[2]) == value ) {
01242                                 if (find(ret.begin(),ret.end(),c) == ret.end()) {
01243                                         if (find(region.begin(),region.end(),c) == region.end()) {
01244                                                 region.push_back(c);
01245                                                 ret.push_back(c);
01246                                         }
01247                                 }
01248                         }
01249                 }
01250         }
01251         return ret;
01252 }
01253 
01254 vector<Vec3i> EMData::mask_contig_region(const float& value, const Vec3i& seed) {
01255         Vec3i coord(seed[0],seed[1],seed[2]);
01256         vector<Vec3i> region;
01257         region.push_back(coord);
01258         vector<Vec3i> find_region_input = region;
01259         while (true) {
01260                 vector<Vec3i> v = find_region(this,find_region_input, value, region);
01261                 if (v.size() == 0 ) break;
01262                 else find_region_input = v;
01263         }
01264         return region;
01265 }

Generated on Mon Mar 7 18:18:30 2011 for EMAN2 by  doxygen 1.3.9.1