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)
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 (rdata != 0) {
00857 rdata = (float*)EMUtil::em_realloc(rdata,size);
00858 } else {
00859
00860 rdata = (float*)EMUtil::em_malloc(size);
00861 }
00862
00863 if ( rdata == 0 )
00864 {
00865 stringstream ss;
00866 string gigs;
00867 ss << (float) size/1000000000.0;
00868 ss >> gigs;
00869 string message = "Cannot allocate " + gigs + " GB - not enough memory.";
00870 throw BadAllocException(message);
00871 }
00872
00873 nx = x;
00874 ny = y;
00875 nz = z;
00876 nxy = nx*ny;
00877 nxyz = (size_t)nx*ny*nz;
00878
00879
00880 #ifdef EMAN2_USING_CUDA
00881 if(cudarwdata) {
00882
00883 rw_free();
00884 rw_alloc();
00885 }
00886 if(cudarodata) {
00887 ro_free();
00888 ro_alloc();
00889 }
00890 #endif // EMAN2_USING_CUDA
00891
00892 if (old_nx == 0) {
00893 EMUtil::em_memset(get_data(),0,size);
00894 }
00895
00896 if (supp) {
00897 EMUtil::em_free(supp);
00898 supp = 0;
00899 }
00900
00901 update();
00902 EXITFUNC;
00903 }
00904
00905 #ifdef EMAN2_USING_CUDA
00906
00907 void EMData::set_size_cuda(int x, int y, int z)
00908 {
00909 ENTERFUNC;
00910
00911 if (x <= 0) {
00912 throw InvalidValueException(x, "x size <= 0");
00913 }
00914 else if (y <= 0) {
00915 throw InvalidValueException(y, "y size <= 0");
00916 }
00917 else if (z <= 0) {
00918 throw InvalidValueException(z, "z size <= 0");
00919 }
00920
00921 nx = x;
00922 ny = y;
00923 nz = z;
00924
00925 nxy = nx*ny;
00926 nxyz = (size_t)nx*ny*nz;
00927
00928
00929
00930 EXITFUNC;
00931 }
00932
00933 #endif // EMAN2_USING_CUDA
00934
00935 MArray2D EMData::get_2dview() const
00936 {
00937 const int ndims = 2;
00938 if (get_ndim() != ndims) {
00939 throw ImageDimensionException("2D only");
00940 }
00941 boost::array<std::size_t,ndims> dims = {{nx, ny}};
00942 MArray2D marray(get_data(), dims, boost::fortran_storage_order());
00943 return marray;
00944 }
00945
00946
00947 MArray3D EMData::get_3dview() const
00948 {
00949 const int ndims = 3;
00950 boost::array<std::size_t,ndims> dims = {{nx, ny, nz}};
00951 MArray3D marray(get_data(), dims, boost::fortran_storage_order());
00952 return marray;
00953 }
00954
00955
00956 MCArray2D EMData::get_2dcview() const
00957 {
00958 const int ndims = 2;
00959 if (get_ndim() != ndims) {
00960 throw ImageDimensionException("2D only");
00961 }
00962 boost::array<std::size_t,ndims> dims = {{nx/2, ny}};
00963 std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
00964 MCArray2D marray(cdata, dims, boost::fortran_storage_order());
00965 return marray;
00966 }
00967
00968
00969 MCArray3D EMData::get_3dcview() const
00970 {
00971 const int ndims = 3;
00972 boost::array<std::size_t,ndims> dims = {{nx/2, ny, nz}};
00973 std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
00974 MCArray3D marray(cdata, dims, boost::fortran_storage_order());
00975 return marray;
00976 }
00977
00978
00979 MCArray3D* EMData::get_3dcviewptr() const
00980 {
00981 const int ndims = 3;
00982 boost::array<std::size_t,ndims> dims = {{nx/2, ny, nz}};
00983 std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
00984 MCArray3D* marray = new MCArray3D(cdata, dims,
00985 boost::fortran_storage_order());
00986 return marray;
00987 }
00988
00989
00990 MArray2D EMData::get_2dview(int x0, int y0) const
00991 {
00992 const int ndims = 2;
00993 if (get_ndim() != ndims) {
00994 throw ImageDimensionException("2D only");
00995 }
00996 boost::array<std::size_t,ndims> dims = {{nx, ny}};
00997 MArray2D marray(get_data(), dims, boost::fortran_storage_order());
00998 boost::array<std::size_t,ndims> bases={{x0, y0}};
00999 marray.reindex(bases);
01000 return marray;
01001 }
01002
01003
01004 MArray3D EMData::get_3dview(int x0, int y0, int z0) const
01005 {
01006 const int ndims = 3;
01007 boost::array<std::size_t,ndims> dims = {{nx, ny, nz}};
01008 MArray3D marray(get_data(), dims, boost::fortran_storage_order());
01009 boost::array<std::size_t,ndims> bases={{x0, y0, z0}};
01010 marray.reindex(bases);
01011 return marray;
01012 }
01013
01014
01015 MCArray2D EMData::get_2dcview(int x0, int y0) const
01016 {
01017 const int ndims = 2;
01018 if (get_ndim() != ndims) {
01019 throw ImageDimensionException("2D only");
01020 }
01021 boost::array<std::size_t,ndims> dims = {{nx/2, ny}};
01022 std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
01023 MCArray2D marray(cdata, dims, boost::fortran_storage_order());
01024 boost::array<std::size_t,ndims> bases={{x0, y0}};
01025 marray.reindex(bases);
01026 return marray;
01027 }
01028
01029
01030 MCArray3D EMData::get_3dcview(int x0, int y0, int z0) const
01031 {
01032 const int ndims = 3;
01033 boost::array<std::size_t,ndims> dims = {{nx/2, ny, nz}};
01034 std::complex<float>* cdata = reinterpret_cast<std::complex<float>*>(get_data());
01035 MCArray3D marray(cdata, dims, boost::fortran_storage_order());
01036 boost::array<std::size_t,ndims> bases={{x0, y0, z0}};
01037 marray.reindex(bases);
01038 return marray;
01039 }
01040
01041 int greaterthan( const void* p1, const void* p2 )
01042 {
01043 float* v1 = (float*) p1;
01044 float* v2 = (float*) p2;
01045
01046 if ( *v1 < *v2 ) return 0;
01047 else return 1;
01048 }
01049
01050
01051 EMObject EMData::get_attr(const string & key) const
01052 {
01053 ENTERFUNC;
01054
01055 if ((flags & EMDATA_NEEDUPD) && (key != "is_fftpad") && (key != "xform.align2d")){update_stat();}
01056
01057 size_t size = (size_t)nx * ny * nz;
01058 if (key == "kurtosis") {
01059 float mean = attr_dict["mean"];
01060 float sigma = attr_dict["sigma"];
01061
01062 float *data = get_data();
01063 double kurtosis_sum = 0;
01064
01065 for (size_t k = 0; k < size; ++k) {
01066 float t = (data[k] - mean) / sigma;
01067 float tt = t * t;
01068 kurtosis_sum += tt * tt;
01069 }
01070
01071 float kurtosis = (float)(kurtosis_sum / size - 3.0);
01072 return kurtosis;
01073 }
01074 else if (key == "skewness") {
01075 float mean = attr_dict["mean"];
01076 float sigma = attr_dict["sigma"];
01077
01078 float *data = get_data();
01079 double skewness_sum = 0;
01080 for (size_t k = 0; k < size; ++k) {
01081 float t = (data[k] - mean) / sigma;
01082 skewness_sum += t * t * t;
01083 }
01084 float skewness = (float)(skewness_sum / size);
01085 return skewness;
01086 }
01087 else if (key == "median")
01088 {
01089 if ( is_complex() ) throw ImageFormatException("Error - can not calculate the median of a complex image");
01090 size_t n = size;
01091 float* tmp = new float[n];
01092 float* d = get_data();
01093 if (tmp == 0 ) throw BadAllocException("Error - could not create deep copy of image data");
01094
01095 std::copy(d, d+n, tmp);
01096 qsort(tmp, n, sizeof(float), &greaterthan);
01097 float median;
01098 if (n%2==1) median = tmp[n/2];
01099 else median = (tmp[n/2-1]+tmp[n/2])/2.0f;
01100 delete [] tmp;
01101 return median;
01102 }
01103 else if (key == "nonzero_median")
01104 {
01105 if ( is_complex() ) throw ImageFormatException("Error - can not calculate the median of a complex image");
01106 vector<float> tmp;
01107 size_t n = size;
01108 float* d = get_data();
01109 for( size_t i = 0; i < n; ++i ) {
01110 if ( d[i] != 0 ) tmp.push_back(d[i]);
01111 }
01112 sort(tmp.begin(), tmp.end());
01113 unsigned int vsize = tmp.size();
01114 float median;
01115 if (vsize%2==1) median = tmp[vsize/2];
01116 else median = (tmp[vsize/2-1]+tmp[vsize/2])/2.0f;
01117 return median;
01118 }
01119 else if (key == "changecount") return EMObject(changecount);
01120 else if (key == "nx") return nx;
01121 else if (key == "ny") return ny;
01122 else if (key == "nz") return nz;
01123
01124 if(attr_dict.has_key(key)) {
01125 return attr_dict[key];
01126 }
01127 else {
01128 throw NotExistingObjectException(key, "The requested key does not exist");
01129 }
01130
01131 EXITFUNC;
01132 }
01133
01134 EMObject EMData::get_attr_default(const string & key, const EMObject & em_obj) const
01135 {
01136 ENTERFUNC;
01137
01138 if(attr_dict.has_key(key)) {
01139 return get_attr(key);
01140 }
01141 else {
01142 return em_obj;
01143 }
01144
01145 EXITFUNC;
01146 }
01147
01148 Dict EMData::get_attr_dict() const
01149 {
01150 update_stat();
01151
01152 Dict tmp=Dict(attr_dict);
01153 tmp["nx"]=nx;
01154 tmp["ny"]=ny;
01155 tmp["nz"]=nz;
01156 tmp["changecount"]=changecount;
01157
01158 return tmp;
01159 }
01160
01161 void EMData::set_attr_dict(const Dict & new_dict)
01162 {
01163
01164
01165 if( ( new_dict.has_key("nx") && nx!=(int)new_dict["nx"] )
01166 || ( new_dict.has_key("ny") && ny!=(int)new_dict["ny"] )
01167 || ( new_dict.has_key("nz") && nz!=(int)new_dict["nz"] ) ) {
01168
01169 int newx, newy, newz;
01170 newx = new_dict.has_key("nx") ? (int)new_dict["nx"] : nx;
01171 newy = new_dict.has_key("ny") ? (int)new_dict["ny"] : ny;
01172 newz = new_dict.has_key("nz") ? (int)new_dict["nz"] : nz;
01173
01174 set_size(newx,newy,newz);
01175
01176
01177
01178
01179
01180
01181
01182 }
01183
01184 vector<string> new_keys = new_dict.keys();
01185 vector<string>::const_iterator it;
01186 for(it = new_keys.begin(); it!=new_keys.end(); ++it) {
01187 this->set_attr(*it, new_dict[*it]);
01188 }
01189 }
01190
01191 void EMData::set_attr_dict_explicit(const Dict & new_dict)
01192 {
01193 attr_dict = new_dict;
01194 }
01195
01196 void EMData::del_attr(const string & attr_name)
01197 {
01198 attr_dict.erase(attr_name);
01199 }
01200
01201 void EMData::del_attr_dict(const vector<string> & del_keys)
01202 {
01203 vector<string>::const_iterator it;
01204 for(it=del_keys.begin(); it!=del_keys.end(); ++it) {
01205 this->del_attr(*it);
01206 }
01207 }
01208
01209 void EMData::set_attr(const string & key, EMObject val)
01210 {
01211 if( key == "nx" && nx != (int)val) { set_size((int)val,ny,nz); return; }
01212 if( key == "ny" && ny != (int)val) { set_size(nx,(int)val,nz); return; }
01213 if( key == "nz" && nz != (int)val) { set_size(nx,ny,(int)val); return; }
01214
01215
01216 if(key == "sigma" ||
01217 key == "sigma_nonzero" ||
01218 key == "square_sum" ||
01219 key == "maximum" ||
01220 key == "minimum" ||
01221 key == "mean" ||
01222 key == "mean_nonzero" )
01223 {
01224 LOGWARN("Ignore setting read only attribute %s", key.c_str());
01225 return;
01226 }
01227
01228 attr_dict[key] = val;
01229
01230
01231
01232 }
01233
01234 void EMData::set_attr_python(const string & key, EMObject val)
01235 {
01236 if( key == "nx" && nx != (int)val) { set_size((int)val,ny,nz); return; }
01237 if( key == "ny" && ny != (int)val) { set_size(nx,(int)val,nz); return; }
01238 if( key == "nz" && nz != (int)val) { set_size(nx,ny,(int)val); return; }
01239
01240
01241 if(key == "sigma" ||
01242 key == "sigma_nonzero" ||
01243 key == "square_sum" ||
01244 key == "maximum" ||
01245 key == "minimum" ||
01246 key == "mean" ||
01247 key == "mean_nonzero" )
01248 {
01249 LOGWARN("Ignore setting read only attribute %s", key.c_str());
01250 return;
01251 }
01252
01253 EMObject::ObjectType argtype = val.get_type();
01254 if (argtype == EMObject::EMDATA) {
01255 EMData* e = (EMData*) val;
01256 e = e->copy();
01257 EMObject v(e);
01258 attr_dict[key] = v;
01259 }
01260 else if (argtype == EMObject::TRANSFORM) {
01261 Transform* t = new Transform(*((Transform*) val));
01262 EMObject v(t);
01263 attr_dict[key] = v;
01264 delete t; t=0;
01265 } else {
01266 attr_dict[key] = val;
01267 }
01268
01269 }
01270
01271 void EMData::scale_pixel(float scale) const
01272 {
01273 attr_dict["apix_x"] = ((float) attr_dict["apix_x"]) * scale;
01274 attr_dict["apix_y"] = ((float) attr_dict["apix_y"]) * scale;
01275 attr_dict["apix_z"] = ((float) attr_dict["apix_z"]) * scale;
01276 if (attr_dict.has_key("ctf")) {
01277 Ctf *ctf=(Ctf *)attr_dict["ctf"];
01278 ctf->apix*=scale;
01279 attr_dict["ctf"]=ctf;
01280 if(ctf) {delete ctf; ctf=0;}
01281 }
01282 }
01283
01284
01285
01286 std::string EMData::get_data_pickle() const
01287 {
01288
01289
01290
01291
01292 std::string vf((const char *)get_data(),nx*ny*nz*sizeof(float));
01293
01294 return vf;
01295 }
01296
01297
01298 void EMData::set_data_pickle(std::string vf)
01299 {
01300
01301
01302
01303 EMUtil::em_memcpy(get_data(),vf.data(),nx*ny*nz*sizeof(float));
01304
01305 }
01306
01307 int EMData::get_supp_pickle() const
01308 {
01309 return 0;
01310 }
01311
01312 void EMData::set_supp_pickle(int)
01313 {
01314 this->supp = 0;
01315 }
01316
01317 float EMData::get_amplitude_thres(float thres)
01318 {
01319
01320 if (thres < 0 || thres > 1){
01321 LOGERR("threshold bust be between 0 and 1.");
01322 throw InvalidValueException(thres, "thres: 0 <= thres <= 1");
01323 }
01324
01325 EMData * amps = get_fft_amplitude();
01326 vector<float> ampvector = amps->get_data_as_vector();
01327
01328 sort (ampvector.begin(), ampvector.end());
01329 int thresidx = int(thres * ampvector.size());
01330 float thresamp = ampvector[thresidx];
01331
01332 return thresamp;
01333 }
01334
01335 vector<Vec3i > find_region(EMData* image,const vector<Vec3i >& coords, const float value, vector<Vec3i >& region)
01336 {
01337 static vector<Vec3i> two_six_connected;
01338 if (two_six_connected.size() == 0) {
01339 for(int i = -1; i <= 1; ++i) {
01340 for(int j = -1; j <= 1; ++j) {
01341 for(int k = -1; k <= 1; ++k) {
01342 if ( j != 0 || i != 0 || k != 0) {
01343 two_six_connected.push_back(Vec3i(i,j,k));
01344 }
01345 }
01346 }
01347 }
01348 }
01349
01350 vector<Vec3i> ret;
01351 for(vector<Vec3i>::const_iterator it = two_six_connected.begin(); it != two_six_connected.end(); ++it ) {
01352 for(vector<Vec3i>::const_iterator it2 = coords.begin(); it2 != coords.end(); ++it2 ) {
01353 if (image->get_value_at((*it2)[0],(*it2)[1],(*it2)[2]) != value) throw;
01354 Vec3i c = (*it)+(*it2);
01355
01356 if ( c[0] < 0 || c[0] >= image->get_xsize()) continue;
01357 if ( c[1] < 0 || c[1] >= image->get_ysize()) continue;
01358 if ( c[2] < 0 || c[2] >= image->get_zsize()) continue;
01359
01360 if( image->get_value_at(c[0],c[1],c[2]) == value ) {
01361 if (find(ret.begin(),ret.end(),c) == ret.end()) {
01362 if (find(region.begin(),region.end(),c) == region.end()) {
01363 region.push_back(c);
01364 ret.push_back(c);
01365 }
01366 }
01367 }
01368 }
01369 }
01370 return ret;
01371 }
01372
01373 vector<Vec3i> EMData::mask_contig_region(const float& value, const Vec3i& seed) {
01374 Vec3i coord(seed[0],seed[1],seed[2]);
01375 vector<Vec3i> region;
01376 region.push_back(coord);
01377 vector<Vec3i> find_region_input = region;
01378 while (true) {
01379 vector<Vec3i> v = find_region(this,find_region_input, value, region);
01380 if (v.size() == 0 ) break;
01381 else find_region_input = v;
01382 }
01383 return region;
01384 }