#include <processor.h>
Inheritance diagram for EMAN::GradientPlaneRemoverProcessor:
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.gradientPlaneRemover" |
mask[in] | optional EMData object to mask the pixels used to fit the plane | |
changeZero[in] | optional bool to specify if the zero value pixels are modified | |
planeParam[out] | optional return parameters [nx, ny, nz, cx, cy, cz] for the fitted plane defined as (x-cx)*nx+(y-cy)*ny+(z-cz)*nz=0 |
Definition at line 3617 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 3631 of file processor.h. 03632 { 03633 return "Remove gradient by least square plane fit"; 03634 }
|
|
Get the processor's name. Each processor is identified by a unique name.
Implements EMAN::Processor. Definition at line 3622 of file processor.h. Referenced by process_inplace(). 03623 {
03624 return NAME;
03625 }
|
|
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 3636 of file processor.h. References EMAN::TypeDict::put(). 03637 { 03638 TypeDict d; 03639 d.put("mask", EMObject::EMDATA, "mask object: nonzero pixel positions will be used to fit plane. default = 0"); 03640 d.put("changeZero", EMObject::INT, "if zero pixels are modified when removing gradient. default = 0"); 03641 d.put("planeParam", EMObject::FLOATARRAY, "fitted plane parameters output"); 03642 return d; 03643 }
|
|
Definition at line 3626 of file processor.h. 03627 { 03628 return new GradientPlaneRemoverProcessor(); 03629 }
|
|
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 2768 of file processor.cpp. References dm, EMAN::EMData::get_data(), get_name(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::Dict::has_key(), ImageDimensionException, LOGERR, LOGWARN, nx, ny, EMAN::EMData::update(), and V. 02769 { 02770 if (!image) { 02771 LOGWARN("NULL Image"); 02772 return; 02773 } 02774 02775 int nz = image->get_zsize(); 02776 if (nz > 1) { 02777 LOGERR("%s Processor doesn't support 3D model", get_name().c_str()); 02778 throw ImageDimensionException("3D map not supported"); 02779 } 02780 02781 int nx = image->get_xsize(); 02782 int ny = image->get_ysize(); 02783 float *d = image->get_data(); 02784 EMData *mask = 0; 02785 float *dm = 0; 02786 if (params.has_key("mask")) { 02787 mask = params["mask"]; 02788 if (nx!=mask->get_xsize() || ny!=mask->get_ysize()) { 02789 LOGERR("%s Processor requires same size mask image", get_name().c_str()); 02790 throw ImageDimensionException("wrong size mask image"); 02791 } 02792 dm = mask->get_data(); 02793 } 02794 int count = 0; 02795 if (dm) { 02796 for(int i=0; i<nx*ny; i++) { 02797 if(dm[i]) count++; 02798 } 02799 } 02800 else { 02801 count = nx * ny; 02802 } 02803 if(count<3) { 02804 LOGERR("%s Processor requires at least 3 pixels to fit a plane", get_name().c_str()); 02805 throw ImageDimensionException("too few usable pixels to fit a plane"); 02806 } 02807 // Allocate the working space 02808 gsl_vector *S=gsl_vector_calloc(3); 02809 gsl_matrix *A=gsl_matrix_calloc(count,3); 02810 gsl_matrix *V=gsl_matrix_calloc(3,3); 02811 02812 double m[3] = {0, 0, 0}; 02813 int index=0; 02814 if (dm) { 02815 for(int j=0; j<ny; j++){ 02816 for(int i=0; i<nx; i++){ 02817 int ij=j*nx+i; 02818 if(dm[ij]) { 02819 m[0]+=i; // x 02820 m[1]+=j; // y 02821 m[2]+=d[ij]; // z 02822 /*printf("index=%d/%d\ti,j=%d,%d\tval=%g\txm,ym,zm=%g,%g,%g\n", \ 02823 index,count,i,j,d[ij],m[0]/(index+1),m[1]/(index+1),m[2]/(index+1));*/ 02824 index++; 02825 } 02826 } 02827 } 02828 } 02829 else { 02830 for(int j=0; j<ny; j++){ 02831 for(int i=0; i<nx; i++){ 02832 int ij=j*nx+i; 02833 m[0]+=i; // x 02834 m[1]+=j; // y 02835 m[2]+=d[ij]; // z 02836 /*printf("index=%d/%d\ti,j=%d,%d\tval=%g\txm,ym,zm=%g,%g,%g\n", \ 02837 index,count,i,j,d[ij],m[0]/(index+1),m[1]/(index+1),m[2]/(index+1));*/ 02838 index++; 02839 } 02840 } 02841 } 02842 02843 for(int i=0; i<3; i++) m[i]/=count; // compute center of the plane 02844 02845 index=0; 02846 if (dm) { 02847 for(int j=0; j<ny; j++){ 02848 for(int i=0; i<nx; i++){ 02849 int ij=j*nx+i; 02850 if(dm[ij]) { 02851 //printf("index=%d/%d\ti,j=%d,%d\tval=%g\n",index,count,i,j,d[index]); 02852 gsl_matrix_set(A,index,0,i-m[0]); 02853 gsl_matrix_set(A,index,1,j-m[1]); 02854 gsl_matrix_set(A,index,2,d[ij]-m[2]); 02855 index++; 02856 } 02857 } 02858 } 02859 mask->update(); 02860 } 02861 else { 02862 for(int j=0; j<ny; j++){ 02863 for(int i=0; i<nx; i++){ 02864 int ij=j*nx+i; 02865 //printf("index=%d/%d\ti,j=%d,%d\tval=%g\n",index,count,i,j,d[index]); 02866 gsl_matrix_set(A,index,0,i-m[0]); 02867 gsl_matrix_set(A,index,1,j-m[1]); 02868 gsl_matrix_set(A,index,2,d[ij]-m[2]); 02869 index++; 02870 } 02871 } 02872 } 02873 02874 // SVD decomposition and use the V vector associated with smallest singular value as the plan normal 02875 gsl_linalg_SV_decomp_jacobi(A, V, S); 02876 02877 double n[3]; 02878 for(int i=0; i<3; i++) n[i] = gsl_matrix_get(V, i, 2); 02879 02880 #ifdef DEBUG 02881 printf("S=%g,%g,%g\n",gsl_vector_get(S,0), gsl_vector_get(S,1), gsl_vector_get(S,2)); 02882 printf("V[0,:]=%g,%g,%g\n",gsl_matrix_get(V,0,0), gsl_matrix_get(V,0,1),gsl_matrix_get(V,0,2)); 02883 printf("V[1,:]=%g,%g,%g\n",gsl_matrix_get(V,1,0), gsl_matrix_get(V,1,1),gsl_matrix_get(V,1,2)); 02884 printf("V[2,:]=%g,%g,%g\n",gsl_matrix_get(V,2,0), gsl_matrix_get(V,2,1),gsl_matrix_get(V,2,2)); 02885 printf("Fitted plane: p0=%g,%g,%g\tn=%g,%g,%g\n",m[0],m[1],m[2],n[0],n[1],n[2]); 02886 #endif 02887 02888 int changeZero = 0; 02889 if (params.has_key("changeZero")) changeZero = params["changeZero"]; 02890 if (changeZero) { 02891 for(int j=0; j<nx; j++){ 02892 for(int i=0; i<ny; i++){ 02893 int ij = j*nx+i; 02894 d[ij]-=static_cast<float>(-((i-m[0])*n[0]+(j-m[1])*n[1])/n[2]+m[2]); 02895 } 02896 } 02897 } 02898 else { 02899 for(int j=0; j<nx; j++){ 02900 for(int i=0; i<ny; i++){ 02901 int ij = j*nx+i; 02902 if(d[ij]) d[ij]-=static_cast<float>(-((i-m[0])*n[0]+(j-m[1])*n[1])/n[2]+m[2]); 02903 } 02904 } 02905 } 02906 image->update(); 02907 // set return plane parameters 02908 vector< float > planeParam; 02909 planeParam.resize(6); 02910 for(int i=0; i<3; i++) planeParam[i] = static_cast<float>(n[i]); 02911 for(int i=0; i<3; i++) planeParam[i+3] = static_cast<float>(m[i]); 02912 params["planeParam"]=EMObject(planeParam); 02913 }
|
|
Definition at line 135 of file processor.cpp. |