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