Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

rsconvolution.cpp File Reference

#include "emdata.h"
#include <algorithm>

Include dependency graph for rsconvolution.cpp:

Include dependency graph

Go to the source code of this file.

Namespaces

namespace  EMAN

Defines

#define SWAP(a, b)   temp=(a);(a)=(b);(b)=temp;

Functions

float mult_internal (EMData &K, EMData &f, int kzmin, int kzmax, int kymin, int kymax, int kxmin, int kxmax, int iz, int iy, int ix)
float mult_circ (EMData &K, EMData &f, int kzmin, int kzmax, int kymin, int kymax, int kxmin, int kxmax, int nzf, int nyf, int nxf, int iz, int iy, int ix)
float select_kth_smallest (float *table, int n, int k)
float median (EMData &f, int nxk, int nyk, int nzk, kernel_shape myshape, int iz, int iy, int ix)


Define Documentation

#define SWAP a,
 )     temp=(a);(a)=(b);(b)=temp;
 

Definition at line 36 of file rsconvolution.cpp.


Function Documentation

float median EMData f,
int  nxk,
int  nyk,
int  nzk,
kernel_shape  myshape,
int  iz,
int  iy,
int  ix
[inline, static]
 

Definition at line 130 of file rsconvolution.cpp.

References EMAN::BLOCK, EMAN::CIRCULAR, EMAN::CROSS, EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), ImageDimensionException, and select_kth_smallest().

Referenced by EMAN::filt_median_().

00130                                                                                                                 {
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         }

float mult_circ EMData K,
EMData f,
int  kzmin,
int  kzmax,
int  kymin,
int  kymax,
int  kxmin,
int  kxmax,
int  nzf,
int  nyf,
int  nxf,
int  iz,
int  iy,
int  ix
[inline, static]
 

Definition at line 60 of file rsconvolution.cpp.

Referenced by EMAN::rsconvolution().

00062                                                                                                       {
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         }

float mult_internal EMData K,
EMData f,
int  kzmin,
int  kzmax,
int  kymin,
int  kymax,
int  kxmin,
int  kxmax,
int  iz,
int  iy,
int  ix
[inline, static]
 

Definition at line 43 of file rsconvolution.cpp.

Referenced by EMAN::rsconvolution().

00046                                                                                    {
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         }

float select_kth_smallest float *  table,
int  n,
int  k
[static]
 

Definition at line 87 of file rsconvolution.cpp.

References flag.

Referenced by median().

00087                                                               {
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         }


Generated on Tue Jun 11 13:41:29 2013 for EMAN2 by  doxygen 1.3.9.1