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