#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 | |
static Processor * | NEW () |
Static Public Attributes | |
static 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 3660 of file processor.h.
string EMAN::GradientPlaneRemoverProcessor::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 3674 of file processor.h.
string EMAN::GradientPlaneRemoverProcessor::get_name | ( | ) | const [inline, virtual] |
Get the processor's name.
Each processor is identified by a unique name.
Implements EMAN::Processor.
Definition at line 3665 of file processor.h.
References NAME.
Referenced by process_inplace().
03666 { 03667 return NAME; 03668 }
TypeDict EMAN::GradientPlaneRemoverProcessor::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 3679 of file processor.h.
References EMAN::EMObject::EMDATA, EMAN::EMObject::FLOATARRAY, EMAN::EMObject::INT, and EMAN::TypeDict::put().
03680 { 03681 TypeDict d; 03682 d.put("mask", EMObject::EMDATA, "mask object: nonzero pixel positions will be used to fit plane. default = 0"); 03683 d.put("changeZero", EMObject::INT, "if zero pixels are modified when removing gradient. default = 0"); 03684 d.put("planeParam", EMObject::FLOATARRAY, "fitted plane parameters output"); 03685 return d; 03686 }
static Processor* EMAN::GradientPlaneRemoverProcessor::NEW | ( | ) | [inline, static] |
Definition at line 3669 of file processor.h.
03670 { 03671 return new GradientPlaneRemoverProcessor(); 03672 }
void GradientPlaneRemoverProcessor::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 2865 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, EMAN::Processor::params, EMAN::EMData::update(), and V.
02866 { 02867 if (!image) { 02868 LOGWARN("NULL Image"); 02869 return; 02870 } 02871 02872 int nz = image->get_zsize(); 02873 if (nz > 1) { 02874 LOGERR("%s Processor doesn't support 3D model", get_name().c_str()); 02875 throw ImageDimensionException("3D map not supported"); 02876 } 02877 02878 int nx = image->get_xsize(); 02879 int ny = image->get_ysize(); 02880 float *d = image->get_data(); 02881 EMData *mask = 0; 02882 float *dm = 0; 02883 if (params.has_key("mask")) { 02884 mask = params["mask"]; 02885 if (nx!=mask->get_xsize() || ny!=mask->get_ysize()) { 02886 LOGERR("%s Processor requires same size mask image", get_name().c_str()); 02887 throw ImageDimensionException("wrong size mask image"); 02888 } 02889 dm = mask->get_data(); 02890 } 02891 int count = 0; 02892 if (dm) { 02893 for(int i=0; i<nx*ny; i++) { 02894 if(dm[i]) count++; 02895 } 02896 } 02897 else { 02898 count = nx * ny; 02899 } 02900 if(count<3) { 02901 LOGERR("%s Processor requires at least 3 pixels to fit a plane", get_name().c_str()); 02902 throw ImageDimensionException("too few usable pixels to fit a plane"); 02903 } 02904 // Allocate the working space 02905 gsl_vector *S=gsl_vector_calloc(3); 02906 gsl_matrix *A=gsl_matrix_calloc(count,3); 02907 gsl_matrix *V=gsl_matrix_calloc(3,3); 02908 02909 double m[3] = {0, 0, 0}; 02910 int index=0; 02911 if (dm) { 02912 for(int j=0; j<ny; j++){ 02913 for(int i=0; i<nx; i++){ 02914 int ij=j*nx+i; 02915 if(dm[ij]) { 02916 m[0]+=i; // x 02917 m[1]+=j; // y 02918 m[2]+=d[ij]; // z 02919 /*printf("index=%d/%d\ti,j=%d,%d\tval=%g\txm,ym,zm=%g,%g,%g\n", \ 02920 index,count,i,j,d[ij],m[0]/(index+1),m[1]/(index+1),m[2]/(index+1));*/ 02921 index++; 02922 } 02923 } 02924 } 02925 } 02926 else { 02927 for(int j=0; j<ny; j++){ 02928 for(int i=0; i<nx; i++){ 02929 int ij=j*nx+i; 02930 m[0]+=i; // x 02931 m[1]+=j; // y 02932 m[2]+=d[ij]; // z 02933 /*printf("index=%d/%d\ti,j=%d,%d\tval=%g\txm,ym,zm=%g,%g,%g\n", \ 02934 index,count,i,j,d[ij],m[0]/(index+1),m[1]/(index+1),m[2]/(index+1));*/ 02935 index++; 02936 } 02937 } 02938 } 02939 02940 for(int i=0; i<3; i++) m[i]/=count; // compute center of the plane 02941 02942 index=0; 02943 if (dm) { 02944 for(int j=0; j<ny; j++){ 02945 for(int i=0; i<nx; i++){ 02946 int ij=j*nx+i; 02947 if(dm[ij]) { 02948 //printf("index=%d/%d\ti,j=%d,%d\tval=%g\n",index,count,i,j,d[index]); 02949 gsl_matrix_set(A,index,0,i-m[0]); 02950 gsl_matrix_set(A,index,1,j-m[1]); 02951 gsl_matrix_set(A,index,2,d[ij]-m[2]); 02952 index++; 02953 } 02954 } 02955 } 02956 mask->update(); 02957 } 02958 else { 02959 for(int j=0; j<ny; j++){ 02960 for(int i=0; i<nx; i++){ 02961 int ij=j*nx+i; 02962 //printf("index=%d/%d\ti,j=%d,%d\tval=%g\n",index,count,i,j,d[index]); 02963 gsl_matrix_set(A,index,0,i-m[0]); 02964 gsl_matrix_set(A,index,1,j-m[1]); 02965 gsl_matrix_set(A,index,2,d[ij]-m[2]); 02966 index++; 02967 } 02968 } 02969 } 02970 02971 // SVD decomposition and use the V vector associated with smallest singular value as the plan normal 02972 gsl_linalg_SV_decomp_jacobi(A, V, S); 02973 02974 double n[3]; 02975 for(int i=0; i<3; i++) n[i] = gsl_matrix_get(V, i, 2); 02976 02977 #ifdef DEBUG 02978 printf("S=%g,%g,%g\n",gsl_vector_get(S,0), gsl_vector_get(S,1), gsl_vector_get(S,2)); 02979 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)); 02980 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)); 02981 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)); 02982 printf("Fitted plane: p0=%g,%g,%g\tn=%g,%g,%g\n",m[0],m[1],m[2],n[0],n[1],n[2]); 02983 #endif 02984 02985 int changeZero = 0; 02986 if (params.has_key("changeZero")) changeZero = params["changeZero"]; 02987 if (changeZero) { 02988 for(int j=0; j<nx; j++){ 02989 for(int i=0; i<ny; i++){ 02990 int ij = j*nx+i; 02991 d[ij]-=static_cast<float>(-((i-m[0])*n[0]+(j-m[1])*n[1])/n[2]+m[2]); 02992 } 02993 } 02994 } 02995 else { 02996 for(int j=0; j<nx; j++){ 02997 for(int i=0; i<ny; i++){ 02998 int ij = j*nx+i; 02999 if(d[ij]) d[ij]-=static_cast<float>(-((i-m[0])*n[0]+(j-m[1])*n[1])/n[2]+m[2]); 03000 } 03001 } 03002 } 03003 image->update(); 03004 // set return plane parameters 03005 vector< float > planeParam; 03006 planeParam.resize(6); 03007 for(int i=0; i<3; i++) planeParam[i] = static_cast<float>(n[i]); 03008 for(int i=0; i<3; i++) planeParam[i+3] = static_cast<float>(m[i]); 03009 params["planeParam"]=EMObject(planeParam); 03010 }
const string GradientPlaneRemoverProcessor::NAME = "filter.gradientPlaneRemover" [static] |