#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 | |
static Processor * | NEW () |
Static Public Attributes | |
static 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 4072 of file processor.h.
string EMAN::BilateralProcessor::get_desc | ( | ) | const [inline, virtual] |
Get the descrition of this specific processor.
This function must be overwritten by a subclass.
Implements EMAN::Processor.
Definition at line 4081 of file processor.h.
04082 { 04083 return "Bilateral processing on 2D or 3D volume data. Bilateral processing does non-linear weighted averaging processing within a certain window. "; 04084 }
string EMAN::BilateralProcessor::get_name | ( | ) | const [inline, virtual] |
Get the processor's name.
Each processor is identified by a unique name.
Implements EMAN::Processor.
Definition at line 4076 of file processor.h.
References NAME.
04077 { 04078 return NAME; 04079 }
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.
Reimplemented from EMAN::Processor.
Definition at line 4091 of file processor.h.
References EMAN::EMObject::FLOAT, EMAN::EMObject::INT, and EMAN::TypeDict::put().
04092 { 04093 TypeDict d; 04094 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."); 04095 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."); 04096 d.put("niter", EMObject::INT, "how many times to apply this processing on your data."); 04097 d.put("half_width", EMObject::INT, "processing window size = (2 * half_widthh + 1) ^ 3."); 04098 return d; 04099 }
static Processor* EMAN::BilateralProcessor::NEW | ( | ) | [inline, static] |
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.
image | The image to be processed. |
Implements EMAN::Processor.
Definition at line 3986 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, EMAN::Processor::params, square, and EMAN::EMData::update().
03987 { 03988 if (!image) { 03989 LOGWARN("NULL Image"); 03990 return; 03991 } 03992 03993 float distance_sigma = params["distance_sigma"]; 03994 float value_sigma = params["value_sigma"]; 03995 int max_iter = params["niter"]; 03996 int half_width = params["half_width"]; 03997 03998 if (half_width < distance_sigma) { 03999 LOGWARN("localwidth(=%d) should be larger than distance_sigma=(%f)\n", 04000 half_width, distance_sigma); 04001 } 04002 04003 distance_sigma *= distance_sigma; 04004 04005 float image_sigma = image->get_attr("sigma"); 04006 if (image_sigma > value_sigma) { 04007 LOGWARN("image sigma(=%f) should be smaller than value_sigma=(%f)\n", 04008 image_sigma, value_sigma); 04009 } 04010 value_sigma *= value_sigma; 04011 04012 int nx = image->get_xsize(); 04013 int ny = image->get_ysize(); 04014 int nz = image->get_zsize(); 04015 04016 if(nz==1) { //for 2D image 04017 int width=nx, height=ny; 04018 04019 int i,j,m,n; 04020 04021 float tempfloat1,tempfloat2,tempfloat3; 04022 int index1,index2,index; 04023 int Iter; 04024 int tempint1,tempint3; 04025 04026 tempint1=width; 04027 tempint3=width+2*half_width; 04028 04029 float* mask=(float*)calloc((2*half_width+1)*(2*half_width+1),sizeof(float)); 04030 float* OrgImg=(float*)calloc((2*half_width+width)*(2*half_width+height),sizeof(float)); 04031 float* NewImg=image->get_data(); 04032 04033 for(m=-(half_width);m<=half_width;m++) 04034 for(n=-(half_width);n<=half_width;n++) { 04035 index=(m+half_width)*(2*half_width+1)+(n+half_width); 04036 mask[index]=exp((float)(-(m*m+n*n)/distance_sigma/2.0)); 04037 } 04038 04039 //printf("entering bilateral filtering process \n"); 04040 04041 Iter=0; 04042 while(Iter<max_iter) { 04043 for(i=0;i<height;i++) 04044 for(j=0;j<width;j++) { 04045 index1=(i+half_width)*tempint3+(j+half_width); 04046 index2=i*tempint1+j; 04047 OrgImg[index1]=NewImg[index2]; 04048 } 04049 04050 // Mirror Padding 04051 for(i=0;i<height;i++){ 04052 for(j=0;j<half_width;j++) OrgImg[(i+half_width)*tempint3+(j)]=OrgImg[(i+half_width)*tempint3+(2*half_width-j)]; 04053 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)]; 04054 } 04055 for(i=0;i<half_width;i++){ 04056 for(j=0;j<(width+2*half_width);j++) OrgImg[i*tempint3+j]=OrgImg[(2*half_width-i)*tempint3+j]; 04057 for(j=0;j<(width+2*half_width);j++) OrgImg[(i+height+half_width)*tempint3+j]=OrgImg[(height+half_width-2-i)*tempint3+j]; 04058 } 04059 04060 //printf("finish mirror padding process \n"); 04061 //now mirror padding have been done 04062 04063 for(i=0;i<height;i++){ 04064 //printf("now processing the %d th row \n",i); 04065 for(j=0;j<width;j++){ 04066 tempfloat1=0.0; tempfloat2=0.0; 04067 for(m=-(half_width);m<=half_width;m++) 04068 for(n=-(half_width);n<=half_width;n++){ 04069 index =(m+half_width)*(2*half_width+1)+(n+half_width); 04070 index1=(i+half_width)*tempint3+(j+half_width); 04071 index2=(i+half_width+m)*tempint3+(j+half_width+n); 04072 tempfloat3=(OrgImg[index1]-OrgImg[index2])*(OrgImg[index1]-OrgImg[index2]); 04073 04074 tempfloat3=mask[index]*(1.0f/(1+tempfloat3/value_sigma)); // Lorentz kernel 04075 //tempfloat3=mask[index]*exp(tempfloat3/Sigma2/(-2.0)); // Guassian kernel 04076 tempfloat1+=tempfloat3; 04077 04078 tempfloat2+=tempfloat3*OrgImg[(i+half_width+m)*tempint3+(j+half_width+n)]; 04079 } 04080 NewImg[i*width+j]=tempfloat2/tempfloat1; 04081 } 04082 } 04083 Iter++; 04084 } 04085 04086 //printf("have finished %d th iteration\n ",Iter); 04087 // doneData(); 04088 free(mask); 04089 free(OrgImg); 04090 // end of BilaFilter routine 04091 04092 } 04093 else { //3D case 04094 int width = nx; 04095 int height = ny; 04096 int slicenum = nz; 04097 04098 int slice_size = width * height; 04099 int new_width = width + 2 * half_width; 04100 int new_slice_size = (width + 2 * half_width) * (height + 2 * half_width); 04101 04102 int width1 = 2 * half_width + 1; 04103 int mask_size = width1 * width1; 04104 int old_img_size = (2 * half_width + width) * (2 * half_width + height); 04105 04106 int zstart = -half_width; 04107 int zend = -half_width; 04108 int is_3d = 0; 04109 if (nz > 1) { 04110 mask_size *= width1; 04111 old_img_size *= (2 * half_width + slicenum); 04112 zend = half_width; 04113 is_3d = 1; 04114 } 04115 04116 float *mask = (float *) calloc(mask_size, sizeof(float)); 04117 float *old_img = (float *) calloc(old_img_size, sizeof(float)); 04118 04119 float *new_img = image->get_data(); 04120 04121 for (int p = zstart; p <= zend; p++) { 04122 int cur_p = (p + half_width) * (2 * half_width + 1) * (2 * half_width + 1); 04123 04124 for (int m = -half_width; m <= half_width; m++) { 04125 int cur_m = (m + half_width) * (2 * half_width + 1) + half_width; 04126 04127 for (int n = -half_width; n <= half_width; n++) { 04128 int l = cur_p + cur_m + n; 04129 mask[l] = exp((float) (-(m * m + n * n + p * p * is_3d) / distance_sigma / 2.0f)); 04130 } 04131 } 04132 } 04133 04134 int iter = 0; 04135 while (iter < max_iter) { 04136 for (int k = 0; k < slicenum; k++) { 04137 size_t cur_k1 = (size_t)(k + half_width) * new_slice_size * is_3d; 04138 int cur_k2 = k * slice_size; 04139 04140 for (int i = 0; i < height; i++) { 04141 int cur_i1 = (i + half_width) * new_width; 04142 int cur_i2 = i * width; 04143 04144 for (int j = 0; j < width; j++) { 04145 size_t k1 = cur_k1 + cur_i1 + (j + half_width); 04146 int k2 = cur_k2 + cur_i2 + j; 04147 old_img[k1] = new_img[k2]; 04148 } 04149 } 04150 } 04151 04152 for (int k = 0; k < slicenum; k++) { 04153 size_t cur_k = (k + half_width) * new_slice_size * is_3d; 04154 04155 for (int i = 0; i < height; i++) { 04156 int cur_i = (i + half_width) * new_width; 04157 04158 for (int j = 0; j < half_width; j++) { 04159 size_t k1 = cur_k + cur_i + j; 04160 size_t k2 = cur_k + cur_i + (2 * half_width - j); 04161 old_img[k1] = old_img[k2]; 04162 } 04163 04164 for (int j = 0; j < half_width; j++) { 04165 size_t k1 = cur_k + cur_i + (width + half_width + j); 04166 size_t k2 = cur_k + cur_i + (width + half_width - j - 2); 04167 old_img[k1] = old_img[k2]; 04168 } 04169 } 04170 04171 04172 for (int i = 0; i < half_width; i++) { 04173 int i2 = i * new_width; 04174 int i3 = (2 * half_width - i) * new_width; 04175 for (int j = 0; j < (width + 2 * half_width); j++) { 04176 size_t k1 = cur_k + i2 + j; 04177 size_t k2 = cur_k + i3 + j; 04178 old_img[k1] = old_img[k2]; 04179 } 04180 04181 i2 = (height + half_width + i) * new_width; 04182 i3 = (height + half_width - 2 - i) * new_width; 04183 for (int j = 0; j < (width + 2 * half_width); j++) { 04184 size_t k1 = cur_k + i2 + j; 04185 size_t k2 = cur_k + i3 + j; 04186 old_img[k1] = old_img[k2]; 04187 } 04188 } 04189 } 04190 04191 size_t idx; 04192 for (int k = 0; k < slicenum; k++) { 04193 size_t cur_k = (k + half_width) * new_slice_size; 04194 04195 for (int i = 0; i < height; i++) { 04196 int cur_i = (i + half_width) * new_width; 04197 04198 for (int j = 0; j < width; j++) { 04199 float f1 = 0; 04200 float f2 = 0; 04201 size_t k1 = cur_k + cur_i + (j + half_width); 04202 04203 for (int p = zstart; p <= zend; p++) { 04204 size_t cur_p1 = (p + half_width) * (2 * half_width + 1) * (2 * half_width + 1); 04205 size_t cur_p2 = (k + half_width + p) * new_slice_size; 04206 04207 for (int m = -half_width; m <= half_width; m++) { 04208 size_t cur_m1 = (m + half_width) * (2 * half_width + 1); 04209 size_t cur_m2 = cur_p2 + cur_i + m * new_width + j + half_width; 04210 04211 for (int n = -half_width; n <= half_width; n++) { 04212 size_t k = cur_p1 + cur_m1 + (n + half_width); 04213 size_t k2 = cur_m2 + n; 04214 float f3 = Util::square(old_img[k1] - old_img[k2]); 04215 04216 f3 = mask[k] * (1.0f / (1 + f3 / value_sigma)); 04217 f1 += f3; 04218 size_t l1 = cur_m2 + n; 04219 f2 += f3 * old_img[l1]; 04220 } 04221 04222 idx = (size_t)k * height * width + i * width + j; 04223 new_img[idx] = f2 / f1; 04224 } 04225 } 04226 } 04227 } 04228 } 04229 iter++; 04230 } 04231 if( mask ) { 04232 free(mask); 04233 mask = 0; 04234 } 04235 04236 if( old_img ) { 04237 free(old_img); 04238 old_img = 0; 04239 } 04240 } 04241 04242 image->update(); 04243 }
const string BilateralProcessor::NAME = "filter.bilateral" [static] |