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

Generated on Sat Nov 21 02:19:15 2009 for EMAN2 by  doxygen 1.5.6