Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

EMAN::GradientPlaneRemoverProcessor Class Reference

Gradient removed by least square plane fit. More...

#include <processor.h>

Inheritance diagram for EMAN::GradientPlaneRemoverProcessor:

Inheritance graph
[legend]
Collaboration diagram for EMAN::GradientPlaneRemoverProcessor:

Collaboration graph
[legend]
List of all members.

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

ProcessorNEW ()

Static Public Attributes

const string NAME = "filter.gradientPlaneRemover"

Detailed Description

Gradient removed by least square plane fit.

Parameters:
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
Author:
Wen Jiang
Date:
2006/7/18

Definition at line 3617 of file processor.h.


Member Function Documentation

string EMAN::GradientPlaneRemoverProcessor::get_desc  )  const [inline, virtual]
 

Get the descrition of this specific processor.

This function must be overwritten by a subclass.

Returns:
The description of this processor.

Implements EMAN::Processor.

Definition at line 3631 of file processor.h.

03632                 {
03633                         return "Remove gradient by least square plane fit";
03634                 }

string EMAN::GradientPlaneRemoverProcessor::get_name  )  const [inline, virtual]
 

Get the processor's name.

Each processor is identified by a unique name.

Returns:
The processor's name.

Implements EMAN::Processor.

Definition at line 3622 of file processor.h.

Referenced by process_inplace().

03623                 {
03624                         return NAME;
03625                 }

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.

Returns:
A dictionary containing the parameter info.

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                 }

Processor* EMAN::GradientPlaneRemoverProcessor::NEW  )  [inline, static]
 

Definition at line 3626 of file processor.h.

03627                 {
03628                         return new GradientPlaneRemoverProcessor();
03629                 }

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.

Parameters:
image The image to be processed.

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 }


Member Data Documentation

const string GradientPlaneRemoverProcessor::NAME = "filter.gradientPlaneRemover" [static]
 

Definition at line 135 of file processor.cpp.


The documentation for this class was generated from the following files:
Generated on Tue Jul 12 13:51:09 2011 for EMAN2 by  doxygen 1.3.9.1