#include <processor.h>
Inheritance diagram for EMAN::BilateralProcessor:


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 | |
| Processor * | NEW () |
Static Public Attributes | |
| const string | NAME = "filter.bilateral" |
Bilateral processing does non-linear weighted averaging processing within a certain window.
| 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 3996 of file processor.h.
|
|
Get the descrition of this specific processor. This function must be overwritten by a subclass.
Implements EMAN::Processor. Definition at line 4005 of file processor.h. 04006 {
04007 return "Bilateral processing on 2D or 3D volume data. Bilateral processing does non-linear weighted averaging processing within a certain window. ";
04008 }
|
|
|
Get the processor's name. Each processor is identified by a unique name.
Implements EMAN::Processor. Definition at line 4000 of file processor.h. 04001 {
04002 return NAME;
04003 }
|
|
|
Get processor parameter information in a dictionary. Each parameter has one record in the dictionary. Each record contains its name, data-type, and description.
Reimplemented from EMAN::Processor. Definition at line 4015 of file processor.h. References EMAN::TypeDict::put(). 04016 {
04017 TypeDict d;
04018 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.");
04019 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.");
04020 d.put("niter", EMObject::INT, "how many times to apply this processing on your data.");
04021 d.put("half_width", EMObject::INT, "processing window size = (2 * half_widthh + 1) ^ 3.");
04022 return d;
04023 }
|
|
|
Definition at line 4010 of file processor.h. 04011 {
04012 return new BilateralProcessor();
04013 }
|
|
|
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.
Implements EMAN::Processor. Definition at line 3895 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(). 03896 {
03897 if (!image) {
03898 LOGWARN("NULL Image");
03899 return;
03900 }
03901
03902 float distance_sigma = params["distance_sigma"];
03903 float value_sigma = params["value_sigma"];
03904 int max_iter = params["niter"];
03905 int half_width = params["half_width"];
03906
03907 if (half_width < distance_sigma) {
03908 LOGWARN("localwidth(=%d) should be larger than distance_sigma=(%f)\n",
03909 half_width, distance_sigma);
03910 }
03911
03912 distance_sigma *= distance_sigma;
03913
03914 float image_sigma = image->get_attr("sigma");
03915 if (image_sigma > value_sigma) {
03916 LOGWARN("image sigma(=%f) should be smaller than value_sigma=(%f)\n",
03917 image_sigma, value_sigma);
03918 }
03919 value_sigma *= value_sigma;
03920
03921 int nx = image->get_xsize();
03922 int ny = image->get_ysize();
03923 int nz = image->get_zsize();
03924
03925 if(nz==1) { //for 2D image
03926 int width=nx, height=ny;
03927
03928 int i,j,m,n;
03929
03930 float tempfloat1,tempfloat2,tempfloat3;
03931 int index1,index2,index;
03932 int Iter;
03933 int tempint1,tempint3;
03934
03935 tempint1=width;
03936 tempint3=width+2*half_width;
03937
03938 float* mask=(float*)calloc((2*half_width+1)*(2*half_width+1),sizeof(float));
03939 float* OrgImg=(float*)calloc((2*half_width+width)*(2*half_width+height),sizeof(float));
03940 float* NewImg=image->get_data();
03941
03942 for(m=-(half_width);m<=half_width;m++)
03943 for(n=-(half_width);n<=half_width;n++) {
03944 index=(m+half_width)*(2*half_width+1)+(n+half_width);
03945 mask[index]=exp((float)(-(m*m+n*n)/distance_sigma/2.0));
03946 }
03947
03948 //printf("entering bilateral filtering process \n");
03949
03950 Iter=0;
03951 while(Iter<max_iter) {
03952 for(i=0;i<height;i++)
03953 for(j=0;j<width;j++) {
03954 index1=(i+half_width)*tempint3+(j+half_width);
03955 index2=i*tempint1+j;
03956 OrgImg[index1]=NewImg[index2];
03957 }
03958
03959 // Mirror Padding
03960 for(i=0;i<height;i++){
03961 for(j=0;j<half_width;j++) OrgImg[(i+half_width)*tempint3+(j)]=OrgImg[(i+half_width)*tempint3+(2*half_width-j)];
03962 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)];
03963 }
03964 for(i=0;i<half_width;i++){
03965 for(j=0;j<(width+2*half_width);j++) OrgImg[i*tempint3+j]=OrgImg[(2*half_width-i)*tempint3+j];
03966 for(j=0;j<(width+2*half_width);j++) OrgImg[(i+height+half_width)*tempint3+j]=OrgImg[(height+half_width-2-i)*tempint3+j];
03967 }
03968
03969 //printf("finish mirror padding process \n");
03970 //now mirror padding have been done
03971
03972 for(i=0;i<height;i++){
03973 //printf("now processing the %d th row \n",i);
03974 for(j=0;j<width;j++){
03975 tempfloat1=0.0; tempfloat2=0.0;
03976 for(m=-(half_width);m<=half_width;m++)
03977 for(n=-(half_width);n<=half_width;n++){
03978 index =(m+half_width)*(2*half_width+1)+(n+half_width);
03979 index1=(i+half_width)*tempint3+(j+half_width);
03980 index2=(i+half_width+m)*tempint3+(j+half_width+n);
03981 tempfloat3=(OrgImg[index1]-OrgImg[index2])*(OrgImg[index1]-OrgImg[index2]);
03982
03983 tempfloat3=mask[index]*(1.0f/(1+tempfloat3/value_sigma)); // Lorentz kernel
03984 //tempfloat3=mask[index]*exp(tempfloat3/Sigma2/(-2.0)); // Guassian kernel
03985 tempfloat1+=tempfloat3;
03986
03987 tempfloat2+=tempfloat3*OrgImg[(i+half_width+m)*tempint3+(j+half_width+n)];
03988 }
03989 NewImg[i*width+j]=tempfloat2/tempfloat1;
03990 }
03991 }
03992 Iter++;
03993 }
03994
03995 //printf("have finished %d th iteration\n ",Iter);
03996 // doneData();
03997 free(mask);
03998 free(OrgImg);
03999 // end of BilaFilter routine
04000
04001 }
04002 else { //3D case
04003 int width = nx;
04004 int height = ny;
04005 int slicenum = nz;
04006
04007 int slice_size = width * height;
04008 int new_width = width + 2 * half_width;
04009 int new_slice_size = (width + 2 * half_width) * (height + 2 * half_width);
04010
04011 int width1 = 2 * half_width + 1;
04012 int mask_size = width1 * width1;
04013 int old_img_size = (2 * half_width + width) * (2 * half_width + height);
04014
04015 int zstart = -half_width;
04016 int zend = -half_width;
04017 int is_3d = 0;
04018 if (nz > 1) {
04019 mask_size *= width1;
04020 old_img_size *= (2 * half_width + slicenum);
04021 zend = half_width;
04022 is_3d = 1;
04023 }
04024
04025 float *mask = (float *) calloc(mask_size, sizeof(float));
04026 float *old_img = (float *) calloc(old_img_size, sizeof(float));
04027
04028 float *new_img = image->get_data();
04029
04030 for (int p = zstart; p <= zend; p++) {
04031 int cur_p = (p + half_width) * (2 * half_width + 1) * (2 * half_width + 1);
04032
04033 for (int m = -half_width; m <= half_width; m++) {
04034 int cur_m = (m + half_width) * (2 * half_width + 1) + half_width;
04035
04036 for (int n = -half_width; n <= half_width; n++) {
04037 int l = cur_p + cur_m + n;
04038 mask[l] = exp((float) (-(m * m + n * n + p * p * is_3d) / distance_sigma / 2.0f));
04039 }
04040 }
04041 }
04042
04043 int iter = 0;
04044 while (iter < max_iter) {
04045 for (int k = 0; k < slicenum; k++) {
04046 size_t cur_k1 = (size_t)(k + half_width) * new_slice_size * is_3d;
04047 int cur_k2 = k * slice_size;
04048
04049 for (int i = 0; i < height; i++) {
04050 int cur_i1 = (i + half_width) * new_width;
04051 int cur_i2 = i * width;
04052
04053 for (int j = 0; j < width; j++) {
04054 size_t k1 = cur_k1 + cur_i1 + (j + half_width);
04055 int k2 = cur_k2 + cur_i2 + j;
04056 old_img[k1] = new_img[k2];
04057 }
04058 }
04059 }
04060
04061 for (int k = 0; k < slicenum; k++) {
04062 size_t cur_k = (k + half_width) * new_slice_size * is_3d;
04063
04064 for (int i = 0; i < height; i++) {
04065 int cur_i = (i + half_width) * new_width;
04066
04067 for (int j = 0; j < half_width; j++) {
04068 size_t k1 = cur_k + cur_i + j;
04069 size_t k2 = cur_k + cur_i + (2 * half_width - j);
04070 old_img[k1] = old_img[k2];
04071 }
04072
04073 for (int j = 0; j < half_width; j++) {
04074 size_t k1 = cur_k + cur_i + (width + half_width + j);
04075 size_t k2 = cur_k + cur_i + (width + half_width - j - 2);
04076 old_img[k1] = old_img[k2];
04077 }
04078 }
04079
04080
04081 for (int i = 0; i < half_width; i++) {
04082 int i2 = i * new_width;
04083 int i3 = (2 * half_width - i) * new_width;
04084 for (int j = 0; j < (width + 2 * half_width); j++) {
04085 size_t k1 = cur_k + i2 + j;
04086 size_t k2 = cur_k + i3 + j;
04087 old_img[k1] = old_img[k2];
04088 }
04089
04090 i2 = (height + half_width + i) * new_width;
04091 i3 = (height + half_width - 2 - i) * new_width;
04092 for (int j = 0; j < (width + 2 * half_width); j++) {
04093 size_t k1 = cur_k + i2 + j;
04094 size_t k2 = cur_k + i3 + j;
04095 old_img[k1] = old_img[k2];
04096 }
04097 }
04098 }
04099
04100 size_t idx;
04101 for (int k = 0; k < slicenum; k++) {
04102 size_t cur_k = (k + half_width) * new_slice_size;
04103
04104 for (int i = 0; i < height; i++) {
04105 int cur_i = (i + half_width) * new_width;
04106
04107 for (int j = 0; j < width; j++) {
04108 float f1 = 0;
04109 float f2 = 0;
04110 size_t k1 = cur_k + cur_i + (j + half_width);
04111
04112 for (int p = zstart; p <= zend; p++) {
04113 size_t cur_p1 = (p + half_width) * (2 * half_width + 1) * (2 * half_width + 1);
04114 size_t cur_p2 = (k + half_width + p) * new_slice_size;
04115
04116 for (int m = -half_width; m <= half_width; m++) {
04117 size_t cur_m1 = (m + half_width) * (2 * half_width + 1);
04118 size_t cur_m2 = cur_p2 + cur_i + m * new_width + j + half_width;
04119
04120 for (int n = -half_width; n <= half_width; n++) {
04121 size_t k = cur_p1 + cur_m1 + (n + half_width);
04122 size_t k2 = cur_m2 + n;
04123 float f3 = Util::square(old_img[k1] - old_img[k2]);
04124
04125 f3 = mask[k] * (1.0f / (1 + f3 / value_sigma));
04126 f1 += f3;
04127 size_t l1 = cur_m2 + n;
04128 f2 += f3 * old_img[l1];
04129 }
04130
04131 idx = (size_t)k * height * width + i * width + j;
04132 new_img[idx] = f2 / f1;
04133 }
04134 }
04135 }
04136 }
04137 }
04138 iter++;
04139 }
04140 if( mask ) {
04141 free(mask);
04142 mask = 0;
04143 }
04144
04145 if( old_img ) {
04146 free(old_img);
04147 old_img = 0;
04148 }
04149 }
04150
04151 image->update();
04152 }
|
|
|
Definition at line 147 of file processor.cpp. |
1.3.9.1