#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 3614 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 3628 of file processor.h. 03629 { 03630 return "Remove gradient by least square plane fit"; 03631 }
|
|
Get the processor's name. Each processor is identified by a unique name.
Implements EMAN::Processor. Definition at line 3619 of file processor.h. Referenced by process_inplace(). 03620 {
03621 return NAME;
03622 }
|
|
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 3633 of file processor.h. References EMAN::TypeDict::put(). 03634 { 03635 TypeDict d; 03636 d.put("mask", EMObject::EMDATA, "mask object: nonzero pixel positions will be used to fit plane. default = 0"); 03637 d.put("changeZero", EMObject::INT, "if zero pixels are modified when removing gradient. default = 0"); 03638 d.put("planeParam", EMObject::FLOATARRAY, "fitted plane parameters output"); 03639 return d; 03640 }
|
|
Definition at line 3623 of file processor.h. 03624 { 03625 return new GradientPlaneRemoverProcessor(); 03626 }
|
|
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 2693 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. 02694 { 02695 if (!image) { 02696 LOGWARN("NULL Image"); 02697 return; 02698 } 02699 02700 int nz = image->get_zsize(); 02701 if (nz > 1) { 02702 LOGERR("%s Processor doesn't support 3D model", get_name().c_str()); 02703 throw ImageDimensionException("3D map not supported"); 02704 } 02705 02706 int nx = image->get_xsize(); 02707 int ny = image->get_ysize(); 02708 float *d = image->get_data(); 02709 EMData *mask = 0; 02710 float *dm = 0; 02711 if (params.has_key("mask")) { 02712 mask = params["mask"]; 02713 if (nx!=mask->get_xsize() || ny!=mask->get_ysize()) { 02714 LOGERR("%s Processor requires same size mask image", get_name().c_str()); 02715 throw ImageDimensionException("wrong size mask image"); 02716 } 02717 dm = mask->get_data(); 02718 } 02719 int count = 0; 02720 if (dm) { 02721 for(int i=0; i<nx*ny; i++) { 02722 if(dm[i]) count++; 02723 } 02724 } 02725 else { 02726 count = nx * ny; 02727 } 02728 if(count<3) { 02729 LOGERR("%s Processor requires at least 3 pixels to fit a plane", get_name().c_str()); 02730 throw ImageDimensionException("too few usable pixels to fit a plane"); 02731 } 02732 // Allocate the working space 02733 gsl_vector *S=gsl_vector_calloc(3); 02734 gsl_matrix *A=gsl_matrix_calloc(count,3); 02735 gsl_matrix *V=gsl_matrix_calloc(3,3); 02736 02737 double m[3] = {0, 0, 0}; 02738 int index=0; 02739 if (dm) { 02740 for(int j=0; j<ny; j++){ 02741 for(int i=0; i<nx; i++){ 02742 int ij=j*nx+i; 02743 if(dm[ij]) { 02744 m[0]+=i; // x 02745 m[1]+=j; // y 02746 m[2]+=d[ij]; // z 02747 /*printf("index=%d/%d\ti,j=%d,%d\tval=%g\txm,ym,zm=%g,%g,%g\n", \ 02748 index,count,i,j,d[ij],m[0]/(index+1),m[1]/(index+1),m[2]/(index+1));*/ 02749 index++; 02750 } 02751 } 02752 } 02753 } 02754 else { 02755 for(int j=0; j<ny; j++){ 02756 for(int i=0; i<nx; i++){ 02757 int ij=j*nx+i; 02758 m[0]+=i; // x 02759 m[1]+=j; // y 02760 m[2]+=d[ij]; // z 02761 /*printf("index=%d/%d\ti,j=%d,%d\tval=%g\txm,ym,zm=%g,%g,%g\n", \ 02762 index,count,i,j,d[ij],m[0]/(index+1),m[1]/(index+1),m[2]/(index+1));*/ 02763 index++; 02764 } 02765 } 02766 } 02767 02768 for(int i=0; i<3; i++) m[i]/=count; // compute center of the plane 02769 02770 index=0; 02771 if (dm) { 02772 for(int j=0; j<ny; j++){ 02773 for(int i=0; i<nx; i++){ 02774 int ij=j*nx+i; 02775 if(dm[ij]) { 02776 //printf("index=%d/%d\ti,j=%d,%d\tval=%g\n",index,count,i,j,d[index]); 02777 gsl_matrix_set(A,index,0,i-m[0]); 02778 gsl_matrix_set(A,index,1,j-m[1]); 02779 gsl_matrix_set(A,index,2,d[ij]-m[2]); 02780 index++; 02781 } 02782 } 02783 } 02784 mask->update(); 02785 } 02786 else { 02787 for(int j=0; j<ny; j++){ 02788 for(int i=0; i<nx; i++){ 02789 int ij=j*nx+i; 02790 //printf("index=%d/%d\ti,j=%d,%d\tval=%g\n",index,count,i,j,d[index]); 02791 gsl_matrix_set(A,index,0,i-m[0]); 02792 gsl_matrix_set(A,index,1,j-m[1]); 02793 gsl_matrix_set(A,index,2,d[ij]-m[2]); 02794 index++; 02795 } 02796 } 02797 } 02798 02799 // SVD decomposition and use the V vector associated with smallest singular value as the plan normal 02800 gsl_linalg_SV_decomp_jacobi(A, V, S); 02801 02802 double n[3]; 02803 for(int i=0; i<3; i++) n[i] = gsl_matrix_get(V, i, 2); 02804 02805 #ifdef DEBUG 02806 printf("S=%g,%g,%g\n",gsl_vector_get(S,0), gsl_vector_get(S,1), gsl_vector_get(S,2)); 02807 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)); 02808 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)); 02809 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)); 02810 printf("Fitted plane: p0=%g,%g,%g\tn=%g,%g,%g\n",m[0],m[1],m[2],n[0],n[1],n[2]); 02811 #endif 02812 02813 int changeZero = 0; 02814 if (params.has_key("changeZero")) changeZero = params["changeZero"]; 02815 if (changeZero) { 02816 for(int j=0; j<nx; j++){ 02817 for(int i=0; i<ny; i++){ 02818 int ij = j*nx+i; 02819 d[ij]-=static_cast<float>(-((i-m[0])*n[0]+(j-m[1])*n[1])/n[2]+m[2]); 02820 } 02821 } 02822 } 02823 else { 02824 for(int j=0; j<nx; j++){ 02825 for(int i=0; i<ny; i++){ 02826 int ij = j*nx+i; 02827 if(d[ij]) d[ij]-=static_cast<float>(-((i-m[0])*n[0]+(j-m[1])*n[1])/n[2]+m[2]); 02828 } 02829 } 02830 } 02831 image->update(); 02832 // set return plane parameters 02833 vector< float > planeParam; 02834 planeParam.resize(6); 02835 for(int i=0; i<3; i++) planeParam[i] = static_cast<float>(n[i]); 02836 for(int i=0; i<3; i++) planeParam[i+3] = static_cast<float>(m[i]); 02837 params["planeParam"]=EMObject(planeParam); 02838 }
|
|
Definition at line 129 of file processor.cpp. |