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 
00048 
00049 using namespace EMAN;
00050 
00051 EMData* EMData::get_fft_amplitude2D()
00052 {
00053         ENTERFUNC;
00054 
00055 //      int ndim = get_ndim();
00056         if (!is_complex()) {
00057                 LOGERR("complex image expected. Input image is real image.");
00058                 throw ImageFormatException("complex image expected. Input image is a real image.");
00059         }
00060         if (nz>1) {
00061                 LOGERR("2D image expected. Input image is 3D");
00062                 throw ImageFormatException("2D odd square complex image"
00063                         " expected Input image is 3D.");
00064         }
00065 
00066         int nx2 = nx/2;
00067 
00068         EMData *dat = copy_head();
00069 
00070         dat->set_size(nx2, ny, nz);
00071         dat->to_zero();
00072 
00073         float temp=0;
00074 
00075         for (int j = 0; j < ny; j++) {
00076                 for (int i = 0; i < nx2; i++) {
00077                         temp = (*this)(2*i,j)*(*this)(2*i,j);
00078                         temp += (*this)(2*i+1,j)*(*this)(2*i+1,j);
00079                         (*dat)(i,j) = std::sqrt(temp);
00080                 }
00081         }
00082 
00083         dat->update();
00084         dat->set_complex(false);
00085         dat->set_ri(false);
00086 
00087         EXITFUNC;
00088         return dat;
00089 }
00090 
00091 
00092 EMData* EMData::get_fft_amplitude()
00093 {
00094         ENTERFUNC;
00095 
00096         if (!is_complex()) {
00097                 LOGERR("complex image expected. Input image is real image.");
00098                 throw ImageFormatException("complex image expected. Input image is a real image.");
00099         }
00100 
00101         ri2ap();
00102 
00103         int nx2 = nx - 2;
00104         EMData *dat = copy_head();
00105         dat->set_size(nx2, ny, nz);
00106         dat->to_zero();
00107 
00108         float *d = dat->get_data();
00109         float *data = get_data();
00110         int ndim = get_ndim();
00111 
00112         size_t idx1, idx2, idx3;
00113         if (ndim == 3) {
00114                 for (int k = 1; k < nz; ++k) {
00115                         for (int j = 1; j < ny; ++j) {
00116                                 for (int i = 0; i < nx2/2; ++i) {
00117                                         idx1 = k*nx2*ny+j*nx2+nx2/2+i;
00118                                         idx2 = k*nx*ny+j*nx+2*i;
00119                                         idx3 = (nz-k)*nx2*ny+(ny-j)*nx2+nx2/2-i;
00120                                         d[idx1] = data[idx2];
00121                                         d[idx3] = data[idx2];
00122                                 }
00123                         }
00124                 }
00125         }
00126         else if (ndim == 2) {
00127                 for (int j = 1; j < ny; ++j) {
00128                         for (int i = 0; i < nx2/2; ++i) {
00129                                 d[j*nx2+nx2/2+i] = data[j*nx+2*i];
00130                                 d[(ny-j)*nx2+nx2/2-i] = data[j*nx+2*i];
00131                         }
00132                 }
00133         }
00134 
00135         dat->update();
00136         dat->set_complex(false);
00137         if(dat->get_ysize()==1 && dat->get_zsize()==1) {
00138                 dat->set_complex_x(false);
00139         }
00140         dat->set_ri(false);
00141 
00142         EXITFUNC;
00143         return dat;
00144 }
00145 
00146 
00147 EMData* EMData::get_fft_phase()
00148 {
00149         ENTERFUNC;
00150 
00151         if (!is_complex()) {
00152                 LOGERR("complex image expected. Input image is real image.");
00153                 throw ImageFormatException("complex image expected. Input image is a real image.");
00154         }
00155 
00156         ri2ap();
00157 
00158         int nx2 = nx - 2;
00159         EMData *dat = copy_head();
00160         dat->set_size(nx2, ny, nz);
00161         dat->to_zero();
00162 
00163         float *d = dat->get_data();
00164         float * data = get_data();
00165 
00166         int ndim = get_ndim();
00167         size_t idx1, idx2, idx3;
00168         if (ndim == 3) {
00169                 for (int k = 1; k < nz; ++k) {
00170                         for (int j = 1; j < ny; ++j) {
00171                                 for (int i = 0; i < nx2/2; ++i) {
00172                                         idx1 = k*nx2*ny+j*nx2+nx2/2+i;
00173                                         idx2 = k*nx*ny+j*nx+2*i+1;
00174                                         idx3 = (nz-k)*nx2*ny+(ny-j)*nx2+nx2/2-i;
00175                                         d[idx1] = data[idx2];
00176                                         d[idx3] = -data[idx2];
00177                                 }
00178                         }
00179                 }
00180         }
00181         else if (ndim == 2) {
00182                 for (int j = 1; j < ny; ++j) {
00183                         for (int i = 0; i < nx2/2; ++i) {
00184                                 d[j*nx2+nx2/2+i] = data[j*nx+2*i+1];
00185                                 d[(ny-j)*nx2+nx2/2-i] = -data[j*nx+2*i+1];
00186                         }
00187                 }
00188         }
00189 
00190         dat->update();
00191         dat->set_complex(false);
00192         if(dat->get_ysize()==1 && dat->get_zsize()==1) {
00193                 dat->set_complex_x(false);
00194         }
00195         dat->set_ri(false);
00196 
00197         EXITFUNC;
00198         return dat;
00199 }
00200 #ifdef EMAN2_USING_CUDA
00201 #include <cuda_runtime_api.h>
00202 
00203 float* EMData::get_data() const
00204 {
00205 // This is inappropriate, should ONLY be in set_size. Was added by David for CUDA work
00206 // I'm taking it back out again  --steve 11/03/2009
00207 /*      if (num_bytes > 0 && rdata == 0) {
00208                 rdata = (float*)EMUtil::em_malloc(num_bytes);
00209                 if ( rdata == 0 )
00210                 {
00211                         stringstream ss;
00212                         string gigs;
00213                         ss << num_bytes/1000000000.0;
00214                         ss >> gigs;
00215                         string message = "Cannot allocate " + gigs + " GB - not enough memory.";
00216                         throw BadAllocException(message);
00217                 }
00218                 if (rdata == 0) throw BadAllocException("The allocation of the raw data failed");
00219         }*/
00220 #ifdef EMAN2_USING_CUDA
00221         size_t num_bytes = nx*ny*nz*sizeof(float);
00222         if ( num_bytes > 0 && gpu_rw_is_current()  && (EMDATA_CPU_NEEDS_UPDATE & flags)) {
00223                 cudaError_t error = cudaMemcpy(rdata,get_cuda_data(),num_bytes,cudaMemcpyDeviceToHost);
00224                 if (error != cudaSuccess ) throw UnexpectedBehaviorException("The device to host cudaMemcpy failed : " + string(cudaGetErrorString(error)));
00225         } else if ( gpu_ro_is_current()  && (EMDATA_CPU_NEEDS_UPDATE & flags)) {
00226                 cout << "Copy ro to cpu" << endl;
00227                 copy_gpu_ro_to_cpu();
00228         }
00229         flags &= ~EMDATA_CPU_NEEDS_UPDATE;
00230 #endif
00231         return rdata;
00232 
00233 }
00234 #endif
00235 
00236 
00237 #include <sys/stat.h>
00238 
00239 void EMData::write_data(string fsp,size_t loc,const Region* area,const int file_nx, const int file_ny, const int file_nz) {
00240 
00241         if (area) {
00242                 struct stat fileinfo;
00243                 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");
00244         }
00245 
00246 
00247         FILE *f = 0;
00248         f=fopen(fsp.c_str(), "rb+");
00249         if (!f) f=fopen(fsp.c_str(), "wb");
00250         if (!f) throw FileAccessException(fsp);
00251         portable_fseek(f,loc,SEEK_SET);
00252         if (!area) {
00253                 if (fwrite(get_data(),nx*ny,nz*4,f)!=(size_t)(nz*4)) throw FileAccessException(fsp);
00254         } else {
00255                 int fnx = nx;
00256                 if (file_nx != 0) fnx = file_nx;
00257                 int fny = ny;
00258                 if (file_ny != 0) fny = file_ny;
00259                 int fnz = nz;
00260                 if (file_nz != 0) fnz = file_nz;
00261 
00262                 EMUtil::process_region_io(get_data(), f, ImageIO::READ_WRITE,
00263                                                                   0, 4,fnx,fny,fnz,area);
00264         }
00265         fclose(f);
00266 }
00267 
00268 void EMData::read_data(string fsp,size_t loc,const Region* area, const int file_nx, const int file_ny, const int file_nz) {
00269         FILE *f = 0;
00270         f=fopen(fsp.c_str(), "rb");
00271         if (!f) throw FileAccessException(fsp);
00272         int fnx = nx;
00273         if (file_nx != 0) fnx = file_nx;
00274         int fny = ny;
00275         if (file_ny != 0) fny = file_ny;
00276         int fnz = nz;
00277         if (file_nz != 0) fnz = file_nz;
00278 
00279         portable_fseek(f,loc,SEEK_SET);
00280         EMUtil::process_region_io(get_data(), f, ImageIO::READ_ONLY,
00281                                                           0, 4,fnx,fny,fnz,area);
00282 //      portable_fseek(f,loc,SEEK_SET);
00283 //      if (fread(get_data(),nx*ny,nz*4,f)!=(size_t)(nz*4)) throw FileAccessException(fsp);
00284         fclose(f);
00285 }
00286 
00287 float EMData::calc_center_density()
00288 {
00289         ENTERFUNC;
00290 
00291         float center = get_attr("mean");
00292         float sigma = get_attr("sigma");
00293         float ds = sigma / 2;
00294         size_t size = nx * ny * nz;
00295         float *d = get_data();
00296         float sigma1 = sigma / 20;
00297         float sigma2 = sigma / 1000;
00298 
00299         while (ds > sigma1) {
00300                 double sum = 0;
00301                 int norm = 0;
00302 
00303                 for (size_t i = 0; i < size; i++) {
00304                         if (fabs(d[i] - center) < ds) {
00305                                 sum += d[i];
00306                                 norm++;
00307                         }
00308                 }
00309                 if (!norm) {
00310                         break;
00311                 }
00312                 float mean = (float)(sum / norm);
00313                 if (fabs(mean - center) < sigma2) {
00314                         ds *= 0.5f;
00315                 }
00316                 center = mean;
00317         }
00318         EXITFUNC;
00319 
00320         return center;
00321 }
00322 
00323 
00324 float EMData::calc_sigma_diff()
00325 {
00326         ENTERFUNC;
00327 
00328         float *d = get_data();
00329         float mean = get_attr("mean");
00330         float sigma = get_attr("sigma");
00331 
00332         double sum_up = 0;
00333         double sum_down = 0;
00334         int nup = 0;
00335         int ndown = 0;
00336 
00337         size_t size = nx * ny * nz;
00338 
00339         for (size_t i = 0; i < size; i++) {
00340                 if (d[i] > mean) {
00341                         sum_up += Util::square(d[i] - mean);
00342                         nup++;
00343                 }
00344                 else {
00345                         sum_down += Util::square(mean - d[i]);
00346                         ndown++;
00347                 }
00348         }
00349 
00350         float sigup = std::sqrt((float)sum_up / nup);
00351         float sigdown = std::sqrt((float)sum_down / ndown);
00352         float sig_diff = fabs(sigup - sigdown) / sigma;
00353 
00354 
00355         EXITFUNC;
00356         return sig_diff;
00357 
00358 }
00359 
00360 
00361 IntPoint EMData::calc_min_location() const
00362 {
00363         ENTERFUNC;
00364 
00365         int di = 1;
00366         if (is_complex() && !is_ri()) {
00367                 di = 2;
00368         }
00369 
00370         float min = FLT_MAX;
00371         int min_x = 0;
00372         int min_y = 0;
00373         int min_z = 0;
00374         int nxy = nx * ny;
00375         float * data = get_data();
00376 
00377         for (int j = 0; j < nz; ++j) {
00378                 size_t cur_z = j * nxy;
00379 
00380                 for (int k = 0; k < ny; ++k) {
00381                         size_t cur_y = k * nx + cur_z;
00382 
00383                         for (int l = 0; l < nx; l += di) {
00384                                 float t = data[l + cur_y];
00385                                 if (t < min) {
00386                                         min_x = l;
00387                                         min_y = k;
00388                                         min_z = j;
00389                                         min = t;
00390                                 }
00391                         }
00392                 }
00393         }
00394 
00395         return IntPoint(min_x, min_y, min_z);
00396 }
00397 
00398 
00399 IntPoint EMData::calc_max_location() const
00400 {
00401         ENTERFUNC;
00402 
00403         int di = 1;
00404         if (is_complex() && !is_ri()) {
00405                 di = 2;
00406         }
00407 
00408         float max = -FLT_MAX;
00409         int max_x = 0;
00410         int max_y = 0;
00411         int max_z = 0;
00412         int nxy = nx * ny;
00413         float * data = get_data();
00414 
00415         for (int j = 0; j < nz; ++j) {
00416                 size_t cur_z = j * nxy;
00417 
00418                 for (int k = 0; k < ny; ++k) {
00419                         size_t cur_y = k * nx + cur_z;
00420 
00421                         for (int l = 0; l < nx; l += di) {
00422                                 float t = data[l + cur_y];
00423                                 if (t > max) {
00424                                         max_x = l;
00425                                         max_y = k;
00426                                         max_z = j;
00427                                         max = t;
00428                                 }
00429                         }
00430                 }
00431         }
00432 
00433         EXITFUNC;
00434         return IntPoint(max_x, max_y, max_z);
00435 }
00436 
00437 
00438 IntPoint EMData::calc_max_location_wrap(const int maxdx, const int maxdy, const int maxdz)
00439 {
00440         int maxshiftx = maxdx, maxshifty = maxdy, maxshiftz = maxdz;
00441         if (maxdx == -1) maxshiftx = get_xsize()/4;
00442         if (maxdy == -1) maxshifty = get_ysize()/4;
00443         if (maxdz == -1) maxshiftz = get_zsize()/4;
00444 
00445         float max_value = -FLT_MAX;
00446 
00447         IntPoint peak(0,0,0);
00448         for (int k = -maxshiftz; k <= maxshiftz; k++) {
00449                 for (int j = -maxshifty; j <= maxshifty; j++) {
00450                         for (int i = -maxshiftx; i <= maxshiftx; i++) {
00451 
00452                                 float value = get_value_at_wrap(i,j,k);
00453 
00454                                 if (value > max_value) {
00455                                         max_value = value;
00456                                         peak[0] = i;
00457                                         peak[1] = j;
00458                                         peak[2] = k;
00459                                 }
00460                         }
00461                 }
00462         }
00463 
00464         return peak;
00465 }
00466 
00467 FloatPoint EMData::calc_center_of_mass()
00468 {
00469         float *data = get_data();
00470 
00471         float sigma = get_attr("sigma");
00472         float mean = get_attr("mean");
00473         float m = 0;
00474 
00475         FloatPoint com(0,0,0);
00476         for (int i = 0; i < nx; ++i) {
00477                 for (int j = 0; j < ny; ++j) {
00478                         int j2 = nx * j;
00479                         for (int k = 0; k < nz; ++k) {
00480                                 size_t l = i + j2 + k * nxy;
00481                                 if (data[l] >= sigma * .75 + mean) {
00482                                         com[0] += i * data[l];
00483                                         com[1] += j * data[l];
00484                                         com[2] += k * data[l];
00485                                         m += data[l];
00486                                 }
00487                         }
00488                 }
00489         }
00490 
00491         com[0] /= m;
00492         com[1] /= m;
00493         com[2] /= m;
00494 
00495         return com;
00496 }
00497 
00498 
00499 int EMData::calc_min_index() const
00500 {
00501         IntPoint min_location = calc_min_location();
00502         int i = min_location[0] + min_location[1] * nx + min_location[2] * nx * ny;
00503         return i;
00504 }
00505 
00506 
00507 int EMData::calc_max_index() const
00508 {
00509         IntPoint max_location = calc_max_location();
00510         int i = max_location[0] + max_location[1] * nx + max_location[2] * nx * ny;
00511         return i;
00512 }
00513 
00514 
00515 vector<Pixel> EMData::calc_highest_locations(float threshold)
00516 {
00517         ENTERFUNC;
00518 
00519         vector<Pixel> result;
00520 
00521         int di = 1;
00522         if (is_complex() && !is_ri()) {
00523                 di = 2;
00524         }
00525 
00526         int nxy = nx * ny;
00527         float * data = get_data();
00528 
00529         for (int j = 0; j < nz; ++j) {
00530                 size_t cur_z = j * nxy;
00531 
00532                 for (int k = 0; k < ny; ++k) {
00533                         size_t cur_y = k * nx + cur_z;
00534 
00535                         for (int l = 0; l < nx; l += di) {
00536                                 float v =data[l + cur_y];
00537                                 if (v > threshold) {
00538                                         result.push_back(Pixel(l, k, j, v));
00539                                 }
00540                         }
00541                 }
00542         }
00543 
00544         std::sort(result.begin(), result.end());
00545 
00546         EXITFUNC;
00547         return result;
00548 }
00549 
00550 
00551 float EMData::get_edge_mean() const
00552 {
00553         ENTERFUNC;
00554 
00555         int di = 0;
00556         double edge_sum = 0;
00557         float edge_mean = 0;
00558         size_t nxy = nx * ny;
00559         float * data = get_data();
00560         if (nz == 1) {
00561                 for (int i = 0, j = (ny - 1) * nx; i < nx; ++i, ++j) {
00562                         edge_sum += data[i] + data[j];
00563                 }
00564                 for (size_t i = 0, j = nx - 1; i < nxy; i += nx, j += nx) {
00565                         edge_sum += data[i] + data[j];
00566                 }
00567                 edge_mean = (float)edge_sum / (nx * 2 + ny * 2);
00568         }
00569         else {
00570                 if (nx == ny && nx == nz * 2 - 1) {
00571                         for (size_t j = (nxy * (nz - 1)); j < nxy * nz; ++j, ++di) {
00572                                 edge_sum += data[j];
00573                         }
00574                 }
00575                 else {
00576                         for (size_t i = 0, j = (nxy * (nz - 1)); i < nxy; ++i, ++j, ++di) {
00577                                 edge_sum += data[i] + data[j];
00578                         }
00579                 }
00580 
00581                 int nxy2 = nx * (ny - 1);
00582                 for (int k = 1; k < nz - 1; ++k) {
00583                         size_t k2 = k * nxy;
00584                         size_t k3 = k2 + nxy2;
00585                         for (int i = 0; i < nx; ++i, ++di) {
00586                                 edge_sum += data[i + k2] + data[i + k3];
00587                         }
00588                 }
00589                 for (int k = 1; k < nz - 1; ++k) {
00590                         size_t k2 = k * nxy;
00591                         size_t k3 = nx - 1 + k2;
00592                         for (int i = 1; i < ny - 1; ++i, ++di) {
00593                                 edge_sum += data[i * nx + k2] + data[i * nx + k3];
00594                         }
00595                 }
00596 
00597                 edge_mean = (float)edge_sum / (di * 2);
00598         }
00599         EXITFUNC;
00600 
00601         return edge_mean;
00602 }
00603 
00604 
00605 float EMData::get_circle_mean()
00606 {
00607         ENTERFUNC;
00608 
00609         static bool busy = false;
00610         static EMData *mask = 0;
00611 
00612         while (busy);
00613         busy = true;
00614 
00615         if (!mask || !EMUtil::is_same_size(this, mask)) {
00616                 if (!mask) {
00617                         mask = new EMData();
00618                 }
00619                 mask->set_size(nx, ny, nz);
00620                 mask->to_one();
00621 
00622                 float radius = (float)(ny / 2 - 2);
00623                 mask->process_inplace("mask.sharp", Dict("inner_radius", radius - 1,
00624                                                                            "outer_radius", radius + 1));
00625 
00626         }
00627         double n = 0,s=0;
00628         float *d = mask->get_data();
00629         float * data = get_data();
00630         size_t size = nx*ny*nz;
00631         for (size_t i = 0; i < size; i++) {
00632                 if (d[i]) { n+=1.0; s+=data[i]; }
00633         }
00634 
00635 
00636         float result = (float)(s/n);
00637         busy = false;
00638 
00639         EXITFUNC;
00640         return result;
00641 }
00642 
00643 
00644 void EMData::set_ctf(Ctf * new_ctf)
00645 {
00646         ENTERFUNC;
00647 
00648         vector<float> vctf = new_ctf->to_vector();
00649         attr_dict["ctf"] = vctf;
00650 
00651         EXITFUNC;
00652 }
00653 
00654 Ctf * EMData::get_ctf() const
00655 {
00656         if(attr_dict.has_key("ctf")) {
00657                 EMAN1Ctf * ctf = new EMAN1Ctf();
00658                 ctf->from_vector(attr_dict["ctf"]);
00659 
00660                 return dynamic_cast<Ctf *>(ctf);
00661         }
00662         else {
00663                 return 0;
00664         }
00665 }
00666 
00667 #include <iostream>
00668 using std::cout;
00669 using std::endl;
00670 
00671 void EMData::set_size(int x, int y, int z)
00672 {
00673         ENTERFUNC;
00674 
00675         if (x <= 0) {
00676                 throw InvalidValueException(x, "x size <= 0");
00677         }
00678         else if (y <= 0) {
00679                 throw InvalidValueException(y, "y size <= 0");
00680         }
00681         else if (z <= 0) {
00682                 throw InvalidValueException(z, "z size <= 0");
00683         }
00684 
00685         int old_nx = nx;
00686 
00687         size_t size = (size_t)(x) * (size_t)y * (size_t)z * sizeof(float);
00688 
00689         if (rdata != 0) {
00690                 rdata = (float*)EMUtil::em_realloc(rdata,size);
00691         } else {
00692                 // Just pass on this for a while....see what happens
00693                 rdata = (float*)EMUtil::em_malloc(size);
00694         }
00695 //      rdata = static_cast < float *>(realloc(rdata, size));
00696         if ( rdata == 0 )
00697         {
00698                 stringstream ss;
00699                 string gigs;
00700                 ss << (float) size/1000000000.0;
00701                 ss >> gigs;
00702                 string message = "Cannot allocate " + gigs + " GB - not enough memory.";
00703                 throw BadAllocException(message);
00704         }
00705 
00706 #ifdef EMAN2_USING_CUDA
00707         // This is important
00708         free_cuda_memory();
00709 #endif // EMAN2_USING_CUDA
00710 
00711         nx = x;
00712         ny = y;
00713         nz = z;
00714         nxy = nx*ny;
00715         nxyz = nx*ny*nz;
00716 
00717         if (old_nx == 0) {
00718                 EMUtil::em_memset(get_data(),0,size);
00719         }
00720 
00721         if (supp) {
00722                 EMUtil::em_free(supp);
00723                 supp = 0;
00724         }
00725 
00726         update();
00727         EXITFUNC;
00728 }
00729 
00730 #ifdef EMAN2_USING_CUDA
00731 
00732 #include <cuda_runtime_api.h>
00733 #include "cuda/cuda_util.h"
00734 
00735 void EMData::set_size_cuda(int x, int y, int z)
00736 {
00737         ENTERFUNC;
00738 
00739         if (x <= 0) {
00740                 throw InvalidValueException(x, "x size <= 0");
00741         }
00742         else if (y <= 0) {
00743                 throw InvalidValueException(y, "y size <= 0");
00744         }
00745         else if (z <= 0) {
00746                 throw InvalidValueException(z, "z size <= 0");
00747         }
00748 
00749         if (cuda_cache_handle!=-1) {
00750                 cuda_cache.clear_item(cuda_cache_handle);
00751                 cuda_cache_handle = -1;
00752         }
00753         nx = x;
00754         ny = y;
00755         nz = z;
00756 
00757         nxy = nx*ny;
00758 
00759         get_cuda_data();
00760 
00761 //      cuda_cache_handle = cuda_rw_cache.cache_data(this,rdata,nx,ny,nz); Let's be lazy
00762         // This is important
00763         free_memory(); // Now release CPU memory, seeing as a GPU resize invalidates it - Actually let's be lazy about it instead
00764 
00765         gpu_update();
00766 
00767         EXITFUNC;
00768 }
00769 
00770 #endif // EMAN2_USING_CUDA
00771 
00772 MArray2D EMData::get_2dview() const
00773 {
00774         const int ndims = 2;
00775         if (get_ndim() != ndims) {
00776                 throw ImageDimensionException("2D only");
00777         }
00778         boost::array<std::size_t,ndims> dims = {{nx, ny}};
00779         MArray2D marray(get_data(), dims, boost::fortran_storage_order());
00780         return marray;
00781 }
00782 
00783 
00784 MArray3D EMData::get_3dview() const
00785 {
00786         const int ndims = 3;
00787         boost::array<std::size_t,ndims> dims = {{nx, ny, nz}};
00788         MArray3D marray(get_data(), dims, boost::fortran_storage_order());
00789         return marray;
00790 }
00791 
00792 
00793 MCArray2D EMData::get_2dcview() const
00794 {
00795         const int ndims = 2;
00796         if (get_ndim() != ndims) {
00797                 throw ImageDimensionException("2D only");
00798         }
00799         boost::array<std::size_t,ndims> dims = {{nx/2, ny}};
00800         std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
00801         MCArray2D marray(cdata, dims, boost::fortran_storage_order());
00802         return marray;
00803 }
00804 
00805 
00806 MCArray3D EMData::get_3dcview() const
00807 {
00808         const int ndims = 3;
00809         boost::array<std::size_t,ndims> dims = {{nx/2, ny, nz}};
00810         std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
00811         MCArray3D marray(cdata, dims, boost::fortran_storage_order());
00812         return marray;
00813 }
00814 
00815 
00816 MCArray3D* EMData::get_3dcviewptr() const
00817 {
00818         const int ndims = 3;
00819         boost::array<std::size_t,ndims> dims = {{nx/2, ny, nz}};
00820         std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
00821         MCArray3D* marray = new MCArray3D(cdata, dims,
00822                                                                           boost::fortran_storage_order());
00823         return marray;
00824 }
00825 
00826 
00827 MArray2D EMData::get_2dview(int x0, int y0) const
00828 {
00829         const int ndims = 2;
00830         if (get_ndim() != ndims) {
00831                 throw ImageDimensionException("2D only");
00832         }
00833         boost::array<std::size_t,ndims> dims = {{nx, ny}};
00834         MArray2D marray(get_data(), dims, boost::fortran_storage_order());
00835         boost::array<std::size_t,ndims> bases={{x0, y0}};
00836         marray.reindex(bases);
00837         return marray;
00838 }
00839 
00840 
00841 MArray3D EMData::get_3dview(int x0, int y0, int z0) const
00842 {
00843         const int ndims = 3;
00844         boost::array<std::size_t,ndims> dims = {{nx, ny, nz}};
00845         MArray3D marray(get_data(), dims, boost::fortran_storage_order());
00846         boost::array<std::size_t,ndims> bases={{x0, y0, z0}};
00847         marray.reindex(bases);
00848         return marray;
00849 }
00850 
00851 
00852 MCArray2D EMData::get_2dcview(int x0, int y0) const
00853 {
00854         const int ndims = 2;
00855         if (get_ndim() != ndims) {
00856                 throw ImageDimensionException("2D only");
00857         }
00858         boost::array<std::size_t,ndims> dims = {{nx/2, ny}};
00859         std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
00860         MCArray2D marray(cdata, dims, boost::fortran_storage_order());
00861         boost::array<std::size_t,ndims> bases={{x0, y0}};
00862         marray.reindex(bases);
00863         return marray;
00864 }
00865 
00866 
00867 MCArray3D EMData::get_3dcview(int x0, int y0, int z0) const
00868 {
00869         const int ndims = 3;
00870         boost::array<std::size_t,ndims> dims = {{nx/2, ny, nz}};
00871         std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
00872         MCArray3D marray(cdata, dims, boost::fortran_storage_order());
00873         boost::array<std::size_t,ndims> bases={{x0, y0, z0}};
00874         marray.reindex(bases);
00875         return marray;
00876 }
00877 
00878 int greaterthan( const void* p1, const void* p2 )
00879 {
00880         float*  v1 = (float*) p1;
00881         float*  v2 = (float*) p2;
00882 
00883         if ( *v1 < *v2 )  return 0;
00884         else return 1;
00885 }
00886 
00887 
00888 EMObject EMData::get_attr(const string & key) const
00889 {
00890         ENTERFUNC;
00891 
00892         size_t size = nx * ny * nz;
00893         update_stat();
00894 
00895         if (key == "kurtosis") {
00896                 float mean = attr_dict["mean"];
00897                 float sigma = attr_dict["sigma"];
00898 
00899                 float *data = get_data();
00900                 double kurtosis_sum = 0;
00901 
00902                 for (size_t k = 0; k < size; k++) {
00903                         float t = (data[k] - mean) / sigma;
00904                         float tt = t * t;
00905                         kurtosis_sum += tt * tt;
00906                 }
00907 
00908                 float kurtosis = (float)(kurtosis_sum / size - 3.0);
00909                 return kurtosis;
00910         }
00911         else if (key == "skewness") {
00912                 float mean = attr_dict["mean"];
00913                 float sigma = attr_dict["sigma"];
00914 
00915                 float *data = get_data();
00916                 double skewness_sum = 0;
00917                 for (size_t k = 0; k < size; k++) {
00918                         float t = (data[k] - mean) / sigma;
00919                         skewness_sum +=  t * t * t;
00920                 }
00921                 float skewness = (float)(skewness_sum / size);
00922                 return skewness;
00923         }
00924         else if (key == "median")
00925         {
00926                 if ( is_complex() ) throw ImageFormatException("Error - can not calculate the median of a complex image");
00927                 size_t n = size;
00928                 float* tmp = new float[n];
00929                 float* d = get_data();
00930                 if (tmp == 0 ) throw BadAllocException("Error - could not create deep copy of image data");
00931                 for(size_t i=0; i < n; ++i) tmp[i] = d[i]; // should just be a memcpy
00932                 qsort(tmp, n, sizeof(float), &greaterthan);
00933                 float median;
00934                 if (n%2==1) median = tmp[n/2];
00935                 else median = (tmp[n/2-1]+tmp[n/2])/2.0f;
00936                 delete [] tmp;
00937                 return median;
00938         }
00939         else if (key == "nonzero_median")
00940         {
00941                 if ( is_complex() ) throw ImageFormatException("Error - can not calculate the median of a complex image");
00942                 vector<float> tmp;
00943                 size_t n = size;
00944                 float* d = get_data();
00945                 for( size_t i = 0; i < n; ++i ) {
00946                         if ( d[i] != 0 ) tmp.push_back(d[i]);
00947                 }
00948                 sort(tmp.begin(), tmp.end());
00949                 unsigned int vsize = tmp.size();
00950                 float median;
00951                 if (vsize%2==1) median = tmp[vsize/2];
00952                 else median = (tmp[vsize/2-1]+tmp[vsize/2])/2.0f;
00953                 return median;
00954         }
00955         else if (key == "changecount") return EMObject(changecount);
00956         else if (key == "nx") return nx;
00957         else if (key == "ny") return ny;
00958         else if (key == "nz") return nz;
00959 
00960         if(attr_dict.has_key(key)) {
00961                 return attr_dict[key];
00962         }
00963         else {
00964                 throw NotExistingObjectException(key, "The requested key does not exist");
00965         }
00966 
00967         EXITFUNC;
00968 }
00969 
00970 EMObject EMData::get_attr_default(const string & key, const EMObject & em_obj) const
00971 {
00972         ENTERFUNC;
00973 
00974         if(attr_dict.has_key(key)) {
00975                 return get_attr(key);
00976         }
00977         else {
00978                 return em_obj;
00979         }
00980 
00981         EXITFUNC;
00982 }
00983 
00984 Dict EMData::get_attr_dict() const
00985 {
00986         update_stat();
00987         
00988         Dict tmp=Dict(attr_dict);
00989         tmp["nx"]=nx;
00990         tmp["ny"]=ny;
00991         tmp["nz"]=nz;
00992         tmp["changecount"]=changecount;
00993 
00994         return tmp;
00995 }
00996 
00997 void EMData::set_attr_dict(const Dict & new_dict)
00998 {
00999         /*set nx, ny nz may resize the image*/
01000         // This wasn't supposed to 'clip' the image, but just redefine the size --steve
01001         if( ( new_dict.has_key("nx") && nx!=(int)new_dict["nx"] )
01002                 || ( new_dict.has_key("ny") && ny!=(int)new_dict["ny"] )
01003                 || ( new_dict.has_key("nz") && nz!=(int)new_dict["nz"] ) ) {
01004 
01005                 int newx, newy, newz;
01006                 newx = new_dict.has_key("nx") ? (int)new_dict["nx"] : nx;
01007                 newy = new_dict.has_key("ny") ? (int)new_dict["ny"] : ny;
01008                 newz = new_dict.has_key("nz") ? (int)new_dict["nz"] : nz;
01009 
01010                 set_size(newx,newy,newz);
01011 
01012 //              EMData * new_image = get_clip(Region((nx-newx)/2, (ny-newy)/2, (nz=newz)/2, newx, newy, newz));
01013 //              if(new_image) {
01014 //                      this->operator=(*new_image);
01015 //                      delete new_image;
01016 //                      new_image = 0;
01017 //              }
01018         }
01019 
01020         vector<string> new_keys = new_dict.keys();
01021         vector<string>::const_iterator it;
01022         for(it = new_keys.begin(); it!=new_keys.end(); ++it) {
01023                 this->set_attr(*it, new_dict[*it]);
01024         }
01025 }
01026 
01027 void EMData::set_attr_dict_explicit(const Dict & new_dict)
01028 {
01029         attr_dict = new_dict;
01030 }
01031 
01032 void EMData::del_attr(const string & attr_name)
01033 {
01034         attr_dict.erase(attr_name);
01035 }
01036 
01037 void EMData::del_attr_dict(const vector<string> & del_keys)
01038 {
01039         vector<string>::const_iterator it;
01040         for(it=del_keys.begin(); it!=del_keys.end(); ++it) {
01041                 this->del_attr(*it);
01042         }
01043 }
01044 
01045 void EMData::set_attr(const string & key, EMObject val)
01046 {
01047         if( key == "nx" && nx != (int)val) { set_size((int)val,ny,nz); return; }
01048         if( key == "ny" && ny != (int)val) { set_size(nx,(int)val,nz); return; }
01049         if( key == "nz" && nz != (int)val) { set_size(nx,ny,(int)val); return; }
01050 
01051         /* Ignore 'read only' attribute. */
01052         if(key == "sigma" ||
01053                 key == "sigma_nonzero" ||
01054                 key == "square_sum" ||
01055                 key == "maximum" ||
01056                 key == "minimum" ||
01057                 key == "mean" ||
01058                 key == "mean_nonzero" )
01059         {
01060                 LOGWARN("Ignore setting read only attribute %s", key.c_str());
01061                 return;
01062         }
01063 
01064         attr_dict[key] = val;
01065 
01066 
01067 
01068 }
01069 
01070 void EMData::set_attr_python(const string & key, EMObject val)
01071 {
01072         if( key == "nx" && nx != (int)val) { set_size((int)val,ny,nz); return; }
01073         if( key == "ny" && ny != (int)val) { set_size(nx,(int)val,nz); return; }
01074         if( key == "nz" && nz != (int)val) { set_size(nx,ny,(int)val); return; }
01075 
01076         /* Ignore 'read only' attribute. */
01077         if(key == "sigma" ||
01078                   key == "sigma_nonzero" ||
01079                   key == "square_sum" ||
01080                   key == "maximum" ||
01081                   key == "minimum" ||
01082                   key == "mean" ||
01083                   key == "mean_nonzero" )
01084         {
01085                 LOGWARN("Ignore setting read only attribute %s", key.c_str());
01086                 return;
01087         }
01088 
01089         EMObject::ObjectType argtype = val.get_type();
01090         if (argtype == EMObject::EMDATA) {
01091                 EMData* e = (EMData*) val;
01092                 e = e->copy();
01093                 EMObject v(e);
01094                 attr_dict[key] = v;
01095         }
01096         else if (argtype == EMObject::TRANSFORM) {
01097                 Transform* t = new Transform(*((Transform*) val));
01098                 EMObject v(t);
01099                 attr_dict[key] = v;
01100                 delete t; t=0;
01101         } else {
01102                 attr_dict[key] = val;
01103         }
01104 
01105 }
01106 
01107 void EMData::scale_pixel(float scale) const
01108 {
01109         attr_dict["apix_x"] = ((float) attr_dict["apix_x"]) * scale;
01110         attr_dict["apix_y"] = ((float) attr_dict["apix_y"]) * scale;
01111         attr_dict["apix_z"] = ((float) attr_dict["apix_z"]) * scale;
01112         if (attr_dict.has_key("ctf")) {
01113                 Ctf *ctf=(Ctf *)attr_dict["ctf"];
01114                 ctf->apix*=scale;
01115                 attr_dict["ctf"]=ctf;
01116                 if(ctf) {delete ctf; ctf=0;}
01117         }
01118 }
01119 
01120 
01121 //vector<float> EMData::get_data_pickle() const
01122 std::string EMData::get_data_pickle() const
01123 {
01124 //      vector<float> vf;
01125 //      vf.resize(nx*ny*nz);
01126 //      std::copy(rdata, rdata+nx*ny*nz, vf.begin());
01127 
01128         std::string vf((const char *)get_data(),nx*ny*nz*sizeof(float));
01129 
01130         return vf;
01131 }
01132 
01133 //void EMData::set_data_pickle(const vector<float>& vf)
01134 void EMData::set_data_pickle(std::string vf)
01135 {
01136 //      if (rdata) printf("rdata exists\n");
01137 //      rdata = (float *)malloc(nx*ny*nz*sizeof(float));
01138 //      std::copy(vf.begin(), vf.end(), rdata);
01139         EMUtil::em_memcpy(get_data(),vf.data(),nx*ny*nz*sizeof(float));
01140 
01141 }
01142 
01143 int EMData::get_supp_pickle() const
01144 {
01145         return 0;
01146 }
01147 
01148 void EMData::set_supp_pickle(int)
01149 {
01150         this->supp = 0;
01151 }
01152 
01153 
01154 vector<Vec3i > find_region(EMData* image,const vector<Vec3i >& coords, const float value, vector<Vec3i >& region)
01155 {
01156         static vector<Vec3i> two_six_connected;
01157         if (two_six_connected.size() == 0) {
01158                 for(int i = -1; i <= 1; ++i) {
01159                         for(int j = -1; j <= 1; ++j) {
01160                                 for(int  k = -1; k <= 1; ++k) {
01161                                         if ( j != 0 || i != 0 || k != 0) {
01162                                                 two_six_connected.push_back(Vec3i(i,j,k));
01163                                         }
01164                                 }
01165                         }
01166                 }
01167         }
01168 
01169         vector<Vec3i> ret;
01170         for(vector<Vec3i>::const_iterator it = two_six_connected.begin(); it != two_six_connected.end(); ++it ) {
01171                 for(vector<Vec3i>::const_iterator it2 = coords.begin(); it2 != coords.end(); ++it2 ) {
01172                         if  (image->get_value_at((*it2)[0],(*it2)[1],(*it2)[2]) != value) throw;
01173                         Vec3i c = (*it)+(*it2);
01174 
01175                         if ( c[0] < 0 || c[0] >= image->get_xsize()) continue;
01176                         if ( c[1] < 0 || c[1] >= image->get_ysize()) continue;
01177                         if ( c[2] < 0 || c[2] >= image->get_zsize()) continue;
01178 
01179                         if( image->get_value_at(c[0],c[1],c[2]) == value ) {
01180                                 if (find(ret.begin(),ret.end(),c) == ret.end()) {
01181                                         if (find(region.begin(),region.end(),c) == region.end()) {
01182                                                 region.push_back(c);
01183                                                 ret.push_back(c);
01184                                         }
01185                                 }
01186                         }
01187                 }
01188         }
01189         return ret;
01190 }
01191 
01192 vector<Vec3i> EMData::mask_contig_region(const float& value, const Vec3i& seed) {
01193         Vec3i coord(seed[0],seed[1],seed[2]);
01194         vector<Vec3i> region;
01195         region.push_back(coord);
01196         vector<Vec3i> find_region_input = region;
01197         while (true) {
01198                 vector<Vec3i> v = find_region(this,find_region_input, value, region);
01199                 if (v.size() == 0 ) break;
01200                 else find_region_input = v;
01201         }
01202         return region;
01203 }
01204 

Generated on Fri Mar 19 02:19:21 2010 for EMAN2 by  doxygen 1.5.6