EMAN2
rsconvolution.cpp
Go to the documentation of this file.
00001 /*
00002  * Author: Pawel A.Penczek, 09/09/2006 (Pawel.A.Penczek@uth.tmc.edu)
00003  * Copyright (c) 2000-2006 The University of Texas - Houston Medical School
00004  *
00005  * This software is issued under a joint BSD/GNU license. You may use the
00006  * source code in this file under either license. However, note that the
00007  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00008  * so you are responsible for compliance with the licenses of these packages
00009  * if you opt to use BSD licensing. The warranty disclaimer below holds
00010  * in either instance.
00011  *
00012  * This complete copyright notice must be included in any revised version of the
00013  * source code. Additional authorship citations may be added, but existing
00014  * author citations must be preserved.
00015  *
00016  * This program is free software; you can redistribute it and/or modify
00017  * it under the terms of the GNU General Public License as published by
00018  * the Free Software Foundation; either version 2 of the License, or
00019  * (at your option) any later version.
00020  *
00021  * This program is distributed in the hope that it will be useful,
00022  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00023  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00024  * GNU General Public License for more details.
00025  *
00026  * You should have received a copy of the GNU General Public License
00027  * along with this program; if not, write to the Free Software
00028  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
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         // K(i,j,k)*f(a-i, b-j, c-k) <-- internal, so no boundary condition issues
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         // K(i,j,k)*f(a-i, b-j, c-k) <-- Circulant boundary conditions
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         // In the future we may want to add other boundary conditions here
00082         
00083 
00084         // This function selects the k-th smallest element in the array table[0..n-1]
00085         // i.e. k = 1 gets the smallest element
00086         //      k = n gets the largest element
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;  // This is necessary because the index of array begins from 0
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) {//Does not work properly in 3D, corners are not done, PAP 07/16/09
00250                 // Kernel should be the smaller image
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                         // whoops, f smaller than K
00255                         swap(f,K); swap(nxf,nxK); swap(nyf,nyK); swap(nzf,nzK);
00256                 } else if ((nxK<=nxf)&&(nyK<=nyf)&&(nzK<=nzf)) {
00257                         // that's what it should be, so do nothing
00258                         ;
00259                 } else {
00260                         // incommensurate sizes
00261                         throw ImageDimensionException("input images are incommensurate");
00262                 }
00263                 // Kernel needs to be _odd_ in size
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                 // kernel corners, need to check for degenerate case
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                 // interior boundaries, need to check for degenerate cases
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                 // interior (no boundary condition issues here)
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                 // corners
00303                 // corner sizes, with checking for degenerate cases
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                 // corner starting locations, with checking for degenerate cases
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                 // corners
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                 // remaining stripes -- should use a more elegant (non-3D-specific) method here
00329                 // ix < ixmin
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                 // ix > ixmax
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                 // iy < iymin
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                 // iy > iymax
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                 // iz < izmin
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                 // iz > izmax
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                 // ix < ixmin, iy < iymin
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                 // ix < ixmin, iy > iymax
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         // ix > ixmax, iy < iymin
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                 // ix > ixmax, iy > iymax
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         // ix < ixmin, iz < izmin
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                  // ix < ixmin, iz > izmax
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          // ix > ixmin, iz < izmin
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                  // ix > ixmin, iz > izmax
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        // iy < iymin, iz < izmin
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        // iy < iymin, iz > izmax
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                 // iy > iymax, iz < izmin
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        // iy > iymax, iz > izmax
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                         // Kernel should be smaller than the size of image
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                         // Kernel needs to be odd in size
00572                         throw ImageDimensionException("Real-space kernel must have odd size so that the center is well-defined.");
00573                 }
00574 
00575                 if ( myshape == CIRCULAR ) {
00576                         // For CIRCULAR kernal, size must be same on all dimensions
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                         // whoops, f smaller than K
00609                         swap(f,K); swap(nxf,nxk); swap(nyf,nyk); swap(nzf,nzk);
00610                 } else if ( nxk > nxf || nyk > nyf || nzk > nzf ) {
00611                         // Incommensurate sizes
00612                         throw ImageDimensionException("Two input images are incommensurate.");
00613                 }
00614 
00615                 if ( nxk % 2 != 1 || nyk % 2 != 1 || nzk % 2 != 1 ) {
00616                         // Kernel needs to be odd in size
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                         // Check whether two images are truly binary.
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 //                                      int kzmin = iz-nzk2 < 0     ?   0   : iz-nzk2 ;
00656 //                                      int kzmax = iz+nzk2 > nzf-1 ? nzf-1 : iz+nzk2 ;
00657 //                                      int kymin = iy-nyk2 < 0     ?   0   : iy-nyk2 ;
00658 //                                      int kymax = iy+nyk2 > nyf-1 ? nyf-1 : iy+nyk2 ;
00659 //                                      int kxmin = ix-nxk2 < 0     ?   0   : ix-nxk2 ;
00660 //                                      int kxmax = ix+nxk2 > nxf-1 ? nxf-1 : ix+nxk2 ;
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                         // whoops, f smaller than K
00716                         swap(f,K); swap(nxf,nxk); swap(nyf,nyk); swap(nzf,nzk);
00717                 } else if ( nxk > nxf || nyk > nyf || nzk > nzf ) {
00718                         // Incommensurate sizes
00719                         throw ImageDimensionException("Two input images are incommensurate.");
00720                 }
00721 
00722                 if ( nxk % 2 != 1 || nyk % 2 != 1 || nzk % 2 != 1 ) {
00723                         // Kernel needs to be odd in size
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                         // Check whether two images are truly binary.
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 namespace {
00809         // K(i,j,k)*f(a-i, b-j, c-k) <-- internal, so no boundary condition issues
00810         inline float kmlt_internal(EMData& K, EMData& f, 
00811                                                        int kzmin, int kzmax, int kymin, int kymax, 
00812                                                            int kxmin, int kxmax, 
00813                                                            int iz, int iy, int ix) {
00814                 float sum = 0.f;
00815                 for (int kz = kzmin; kz <= kzmax; kz++) {
00816                         for (int ky = kymin; ky <= kymax; ky++) {
00817                                 for (int kx = kxmin; kx <= kxmax; kx++) {
00818                                         float Kp = K(kx,ky,kz);
00819                                         float fp = f(ix-kx,iy-ky,iz-kz);
00820                                         sum += Kp*fp;
00821                                 }
00822                         }
00823                 }
00824                 return sum;
00825         }
00826         // K(i,j,k)*f(a-i, b-j, c-k) <-- Circulant boundary conditions
00827         inline float kmlt_circ(EMData& K, EMData& f, int kzmin, int kzmax,
00828                                        int kymin, int kymax, int kxmin, int kxmax,
00829                                                    int nzf, int nyf, int nxf, int iz, int iy, int ix) {
00830                 float sum = 0.f;
00831                 for (int kz = kzmin; kz <= kzmax; kz++) {
00832                         int jz = (iz - kz) % nzf;
00833                         if (jz < 0) jz += nzf;
00834                         for (int ky = kymin; ky <= kymax; ky++) {
00835                                 int jy = (iy - ky) % nyf; 
00836                                 if (jy < 0) jy += nyf;
00837                                 for (int kx = kxmin; kx <= kxmax; kx++) {
00838                                         int jx = (ix - kx) % nxf; 
00839                                         if (jx < 0) jx += nxf;
00840                                         float Kp = K(kx,ky,kz);
00841                                         float fp = f(jx,jy,jz);
00842                                         sum += Kp*fp;
00843                                 }
00844                         }
00845                 }
00846                 return sum;
00847         }
00848         // In the future we may want to add other boundary conditions here
00849 }
00850 namespace EMAN {
00851 
00852     EMData* rscp(EMData* f) {
00853                 // Kernel should be the smaller image
00854                 int nxf=f->get_xsize(); int nyf=f->get_ysize(); int nzf=f->get_zsize();
00855                 const int npad = 2;
00856                 const int m = Util::get_min(nxf,nyf,nzf);
00857                 const int n = m*npad;
00858 
00859                 const int K = 6;  //params["kb_K"];
00860                 const float alpha = 1.75;  //params["kb_alpha"];
00861                 Util::KaiserBessel kb(alpha, K, m/2,K/(2.*n),n);
00862 
00863                 int nxK = K/2+1; nyK=nxK; nzK=nxK;
00864 
00865                 EMData* result = new EMData();
00866                 result->set_size(nxf, nyf, nzf);
00867                 result->to_zero();
00868                 // kernel corners, need to check for degenerate case
00869                 int kxmin = -nxK/2; int kymin = -nyK/2; int kzmin = -nzK/2;
00870                 int kxmax = (1 == nxK % 2) ? -kxmin : -kxmin - 1;
00871                 int kymax = (1 == nyK % 2) ? -kymin : -kymin - 1;
00872                 int kzmax = (1 == nzK % 2) ? -kzmin : -kzmin - 1;
00873                 // interior boundaries, need to check for degenerate cases
00874                 int izmin = 0, izmax = 0, iymin = 0, iymax = 0, ixmin = 0, ixmax = 0;
00875                 if (1 != nzf) {
00876                         izmin = -kzmin;
00877                         izmax = nzf - 1 - kzmax;
00878                 }
00879                 if (1 != nyf) {
00880                         iymin = -kymin;
00881                         iymax = nyf - 1 - kymax;
00882                 }
00883                 if (1 != nxf) {
00884                         ixmin = -kxmin;
00885                         ixmax = nxf - 1 - kxmax;
00886                 }
00887                 // interior (no boundary condition issues here)
00888                 for (int iz = izmin; iz <= izmax; iz++) {
00889                         for (int iy = iymin; iy <= iymax; iy++) {
00890                                 for (int ix = ixmin; ix <= ixmax; ix++) {
00891                                         (*result)(ix,iy,iz) =
00892                                                 mult_internal(*K, *f, 
00893                                                                       kzmin, kzmax, kymin, kymax, kxmin, kxmax,
00894                                                                           iz, iy, ix);
00895                                 }
00896                         }
00897                 }
00898                 //   INITIALLY SKIP IT / corners  
00899                 // corner sizes, with checking for degenerate cases
00900                 int sz = (1 == nzK) ? 1 : -kzmin + kzmax;
00901                 int sy = (1 == nyK) ? 1 : -kymin + kymax;
00902                 int sx = (1 == nxK) ? 1 : -kxmin + kxmax;
00903                 // corner starting locations, with checking for degenerate cases
00904                 int zstart = (0 == izmin) ? 0 : izmin - 1;
00905                 int ystart = (0 == iymin) ? 0 : iymin - 1;
00906                 int xstart = (0 == ixmin) ? 0 : ixmin - 1;
00907                 // corners
00908                 for (int cz = 0; cz < sz; cz++) {
00909                         int iz = (zstart - cz) % nzf;
00910                         if (iz < 0) iz += nzf;
00911                         for (int cy = 0; cy < sy; cy++) {
00912                                 int iy = (ystart - cy) % nyf;
00913                                 if (iy < 0) iy += nyf;
00914                                 for (int cx=0; cx < sx; cx++) {
00915                                         int ix = (xstart - cx) % nxf;
00916                                         if (ix < 0) ix += nxf;
00917                                         (*result)(ix,iy,iz) =
00918                                                 mult_circ(*K, *f, kzmin, kzmax, kymin, 
00919                                                                  kymax, kxmin, kxmax,
00920                                                                  nzf, nyf, nxf, iz, iy, ix);
00921                                 }
00922                         }
00923                 }
00924                 // remaining stripes -- should use a more elegant (non-3D-specific) method here
00925                 // ix < ixmin
00926                 for (int ix = 0; ix < ixmin; ix++) {
00927                         for (int iy = iymin; iy <= iymax; iy++) {
00928                                 for (int iz = izmin; iz <= izmax; iz++) {
00929                                         (*result)(ix,iy,iz) =
00930                                                 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax, 
00931                                                                  kxmin, kxmax,
00932                                                                  nzf, nyf, nxf, iz, iy, ix);
00933                                 }
00934                         }
00935                 }
00936                 // ix > ixmax
00937                 for (int ix = ixmax+1; ix < nxf; ix++) {
00938                         for (int iy = iymin; iy <= iymax; iy++) {
00939                                 for (int iz = izmin; iz <= izmax; iz++) {
00940                                         (*result)(ix,iy,iz) =
00941                                                 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax, 
00942                                                                  kxmin, kxmax,
00943                                                                  nzf, nyf, nxf, iz, iy, ix);
00944                                 }
00945                         }
00946                 }
00947                 // iy < iymin
00948                 for (int iy = 0; iy < iymin; iy++) {
00949                         for (int ix = ixmin; ix <= ixmax; ix++) {
00950                                 for (int iz = izmin; iz <= izmax; iz++) {
00951                                         (*result)(ix,iy,iz) =
00952                                                 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax, 
00953                                                                  kxmin, kxmax,
00954                                                                  nzf, nyf, nxf, iz, iy, ix);
00955                                 }
00956                         }
00957                 }
00958                 // iy > iymax
00959                 for (int iy = iymax+1; iy < nyf; iy++) {
00960                         for (int ix = ixmin; ix <= ixmax; ix++) {
00961                                 for (int iz = izmin; iz <= izmax; iz++) {
00962                                         (*result)(ix,iy,iz) =
00963                                                 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax, 
00964                                                                  kxmin, kxmax,
00965                                                                  nzf, nyf, nxf, iz, iy, ix);
00966                                 }
00967                         }
00968                 }
00969                 // iz < izmin
00970                 for (int iz = 0; iz < izmin; iz++) {
00971                         for (int ix = ixmin; ix <= ixmax; ix++) {
00972                                 for (int iy = iymin; iy <= iymax; iy++) {
00973                                         (*result)(ix,iy,iz) =
00974                                                 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax, 
00975                                                                  kxmin, kxmax,
00976                                                                  nzf, nyf, nxf, iz, iy, ix);
00977                                 }
00978                         }
00979                 }
00980                 // iz > izmax
00981                 for (int iz = izmax+1; iz < nzf; iz++) {
00982                         for (int ix = ixmin; ix <= ixmax; ix++) {
00983                                 for (int iy = iymin; iy <= iymax; iy++) {
00984                                         (*result)(ix,iy,iz) =
00985                                                 mult_circ(*K, *f, kzmin, kzmax, kymin, kymax, 
00986                                                                  kxmin, kxmax,
00987                                                                  nzf, nyf, nxf, iz, iy, ix);
00988                                 }
00989                         }
00990                 }
00991                 //K->set_array_offsets(K_saved_offsets);
00992                 result->update();
00993                 return result;
00994         }
00995 
00996 }
00997 */        
00998 /* vim: set ts=4 noet: */