00001
00002
00003
00004
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 #include "emdata.h"
00033
00034 #include <algorithm>
00035
00036 #define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
00037
00038 using namespace EMAN;
00039 using std::swap;
00040
00041 namespace {
00042
00043 inline float mult_internal(EMData& K, EMData& f,
00044 int kzmin, int kzmax, int kymin, int kymax,
00045 int kxmin, int kxmax,
00046 int iz, int iy, int ix) {
00047 float sum = 0.f;
00048 for (int kz = kzmin; kz <= kzmax; kz++) {
00049 for (int ky = kymin; ky <= kymax; ky++) {
00050 for (int kx = kxmin; kx <= kxmax; kx++) {
00051 float Kp = K(kx,ky,kz);
00052 float fp = f(ix-kx,iy-ky,iz-kz);
00053 sum += Kp*fp;
00054 }
00055 }
00056 }
00057 return sum;
00058 }
00059
00060 inline float mult_circ(EMData& K, EMData& f, int kzmin, int kzmax,
00061 int kymin, int kymax, int kxmin, int kxmax,
00062 int nzf, int nyf, int nxf, int iz, int iy, int ix) {
00063 float sum = 0.f;
00064 for (int kz = kzmin; kz <= kzmax; kz++) {
00065 int jz = (iz - kz) % nzf;
00066 if (jz < 0) jz += nzf;
00067 for (int ky = kymin; ky <= kymax; ky++) {
00068 int jy = (iy - ky) % nyf;
00069 if (jy < 0) jy += nyf;
00070 for (int kx = kxmin; kx <= kxmax; kx++) {
00071 int jx = (ix - kx) % nxf;
00072 if (jx < 0) jx += nxf;
00073 float Kp = K(kx,ky,kz);
00074 float fp = f(jx,jy,jz);
00075 sum += Kp*fp;
00076 }
00077 }
00078 }
00079 return sum;
00080 }
00081
00082
00083
00084
00085
00086
00087 float select_kth_smallest(float *table, int n, int k) {
00088
00089 int i,j,left,middle,right;
00090 float temp;
00091 bool flag = 0;
00092
00093 left = 0;
00094 right = n-1;
00095 k = k-1;
00096 while (flag == 0) {
00097 if ( left+1 < right ) {
00098 middle = (left+right)/2;
00099 swap(table[middle],table[left+1]);
00100 if ( table[left+1] > table [right] )
00101 swap(table[left+1], table[right]);
00102 if ( table[left] > table[right] )
00103 swap(table[left], table[right]);
00104 if ( table[left+1] > table[left] )
00105 swap(table[left+1], table[left]);
00106 i = left+1;
00107 j = right;
00108 temp = table[left];
00109 do {
00110 i++;
00111 while (table[i] < temp) i++;
00112 j--;
00113 while (table[j] > temp) j--;
00114 if (j >= i)
00115 swap(table[i], table[j]);
00116 } while (j >= i);
00117 table[left] = table[j];
00118 table[j] = temp;
00119 if (j >= k) right = j-1;
00120 if (j <= k) left = i;
00121 } else {
00122 if ( right == left+1 && table[right] < table[left] )
00123 swap(table[left], table[right]);
00124 flag = 1;
00125 }
00126 }
00127 return table[k];
00128 }
00129
00130 inline float median(EMData& f, int nxk, int nyk, int nzk, kernel_shape myshape, int iz, int iy, int ix) {
00131 size_t index = 0;
00132 int dimension = 3;
00133 float median_value = 0.f;
00134 float *table = 0;
00135
00136 int nxf = (&f)->get_xsize();
00137 int nyf = (&f)->get_ysize();
00138 int nzf = (&f)->get_zsize();
00139
00140 int nxk2 = (nxk-1)/2;
00141 int nyk2 = (nyk-1)/2;
00142 int nzk2 = (nzk-1)/2;
00143
00144 int kzmin = iz-nzk2;
00145 int kzmax = iz+nzk2;
00146 int kymin = iy-nyk2;
00147 int kymax = iy+nyk2;
00148 int kxmin = ix-nxk2;
00149 int kxmax = ix+nxk2;
00150
00151 if ( nzf == 1 ) {
00152 dimension--;
00153 if ( nyf == 1 ) dimension--;
00154 }
00155
00156 switch (myshape) {
00157 case BLOCK:
00158 switch (dimension) {
00159 case 1:
00160 table = (float*)malloc(nxk*sizeof(float));
00161 break;
00162 case 2: table = (float*)malloc(nxk*nyk*sizeof(float));
00163 break;
00164 case 3: table = (float*)malloc(nxk*nyk*nzk*sizeof(float));
00165 break;
00166 }
00167 for (int kz = kzmin; kz <= kzmax; kz++) {
00168 int jz = kz < 0 ? kz+nzf : kz % nzf;
00169 for (int ky = kymin; ky <= kymax; ky++) {
00170 int jy = ky < 0 ? ky+nyf : ky % nyf;
00171 for (int kx = kxmin; kx <= kxmax; kx++) {
00172 int jx = kx < 0 ? kx+nxf : kx % nxf;
00173 table[index] = f(jx,jy,jz);
00174 index++;
00175 }
00176 }
00177 }
00178 break;
00179 case CIRCULAR:
00180 switch (dimension) {
00181 case 1:
00182 table = (float*)malloc(nxk*sizeof(float));
00183 break;
00184 case 2: table = (float*)malloc(nxk*nxk*sizeof(float));
00185 break;
00186 case 3: table = (float*)malloc(nxk*nxk*nxk*sizeof(float));
00187 break;
00188 }
00189 for (int kz = kzmin; kz <= kzmax; kz++) {
00190 int jz = kz < 0 ? kz+nzf : kz % nzf;
00191 for (int ky = kymin; ky <= kymax; ky++) {
00192 int jy = ky < 0 ? ky+nyf : ky % nyf;
00193 for (int kx = kxmin; kx <= kxmax; kx++) {
00194 int jx = kx < 0 ? kx+nxf : kx % nxf;
00195 if ( (kz-iz)*(kz-iz)+(ky-iy)*(ky-iy)+(kx-ix)*(kx-ix) <= nxk2*nxk2 ) {
00196 table[index] = f(jx,jy,jz);
00197 index++;
00198 }
00199 }
00200 }
00201 }
00202 break;
00203 case CROSS:
00204 if ( nzf != 1 ) {
00205 table = (float*)malloc((nxk+nyk+nzk-2)*sizeof(float));
00206 for (int kz = kzmin; kz <= kzmax; kz++) {
00207 int jz = kz < 0 ? kz+nzf : kz % nzf;
00208 if ( kz != iz ) { table[index] = f(ix,iy,jz); index++; }
00209 }
00210 for (int ky = kymin; ky <= kymax; ky++) {
00211 int jy = ky < 0 ? ky+nyf : ky % nyf;
00212 if ( ky != iy ) { table[index] = f(ix,jy,iz); index++; }
00213 }
00214 for (int kx = kxmin; kx <= kxmax; kx++) {
00215 int jx = kx < 0 ? kx+nxf : kx % nxf;
00216 table[index] = f(jx,iy,iz);
00217 index++;
00218 }
00219 } else if ( nyf != 1 ) {
00220 table = (float*)malloc((nxk+nyk-1)*sizeof(float));
00221 for (int ky = kymin; ky <= kymax; ky++) {
00222 int jy = ky < 0 ? ky+nyf : ky % nyf;
00223 if ( ky != iy ) { table[index] = f(ix,jy,iz); index++; }
00224 }
00225 for (int kx = kxmin; kx <= kxmax; kx++) {
00226 int jx = kx < 0 ? kx+nxf : kx % nxf;
00227 table[index] = f(jx,iy,iz);
00228 index++;
00229 }
00230 } else {
00231 table = (float*)malloc(nxk*sizeof(float));
00232 for (int kx = kxmin; kx <= kxmax; kx++) {
00233 int jx = kx < 0 ? kx+nxf : kx % nxf;
00234 table[index] = f(jx,iy,iz);
00235 index++;
00236 }
00237 }
00238 break;
00239 default: throw ImageDimensionException("Illegal Kernal Shape!");
00240 }
00241 median_value=select_kth_smallest(table, index, (index+1)/2);
00242 free((void *)table);
00243 return median_value;
00244 }
00245 }
00246
00247 namespace EMAN {
00248
00249 EMData* rsconvolution(EMData* f, EMData* K) {
00250
00251 int nxf=f->get_xsize(); int nyf=f->get_ysize(); int nzf=f->get_zsize();
00252 int nxK=K->get_xsize(); int nyK=K->get_ysize(); int nzK=K->get_zsize();
00253 if ((nxf<nxK)&&(nyf<nyK)&&(nzf<nzK)) {
00254
00255 swap(f,K); swap(nxf,nxK); swap(nyf,nyK); swap(nzf,nzK);
00256 } else if ((nxK<=nxf)&&(nyK<=nyf)&&(nzK<=nzf)) {
00257
00258 ;
00259 } else {
00260
00261 throw ImageDimensionException("input images are incommensurate");
00262 }
00263
00264 if ((nxK % 2 != 1) || (nyK % 2 != 1) || (nzK % 2 != 1))
00265 throw ImageDimensionException("Real-space convolution kernel"
00266 " must have odd nx,ny,nz (so the center is well-defined).");
00267 EMData* result = new EMData();
00268 result->set_size(nxf, nyf, nzf);
00269 result->to_zero();
00270
00271 int kxmin = -nxK/2; int kymin = -nyK/2; int kzmin = -nzK/2;
00272 int kxmax = (1 == nxK % 2) ? -kxmin : -kxmin - 1;
00273 int kymax = (1 == nyK % 2) ? -kymin : -kymin - 1;
00274 int kzmax = (1 == nzK % 2) ? -kzmin : -kzmin - 1;
00275 vector<int> K_saved_offsets = K->get_array_offsets();
00276 K->set_array_offsets(kxmin,kymin,kzmin);
00277
00278 int izmin = 0, izmax = 0, iymin = 0, iymax = 0, ixmin = 0, ixmax = 0;
00279 if (1 != nzf) {
00280 izmin = -kzmin;
00281 izmax = nzf - 1 - kzmax;
00282 }
00283 if (1 != nyf) {
00284 iymin = -kymin;
00285 iymax = nyf - 1 - kymax;
00286 }
00287 if (1 != nxf) {
00288 ixmin = -kxmin;
00289 ixmax = nxf - 1 - kxmax;
00290 }
00291
00292 for (int iz = izmin; iz <= izmax; iz++) {
00293 for (int iy = iymin; iy <= iymax; iy++) {
00294 for (int ix = ixmin; ix <= ixmax; ix++) {
00295 (*result)(ix,iy,iz) =
00296 mult_internal(*K, *f,
00297 kzmin, kzmax, kymin, kymax, kxmin, kxmax,
00298 iz, iy, ix);
00299 }
00300 }
00301 }
00302
00303
00304 int sz = (1 == nzK) ? 1 : -kzmin + kzmax;
00305 int sy = (1 == nyK) ? 1 : -kymin + kymax;
00306 int sx = (1 == nxK) ? 1 : -kxmin + kxmax;
00307
00308 int zstart = (0 == izmin) ? 0 : izmin - 1;
00309 int ystart = (0 == iymin) ? 0 : iymin - 1;
00310 int xstart = (0 == ixmin) ? 0 : ixmin - 1;
00311
00312 for (int cz = 0; cz < sz; cz++) {
00313 int iz = (zstart - cz) % nzf;
00314 if (iz < 0) iz += nzf;
00315 for (int cy = 0; cy < sy; cy++) {
00316 int iy = (ystart - cy) % nyf;
00317 if (iy < 0) iy += nyf;
00318 for (int cx=0; cx < sx; cx++) {
00319 int ix = (xstart - cx) % nxf;
00320 if (ix < 0) ix += nxf;
00321 (*result)(ix,iy,iz) =
00322 mult_circ(*K, *f, kzmin, kzmax, kymin,
00323 kymax, kxmin, kxmax,
00324 nzf, nyf, nxf, iz, iy, ix);
00325 }
00326 }
00327 }
00328
00329
00330 for (int ix = 0; ix < ixmin; ix++) {
00331 for (int iy = iymin; iy <= iymax; iy++) {
00332 for (int iz = izmin; iz <= izmax; iz++) {
00333 (*result)(ix,iy,iz) =
00334 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax,
00335 kxmin, kxmax,
00336 nzf, nyf, nxf, iz, iy, ix);
00337 }
00338 }
00339 }
00340
00341 for (int ix = ixmax+1; ix < nxf; ix++) {
00342 for (int iy = iymin; iy <= iymax; iy++) {
00343 for (int iz = izmin; iz <= izmax; iz++) {
00344 (*result)(ix,iy,iz) =
00345 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax,
00346 kxmin, kxmax,
00347 nzf, nyf, nxf, iz, iy, ix);
00348 }
00349 }
00350 }
00351
00352 for (int iy = 0; iy < iymin; iy++) {
00353 for (int ix = ixmin; ix <= ixmax; ix++) {
00354 for (int iz = izmin; iz <= izmax; iz++) {
00355 (*result)(ix,iy,iz) =
00356 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax,
00357 kxmin, kxmax,
00358 nzf, nyf, nxf, iz, iy, ix);
00359 }
00360 }
00361 }
00362
00363 for (int iy = iymax+1; iy < nyf; iy++) {
00364 for (int ix = ixmin; ix <= ixmax; ix++) {
00365 for (int iz = izmin; iz <= izmax; iz++) {
00366 (*result)(ix,iy,iz) =
00367 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax,
00368 kxmin, kxmax,
00369 nzf, nyf, nxf, iz, iy, ix);
00370 }
00371 }
00372 }
00373
00374 for (int iz = 0; iz < izmin; iz++) {
00375 for (int ix = ixmin; ix <= ixmax; ix++) {
00376 for (int iy = iymin; iy <= iymax; iy++) {
00377 (*result)(ix,iy,iz) =
00378 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax,
00379 kxmin, kxmax,
00380 nzf, nyf, nxf, iz, iy, ix);
00381 }
00382 }
00383 }
00384
00385 for (int iz = izmax+1; iz < nzf; iz++) {
00386 for (int ix = ixmin; ix <= ixmax; ix++) {
00387 for (int iy = iymin; iy <= iymax; iy++) {
00388 (*result)(ix,iy,iz) =
00389 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax,
00390 kxmin, kxmax,
00391 nzf, nyf, nxf, iz, iy, ix);
00392 }
00393 }
00394 }
00395
00396
00397
00398 for (int ix = 0; ix < ixmin; ix++) {
00399 for (int iy = 0; iy < iymin; iy++) {
00400 for (int iz = izmin; iz <= izmax; iz++) {
00401 (*result)(ix,iy,iz) =
00402 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax,
00403 kxmin, kxmax,
00404 nzf, nyf, nxf, iz, iy, ix);
00405 }
00406 }
00407 }
00408
00409
00410 for (int ix = 0; ix < ixmin; ix++) {
00411 for (int iy = iymax+1; iy < nyf; iy++) {
00412 for (int iz = izmin; iz <= izmax; iz++) {
00413 (*result)(ix,iy,iz) =
00414 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax,
00415 kxmin, kxmax,
00416 nzf, nyf, nxf, iz, iy, ix);
00417 }
00418 }
00419 }
00420
00421
00422 for (int ix = ixmax+1; ix < nxf; ix++) {
00423 for (int iy = 0; iy < iymin; iy++) {
00424 for (int iz = izmin; iz <= izmax; iz++) {
00425 (*result)(ix,iy,iz) =
00426 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax,
00427 kxmin, kxmax,
00428 nzf, nyf, nxf, iz, iy, ix);
00429 }
00430 }
00431 }
00432
00433
00434 for (int ix = ixmax+1; ix < nxf; ix++) {
00435 for (int iy = iymax+1; iy < nyf; iy++) {
00436 for (int iz = izmin; iz <= izmax; iz++) {
00437 (*result)(ix,iy,iz) =
00438 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax,
00439 kxmin, kxmax,
00440 nzf, nyf, nxf, iz, iy, ix);
00441 }
00442 }
00443 }
00444
00445
00446
00447
00448 for (int ix = 0; ix < ixmin; ix++) {
00449 for (int iy = iymin; iy <= iymax; iy++) {
00450 for (int iz = 0; iz < izmin; iz++) {
00451 (*result)(ix,iy,iz) =
00452 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax,
00453 kxmin, kxmax,
00454 nzf, nyf, nxf, iz, iy, ix);
00455 }
00456 }
00457 }
00458
00459
00460 for (int ix = 0; ix < ixmin; ix++) {
00461 for (int iy = iymin; iy <= iymax; iy++) {
00462 for (int iz = izmax+1; iz < nzf; iz++) {
00463 (*result)(ix,iy,iz) =
00464 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax,
00465 kxmin, kxmax,
00466 nzf, nyf, nxf, iz, iy, ix);
00467 }
00468 }
00469 }
00470
00471
00472
00473 for (int ix = ixmax+1; ix < nxf; ix++) {
00474 for (int iy = iymin; iy <= iymax; iy++) {
00475 for (int iz = 0; iz < izmin; iz++) {
00476 (*result)(ix,iy,iz) =
00477 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax,
00478 kxmin, kxmax,
00479 nzf, nyf, nxf, iz, iy, ix);
00480 }
00481 }
00482 }
00483
00484
00485 for (int ix = ixmax+1; ix < nxf; ix++) {
00486 for (int iy = iymin; iy <= iymax; iy++) {
00487 for (int iz = izmax+1; iz < nzf; iz++) {
00488 (*result)(ix,iy,iz) =
00489 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax,
00490 kxmin, kxmax,
00491 nzf, nyf, nxf, iz, iy, ix);
00492 }
00493 }
00494 }
00495
00496
00497
00498
00499
00500 for (int iz = 0; iz < izmin; iz++) {
00501 for (int ix = ixmin; ix <= ixmax; ix++) {
00502 for (int iy = 0; iy < iymin; iy++) {
00503 (*result)(ix,iy,iz) =
00504 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax,
00505 kxmin, kxmax,
00506 nzf, nyf, nxf, iz, iy, ix);
00507 }
00508 }
00509 }
00510
00511
00512
00513
00514 for (int iz = izmax+1; iz < nzf; iz++) {
00515 for (int ix = ixmin; ix <= ixmax; ix++) {
00516 for (int iy = 0; iy < iymin; iy++) {
00517 (*result)(ix,iy,iz) =
00518 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax,
00519 kxmin, kxmax,
00520 nzf, nyf, nxf, iz, iy, ix);
00521 }
00522 }
00523 }
00524
00525
00526
00527
00528 for (int iz = 0; iz < izmin; iz++) {
00529 for (int ix = ixmin; ix <= ixmax; ix++) {
00530 for (int iy = iymax+1; iy < nyf; iy++) {
00531 (*result)(ix,iy,iz) =
00532 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax,
00533 kxmin, kxmax,
00534 nzf, nyf, nxf, iz, iy, ix);
00535 }
00536 }
00537 }
00538
00539
00540
00541
00542 for (int iz = izmax+1; iz < nzf; iz++) {
00543 for (int ix = ixmin; ix <= ixmax; ix++) {
00544 for (int iy = iymax+1; iy < nyf; iy++) {
00545 (*result)(ix,iy,iz) =
00546 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax,
00547 kxmin, kxmax,
00548 nzf, nyf, nxf, iz, iy, ix);
00549 }
00550 }
00551 }
00552
00553
00554 K->set_array_offsets(K_saved_offsets);
00555 result->update();
00556 return result;
00557 }
00558
00559 EMData* filt_median_(EMData* f, int nxk, int nyk, int nzk, kernel_shape myshape) {
00560
00561 int nxf = f->get_xsize();
00562 int nyf = f->get_ysize();
00563 int nzf = f->get_zsize();
00564
00565 if ( nxk > nxf || nyk > nyf || nzk > nzf ) {
00566
00567 throw ImageDimensionException("Kernel should be smaller than the size of image.");
00568 }
00569
00570 if ( nxk % 2 != 1 || nyk % 2 != 1 || nzk % 2 != 1 ) {
00571
00572 throw ImageDimensionException("Real-space kernel must have odd size so that the center is well-defined.");
00573 }
00574
00575 if ( myshape == CIRCULAR ) {
00576
00577 if ( (nzf != 1 && ( nxk != nyk || nxk != nzk )) || (nzf == 1 && nyf != 1 && nxk != nyk) ) {
00578 throw ImageDimensionException("For CIRCULAR kernal, size must be same on all dimensions.");
00579 }
00580 }
00581
00582 EMData* result = new EMData();
00583 result->set_size(nxf, nyf, nzf);
00584 result->to_zero();
00585
00586 for (int iz = 0; iz <= nzf-1; iz++) {
00587 for (int iy = 0; iy <= nyf-1; iy++) {
00588 for (int ix = 0; ix <= nxf-1; ix++) {
00589 (*result)(ix,iy,iz) = median (*f, nxk, nyk, nzk, myshape, iz, iy, ix);
00590 }
00591 }
00592 }
00593
00594 return result;
00595 }
00596
00597 EMData* filt_dilation_(EMData* f, EMData* K, morph_type mydilation) {
00598
00599 int nxf = f->get_xsize();
00600 int nyf = f->get_ysize();
00601 int nzf = f->get_zsize();
00602
00603 int nxk = K->get_xsize();
00604 int nyk = K->get_ysize();
00605 int nzk = K->get_zsize();
00606
00607 if ( nxf < nxk && nyf < nyk && nzf < nzk ) {
00608
00609 swap(f,K); swap(nxf,nxk); swap(nyf,nyk); swap(nzf,nzk);
00610 } else if ( nxk > nxf || nyk > nyf || nzk > nzf ) {
00611
00612 throw ImageDimensionException("Two input images are incommensurate.");
00613 }
00614
00615 if ( nxk % 2 != 1 || nyk % 2 != 1 || nzk % 2 != 1 ) {
00616
00617 throw ImageDimensionException("Kernel should have odd nx,ny,nz so that the center is well-defined.");
00618 }
00619
00620 int nxk2 = (nxk-1)/2;
00621 int nyk2 = (nyk-1)/2;
00622 int nzk2 = (nzk-1)/2;
00623
00624 if ( mydilation == BINARY ) {
00625
00626 for (int iz = 0; iz <= nzf-1; iz++) {
00627 for (int iy = 0; iy <= nyf-1; iy++) {
00628 for (int ix = 0; ix <= nxf-1; ix++) {
00629 int fxyz=(int)(*f)(ix,iy,iz);
00630 if ( fxyz != 0 && fxyz != 1 ) {
00631 throw ImageDimensionException("One of the two images is not binary.");
00632 }
00633 }
00634 }
00635 }
00636 for (int iz = 0; iz <= nzk-1; iz++) {
00637 for (int iy = 0; iy <= nyk-1; iy++) {
00638 for (int ix = 0; ix <= nxk-1; ix++) {
00639 int kxyz=(int)(*K)(ix,iy,iz);
00640 if ( kxyz != 0 && kxyz != 1 ) {
00641 throw ImageDimensionException("One of the two images is not binary.");
00642 }
00643 }
00644 }
00645 }
00646 }
00647
00648 EMData* result = new EMData();
00649 result->set_size(nxf, nyf, nzf);
00650 result->to_zero();
00651
00652 for (int iz = 0; iz <= nzf-1; iz++) {
00653 for (int iy = 0; iy <= nyf-1; iy++) {
00654 for (int ix = 0; ix <= nxf-1; ix++) {
00655
00656
00657
00658
00659
00660
00661 if ( mydilation == BINARY ) {
00662 int fxyz = (int)(*f)(ix,iy,iz);
00663 if ( fxyz == 1 ) {
00664 for (int jz = -nzk2; jz <= nzk2; jz++) {
00665 for (int jy = -nyk2; jy <= nyk2; jy++) {
00666 for (int jx= -nxk2; jx <= nxk2; jx++) {
00667 if ( (int)(*K)(jx+nxk2,jy+nyk2,jz+nzk2) == 1 ) {
00668 int fz = iz+jz;
00669 int fy = iy+jy;
00670 int fx = ix+jx;
00671 if ( fz >= 0 && fz <= nzf-1 && fy >= 0 && fy <= nyf-1 && fx >= 0 && fx <= nxf-1 )
00672 (*result)(fx,fy,fz) = 1;
00673 }
00674 }
00675 }
00676 }
00677 }
00678 } else if ( mydilation == GRAYLEVEL ) {
00679 float pmax = (*f)(ix,iy,iz)+(*K)(nxk2,nyk2,nzk2);
00680 for (int jz = -nzk2; jz <= nzk2; jz++) {
00681 for (int jy = -nyk2; jy <= nyk2; jy++) {
00682 for (int jx = -nxk2; jx <= nxk2; jx++) {
00683 int fz = iz+jz;
00684 int fy = iy+jy;
00685 int fx = ix+jx;
00686 if ( fz >= 0 && fz <= nzf-1 && fy >= 0 && fy <= nyf-1 && fx >= 0 && fx <= nxf-1 ) {
00687 float kxyz = (*K)(jx+nxk2,jy+nyk2,jz+nzk2);
00688 float fxyz = (*f)(fx,fy,fz);
00689 if ( kxyz+fxyz > pmax ) pmax = kxyz+fxyz;
00690 }
00691 }
00692 }
00693 }
00694 (*result)(ix,iy,iz) = pmax;
00695 } else {
00696 throw ImageDimensionException("Illegal dilation type!");
00697 }
00698 }
00699 }
00700 }
00701 return result;
00702 }
00703
00704 EMData* filt_erosion_(EMData* f, EMData* K, morph_type myerosion) {
00705
00706 int nxf = f->get_xsize();
00707 int nyf = f->get_ysize();
00708 int nzf = f->get_zsize();
00709
00710 int nxk = K->get_xsize();
00711 int nyk = K->get_ysize();
00712 int nzk = K->get_zsize();
00713
00714 if ( nxf < nxk && nyf < nyk && nzf < nzk ) {
00715
00716 swap(f,K); swap(nxf,nxk); swap(nyf,nyk); swap(nzf,nzk);
00717 } else if ( nxk > nxf || nyk > nyf || nzk > nzf ) {
00718
00719 throw ImageDimensionException("Two input images are incommensurate.");
00720 }
00721
00722 if ( nxk % 2 != 1 || nyk % 2 != 1 || nzk % 2 != 1 ) {
00723
00724 throw ImageDimensionException("Kernel should have odd nx,ny,nz so that the center is well-defined.");
00725 }
00726
00727 int nxk2 = (nxk-1)/2;
00728 int nyk2 = (nyk-1)/2;
00729 int nzk2 = (nzk-1)/2;
00730
00731 if ( myerosion == BINARY ) {
00732
00733 for (int iz = 0; iz <= nzf-1; iz++) {
00734 for (int iy = 0; iy <= nyf-1; iy++) {
00735 for (int ix = 0; ix <= nxf-1; ix++) {
00736 int fxyz=(int)(*f)(ix,iy,iz);
00737 if ( fxyz != 0 && fxyz != 1 ) {
00738 throw ImageDimensionException("One of the two images is not binary.");
00739 }
00740 }
00741 }
00742 }
00743 for (int iz = 0; iz <= nzk-1; iz++) {
00744 for (int iy = 0; iy <= nyk-1; iy++) {
00745 for (int ix = 0; ix <= nxk-1; ix++) {
00746 int kxyz=(int)(*K)(ix,iy,iz);
00747 if ( kxyz != 0 && kxyz != 1 ) {
00748 throw ImageDimensionException("One of the two images is not binary.");
00749 }
00750 }
00751 }
00752 }
00753 }
00754
00755 EMData* result = new EMData();
00756 result->set_size(nxf, nyf, nzf);
00757 result->to_one();
00758
00759 for (int iz = 0; iz <= nzf-1; iz++) {
00760 for (int iy = 0; iy <= nyf-1; iy++) {
00761 for (int ix = 0; ix <= nxf-1; ix++) {
00762 if ( myerosion == BINARY ) {
00763 int fxyz = (int)(*f)(ix,iy,iz);
00764 if ( fxyz == 0 ) {
00765 for (int jz = -nzk2; jz <= nzk2; jz++) {
00766 for (int jy = -nyk2; jy <= nyk2; jy++) {
00767 for (int jx= -nxk2; jx <= nxk2; jx++) {
00768 if ( (int)(*K)(jx+nxk2,jy+nyk2,jz+nzk2) == 1 ) {
00769 int fz = iz+jz;
00770 int fy = iy+jy;
00771 int fx = ix+jx;
00772 if ( fz >= 0 && fz <= nzf-1 && fy >= 0 && fy <= nyf-1 && fx >= 0 && fx <= nxf-1 )
00773 (*result)(fx,fy,fz) = 0;
00774 }
00775 }
00776 }
00777 }
00778 }
00779 } else if ( myerosion == GRAYLEVEL ) {
00780 float pmin = (*f)(ix,iy,iz)-(*K)(nxk2,nyk2,nzk2);
00781 for (int jz = -nzk2; jz <= nzk2; jz++) {
00782 for (int jy = -nyk2; jy <= nyk2; jy++) {
00783 for (int jx = -nxk2; jx <= nxk2; jx++) {
00784 int fz = iz+jz;
00785 int fy = iy+jy;
00786 int fx = ix+jx;
00787 if ( fz >= 0 && fz <= nzf-1 && fy >= 0 && fy <= nyf-1 && fx >= 0 && fx <= nxf-1 ) {
00788 float kxyz = (*K)(jx+nxk2,jy+nyk2,jz+nzk2);
00789 float fxyz = (*f)(fx,fy,fz);
00790 if ( fxyz-kxyz < pmin ) pmin = fxyz-kxyz;
00791 }
00792 }
00793 }
00794 }
00795 (*result)(ix,iy,iz) = pmin;
00796 } else {
00797 throw ImageDimensionException("Illegal dilation type!");
00798 }
00799 }
00800 }
00801 }
00802 return result;
00803 }
00804
00805 }
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998