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