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

EMAN::BilateralProcessor Class Reference

Bilateral processing on 2D or 3D volume data. More...

#include <processor.h>

Inheritance diagram for EMAN::BilateralProcessor:

Inheritance graph
[legend]
Collaboration diagram for EMAN::BilateralProcessor:

Collaboration graph
[legend]
List of all members.

Public Member Functions

void process_inplace (EMData *image)
 To process an image in-place.
string get_name () const
 Get the processor's name.
string get_desc () const
 Get the descrition of this specific processor.
TypeDict get_param_types () const
 Get processor parameter information in a dictionary.

Static Public Member Functions

ProcessorNEW ()

Static Public Attributes

const string NAME = "bilateral"

Detailed Description

Bilateral processing on 2D or 3D volume data.

Bilateral processing does non-linear weighted averaging processing within a certain window.

Parameters:
distance_sigma means how large the voxel has impact on its neighbors in spatial domain. The larger it is, the more blurry the resulting image.
value_sigma eans how large the voxel has impact on its in range domain. The larger it is, the more blurry the resulting image.
niter how many times to apply this processing on your data.
half_width processing window size = (2 * half_widthh + 1) ^ 3.

Definition at line 3902 of file processor.h.


Member Function Documentation

string EMAN::BilateralProcessor::get_desc  )  const [inline, virtual]
 

Get the descrition of this specific processor.

This function must be overwritten by a subclass.

Returns:
The description of this processor.

Implements EMAN::Processor.

Definition at line 3911 of file processor.h.

03912                 {
03913                         return "Bilateral processing on 2D or 3D volume data. Bilateral processing does non-linear weighted averaging processing within a certain window. ";
03914                 }

string EMAN::BilateralProcessor::get_name  )  const [inline, virtual]
 

Get the processor's name.

Each processor is identified by a unique name.

Returns:
The processor's name.

Implements EMAN::Processor.

Definition at line 3906 of file processor.h.

03907                 {
03908                         return NAME;
03909                 }

TypeDict EMAN::BilateralProcessor::get_param_types  )  const [inline, virtual]
 

Get processor parameter information in a dictionary.

Each parameter has one record in the dictionary. Each record contains its name, data-type, and description.

Returns:
A dictionary containing the parameter info.

Reimplemented from EMAN::Processor.

Definition at line 3921 of file processor.h.

References EMAN::TypeDict::put().

03922                 {
03923                         TypeDict d;
03924                         d.put("distance_sigma", EMObject::FLOAT, "means how large the voxel has impact on its neighbors in spatial domain. The larger it is, the more blurry the resulting image.");
03925                         d.put("value_sigma", EMObject::FLOAT, "means how large the voxel has impact on its in  range domain. The larger it is, the more blurry the resulting image.");
03926                         d.put("niter", EMObject::INT, "how many times to apply this processing on your data.");
03927                         d.put("half_width", EMObject::INT, "processing window size = (2 * half_widthh + 1) ^ 3.");
03928                         return d;
03929                 }

Processor* EMAN::BilateralProcessor::NEW  )  [inline, static]
 

Definition at line 3916 of file processor.h.

03917                 {
03918                         return new BilateralProcessor();
03919                 }

void BilateralProcessor::process_inplace EMData image  )  [virtual]
 

To process an image in-place.

For those processors which can only be processed out-of-place, override this function to just print out some error message to remind user call the out-of-place version.

Parameters:
image The image to be processed.

Implements EMAN::Processor.

Definition at line 3869 of file processor.cpp.

References EMAN::EMData::get_attr(), EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), LOGWARN, nx, ny, square, and EMAN::EMData::update().

03870 {
03871         if (!image) {
03872                 LOGWARN("NULL Image");
03873                 return;
03874         }
03875 
03876         float distance_sigma = params["distance_sigma"];
03877         float value_sigma = params["value_sigma"];
03878         int max_iter = params["niter"];
03879         int half_width = params["half_width"];
03880 
03881         if (half_width < distance_sigma) {
03882                 LOGWARN("localwidth(=%d) should be larger than distance_sigma=(%f)\n",
03883                                                         half_width, distance_sigma);
03884         }
03885 
03886         distance_sigma *= distance_sigma;
03887 
03888         float image_sigma = image->get_attr("sigma");
03889         if (image_sigma > value_sigma) {
03890                 LOGWARN("image sigma(=%f) should be smaller than value_sigma=(%f)\n",
03891                                                         image_sigma, value_sigma);
03892         }
03893         value_sigma *= value_sigma;
03894 
03895         int nx = image->get_xsize();
03896         int ny = image->get_ysize();
03897         int nz = image->get_zsize();
03898 
03899         if(nz==1) {     //for 2D image
03900                 int width=nx, height=ny;
03901 
03902                 int i,j,m,n;
03903 
03904                 float tempfloat1,tempfloat2,tempfloat3;
03905                 int   index1,index2,index;
03906                 int   Iter;
03907                 int   tempint1,tempint3;
03908 
03909                 tempint1=width;
03910                 tempint3=width+2*half_width;
03911 
03912                 float* mask=(float*)calloc((2*half_width+1)*(2*half_width+1),sizeof(float));
03913                 float* OrgImg=(float*)calloc((2*half_width+width)*(2*half_width+height),sizeof(float));
03914                 float* NewImg=image->get_data();
03915 
03916                 for(m=-(half_width);m<=half_width;m++)
03917                         for(n=-(half_width);n<=half_width;n++) {
03918                    index=(m+half_width)*(2*half_width+1)+(n+half_width);
03919                    mask[index]=exp((float)(-(m*m+n*n)/distance_sigma/2.0));
03920                 }
03921 
03922                 //printf("entering bilateral filtering process \n");
03923 
03924                 Iter=0;
03925                 while(Iter<max_iter) {
03926                         for(i=0;i<height;i++)
03927                         for(j=0;j<width;j++) {
03928                                 index1=(i+half_width)*tempint3+(j+half_width);
03929                                         index2=i*tempint1+j;
03930                                 OrgImg[index1]=NewImg[index2];
03931                         }
03932 
03933                         // Mirror Padding
03934                         for(i=0;i<height;i++){
03935                                 for(j=0;j<half_width;j++) OrgImg[(i+half_width)*tempint3+(j)]=OrgImg[(i+half_width)*tempint3+(2*half_width-j)];
03936                                 for(j=0;j<half_width;j++) OrgImg[(i+half_width)*tempint3+(j+width+half_width)]=OrgImg[(i+half_width)*tempint3+(width+half_width-j-2)];
03937                         }
03938                         for(i=0;i<half_width;i++){
03939                                 for(j=0;j<(width+2*half_width);j++) OrgImg[i*tempint3+j]=OrgImg[(2*half_width-i)*tempint3+j];
03940                                 for(j=0;j<(width+2*half_width);j++) OrgImg[(i+height+half_width)*tempint3+j]=OrgImg[(height+half_width-2-i)*tempint3+j];
03941                         }
03942 
03943                         //printf("finish mirror padding process \n");
03944                         //now mirror padding have been done
03945 
03946                         for(i=0;i<height;i++){
03947                                 //printf("now processing the %d th row \n",i);
03948                                 for(j=0;j<width;j++){
03949                                         tempfloat1=0.0; tempfloat2=0.0;
03950                                         for(m=-(half_width);m<=half_width;m++)
03951                                                 for(n=-(half_width);n<=half_width;n++){
03952                                                         index =(m+half_width)*(2*half_width+1)+(n+half_width);
03953                                                         index1=(i+half_width)*tempint3+(j+half_width);
03954                                                         index2=(i+half_width+m)*tempint3+(j+half_width+n);
03955                                                         tempfloat3=(OrgImg[index1]-OrgImg[index2])*(OrgImg[index1]-OrgImg[index2]);
03956 
03957                                                         tempfloat3=mask[index]*(1.0f/(1+tempfloat3/value_sigma));       // Lorentz kernel
03958                                                         //tempfloat3=mask[index]*exp(tempfloat3/Sigma2/(-2.0)); // Guassian kernel
03959                                                         tempfloat1+=tempfloat3;
03960 
03961                                                         tempfloat2+=tempfloat3*OrgImg[(i+half_width+m)*tempint3+(j+half_width+n)];
03962                                         }
03963                                         NewImg[i*width+j]=tempfloat2/tempfloat1;
03964                                 }
03965                         }
03966                         Iter++;
03967             }
03968 
03969             //printf("have finished %d  th iteration\n ",Iter);
03970 //              doneData();
03971                 free(mask);
03972                 free(OrgImg);
03973                 // end of BilaFilter routine
03974 
03975         }
03976         else {  //3D case
03977                 int width = nx;
03978                 int height = ny;
03979                 int slicenum = nz;
03980 
03981                 int slice_size = width * height;
03982                 int new_width = width + 2 * half_width;
03983                 int new_slice_size = (width + 2 * half_width) * (height + 2 * half_width);
03984 
03985                 int width1 = 2 * half_width + 1;
03986                 int mask_size = width1 * width1;
03987                 int old_img_size = (2 * half_width + width) * (2 * half_width + height);
03988 
03989                 int zstart = -half_width;
03990                 int zend = -half_width;
03991                 int is_3d = 0;
03992                 if (nz > 1) {
03993                         mask_size *= width1;
03994                         old_img_size *= (2 * half_width + slicenum);
03995                         zend = half_width;
03996                         is_3d = 1;
03997                 }
03998 
03999                 float *mask = (float *) calloc(mask_size, sizeof(float));
04000                 float *old_img = (float *) calloc(old_img_size, sizeof(float));
04001 
04002                 float *new_img = image->get_data();
04003 
04004                 for (int p = zstart; p <= zend; p++) {
04005                         int cur_p = (p + half_width) * (2 * half_width + 1) * (2 * half_width + 1);
04006 
04007                         for (int m = -half_width; m <= half_width; m++) {
04008                                 int cur_m = (m + half_width) * (2 * half_width + 1) + half_width;
04009 
04010                                 for (int n = -half_width; n <= half_width; n++) {
04011                                         int l = cur_p + cur_m + n;
04012                                         mask[l] = exp((float) (-(m * m + n * n + p * p * is_3d) / distance_sigma / 2.0f));
04013                                 }
04014                         }
04015                 }
04016 
04017                 int iter = 0;
04018                 while (iter < max_iter) {
04019                         for (int k = 0; k < slicenum; k++) {
04020                                 size_t cur_k1 = (size_t)(k + half_width) * new_slice_size * is_3d;
04021                                 int cur_k2 = k * slice_size;
04022 
04023                                 for (int i = 0; i < height; i++) {
04024                                         int cur_i1 = (i + half_width) * new_width;
04025                                         int cur_i2 = i * width;
04026 
04027                                         for (int j = 0; j < width; j++) {
04028                                                 size_t k1 = cur_k1 + cur_i1 + (j + half_width);
04029                                                 int k2 = cur_k2 + cur_i2 + j;
04030                                                 old_img[k1] = new_img[k2];
04031                                         }
04032                                 }
04033                         }
04034 
04035                         for (int k = 0; k < slicenum; k++) {
04036                                 size_t cur_k = (k + half_width) * new_slice_size * is_3d;
04037 
04038                                 for (int i = 0; i < height; i++) {
04039                                         int cur_i = (i + half_width) * new_width;
04040 
04041                                         for (int j = 0; j < half_width; j++) {
04042                                                 size_t k1 = cur_k + cur_i + j;
04043                                                 size_t k2 = cur_k + cur_i + (2 * half_width - j);
04044                                                 old_img[k1] = old_img[k2];
04045                                         }
04046 
04047                                         for (int j = 0; j < half_width; j++) {
04048                                                 size_t k1 = cur_k + cur_i + (width + half_width + j);
04049                                                 size_t k2 = cur_k + cur_i + (width + half_width - j - 2);
04050                                                 old_img[k1] = old_img[k2];
04051                                         }
04052                                 }
04053 
04054 
04055                                 for (int i = 0; i < half_width; i++) {
04056                                         int i2 = i * new_width;
04057                                         int i3 = (2 * half_width - i) * new_width;
04058                                         for (int j = 0; j < (width + 2 * half_width); j++) {
04059                                                 size_t k1 = cur_k + i2 + j;
04060                                                 size_t k2 = cur_k + i3 + j;
04061                                                 old_img[k1] = old_img[k2];
04062                                         }
04063 
04064                                         i2 = (height + half_width + i) * new_width;
04065                                         i3 = (height + half_width - 2 - i) * new_width;
04066                                         for (int j = 0; j < (width + 2 * half_width); j++) {
04067                                                 size_t k1 = cur_k + i2 + j;
04068                                                 size_t k2 = cur_k + i3 + j;
04069                                                 old_img[k1] = old_img[k2];
04070                                         }
04071                                 }
04072                         }
04073 
04074                         size_t idx;
04075                         for (int k = 0; k < slicenum; k++) {
04076                                 size_t cur_k = (k + half_width) * new_slice_size;
04077 
04078                                 for (int i = 0; i < height; i++) {
04079                                         int cur_i = (i + half_width) * new_width;
04080 
04081                                         for (int j = 0; j < width; j++) {
04082                                                 float f1 = 0;
04083                                                 float f2 = 0;
04084                                                 size_t k1 = cur_k + cur_i + (j + half_width);
04085 
04086                                                 for (int p = zstart; p <= zend; p++) {
04087                                                         size_t cur_p1 = (p + half_width) * (2 * half_width + 1) * (2 * half_width + 1);
04088                                                         size_t cur_p2 = (k + half_width + p) * new_slice_size;
04089 
04090                                                         for (int m = -half_width; m <= half_width; m++) {
04091                                                                 size_t cur_m1 = (m + half_width) * (2 * half_width + 1);
04092                                                                 size_t cur_m2 = cur_p2 + cur_i + m * new_width + j + half_width;
04093 
04094                                                                 for (int n = -half_width; n <= half_width; n++) {
04095                                                                         size_t k = cur_p1 + cur_m1 + (n + half_width);
04096                                                                         size_t k2 = cur_m2 + n;
04097                                                                         float f3 = Util::square(old_img[k1] - old_img[k2]);
04098 
04099                                                                         f3 = mask[k] * (1.0f / (1 + f3 / value_sigma));
04100                                                                         f1 += f3;
04101                                                                         size_t l1 = cur_m2 + n;
04102                                                                         f2 += f3 * old_img[l1];
04103                                                                 }
04104 
04105                                                                 idx = (size_t)k * height * width + i * width + j;
04106                                                                 new_img[idx] = f2 / f1;
04107                                                         }
04108                                                 }
04109                                         }
04110                                 }
04111                         }
04112                         iter++;
04113                 }
04114                 if( mask ) {
04115                         free(mask);
04116                         mask = 0;
04117                 }
04118 
04119                 if( old_img ) {
04120                         free(old_img);
04121                         old_img = 0;
04122                 }
04123         }
04124 
04125         image->update();
04126 }


Member Data Documentation

const string BilateralProcessor::NAME = "bilateral" [static]
 

Definition at line 139 of file processor.cpp.


The documentation for this class was generated from the following files:
Generated on Mon Mar 7 18:20:59 2011 for EMAN2 by  doxygen 1.3.9.1