ProgVolumeEnhanceContrast Class Reference

#include <volume_enhance_contrast.h>

Inheritance diagram for ProgVolumeEnhanceContrast:

Collaboration diagram for ProgVolumeEnhanceContrast:

Public Member Functions
void	show ()
void	enhance (MultidimArray< double > &vol)
Public Member Functions inherited from XmippProgram
const char *	getParam (const char *param, int arg=0)
const char *	getParam (const char *param, const char *subparam, int arg=0)
int	getIntParam (const char *param, int arg=0)
int	getIntParam (const char *param, const char *subparam, int arg=0)
double	getDoubleParam (const char *param, int arg=0)
double	getDoubleParam (const char *param, const char *subparam, int arg=0)
float	getFloatParam (const char *param, int arg=0)
float	getFloatParam (const char *param, const char *subparam, int arg=0)
void	getListParam (const char *param, StringVector &list)
int	getCountParam (const char *param)
bool	checkParam (const char *param)
bool	existsParam (const char *param)
void	addParamsLine (const String &line)
void	addParamsLine (const char *line)
ParamDef *	getParamDef (const char *param) const
virtual void	quit (int exit_code=0) const
virtual int	tryRun ()
void	initProgress (size_t total, size_t stepBin=60)
void	setProgress (size_t value=0)
void	endProgress ()
void	processDefaultComment (const char *param, const char *left)
void	setDefaultComment (const char *param, const char *comment)
virtual void	initComments ()
void	setProgramName (const char *name)
void	addUsageLine (const char *line, bool verbatim=false)
void	clearUsage ()
void	addExampleLine (const char *example, bool verbatim=true)
void	addSeeAlsoLine (const char *seeAlso)
void	addKeywords (const char *keywords)
const char *	name () const
virtual void	usage (int verb=0) const
virtual void	usage (const String &param, int verb=2)
int	version () const
virtual void	show () const
virtual void	read (int argc, const char **argv, bool reportErrors=true)
virtual void	read (int argc, char **argv, bool reportErrors=true)
void	read (const String &argumentsLine)
	XmippProgram ()
	XmippProgram (int argc, const char **argv)
virtual	~XmippProgram ()

Public Attributes
FileName	fnIn
FileName	fnOut
double	alpha
double	lowerIntensity
bool	removeBg
FileName	fnMask
Public Attributes inherited from XmippProgram
bool	doRun
bool	runWithoutArgs
int	verbose
	Verbosity level. More...
int	debug

Protected Member Functions
void	defineParams ()
void	readParams ()
void	run ()
Protected Member Functions inherited from XmippProgram
void	defineCommons ()

Additional Inherited Members
Protected Attributes inherited from XmippProgram
int	errorCode
ProgramDef *	progDef
	Program definition and arguments parser. More...
std::map< String, CommentList >	defaultComments
int	argc
	Original command line arguments. More...
const char **	argv

Detailed Description

Parameters for enhance contrast program

Definition at line 47 of file volume_enhance_contrast.h.

Member Function Documentation

defineParams()

void ProgVolumeEnhanceContrast::defineParams ( )

protectedvirtual

Params definition

Reimplemented from XmippProgram.

Definition at line 32 of file volume_enhance_contrast.cpp.

{
    addUsageLine("Enhance volume contrast by applying a non-linear transformation to gray levels");
    addUsageLine("+The contrast enhancement is based on an evolution of an article of");
    addUsageLine("+[[http://www.ncbi.nlm.nih.gov/pubmed/16579377][Sean Matz]] and is further ");
    addUsageLine("+explained [[http://biocomp.cnb.csic.es/~coss/Articulos/Fuentes2010.pdf][here]]");
    addUsageLine("+The algorithm detects the molecule borders and applies a nonlinear transformation ");
    addUsageLine("+of the gray values");
    addSeeAlsoLine("volume_correct_bfactor");
    addParamsLine("   -i <file>                  : Input volume");
    addParamsLine("  [-o <file=\"\">]                 : Output volume");
    addParamsLine("  [--removeBackground]        : Remove the noise of the background");
    addParamsLine("  [--alpha+ <a=0.01>]         : Confidence interval for background identification");
    addParamsLine("  [--lowerIntensity+ <intensity=0>]     : Only process if the gray value is higher than this value");
    addParamsLine("  [--saveMask+ <filename=\"\">]  : Filename for the background mask");
    addExampleLine("xmipp_volume_enhance_contrast -i volume.vol -o volumeEnhanced.vol");
}

enhance()

void ProgVolumeEnhanceContrast::enhance

(

MultidimArray< double > &

vol

)

Enhance contrast of a volume.

Definition at line 86 of file volume_enhance_contrast.cpp.

{
    // 1.-Scale volume between 0 and 255----------------------------------
    double minVal=0., maxVal=0.;
    vol.computeDoubleMinMax(minVal,maxVal);
    vol.rangeAdjust(0,255);

    double vol_avg = vol.computeAvg();

    // Padd volume with two new rows and cols
    if( FINISHINGZ(vol) - STARTINGZ(vol) > 0 )
    {
        vol.selfWindow( STARTINGZ(vol)-2,STARTINGY(vol)-2,STARTINGX(vol)-2,
                    FINISHINGZ(vol)+2, FINISHINGY(vol)+2,FINISHINGX(vol)+2, vol_avg );
    }
    else
    {
        vol.selfWindow( STARTINGZ(vol),STARTINGY(vol)-2,STARTINGX(vol)-2,
                    FINISHINGZ(vol), FINISHINGY(vol)+2,FINISHINGX(vol)+2, vol_avg );
    }

    Image<double> save;
#ifdef DEBUG

    save()=vol;
    save.write("PPP1_init.vol");
#endif

    // 2.-Elimination of the background-----------------------------------
    MultidimArray<double> mask;
    double bg_mean;
    detectBackground(vol,mask,0.01,bg_mean);
#ifdef DEBUG

    save()=mask;
    save.write("PPPmask.vol");
#endif

    // We change 0<->1
    FOR_ALL_ELEMENTS_IN_ARRAY3D(mask)
    {
        if (A3D_ELEM(mask,k,i,j)==1)
        {
            A3D_ELEM(mask,k,i,j)=0;
        }
        else
        {
            A3D_ELEM(mask,k,i,j)=1;
        }
    } // Now if 0:background and if 1:mol
#ifdef DEBUG
    save()=mask;
    save.write("PPPmask_c.vol");
#endif

    if (fnMask!="")
    {
        save()=mask;
        save.write(fnMask);
    }
    if (removeBg==true)
    {
        // We put all the background with the same value
        FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(vol)
        if (DIRECT_A3D_ELEM(mask,k,i,j)==0)
            DIRECT_A3D_ELEM(vol,k,i,j)=bg_mean;
    }
#ifdef DEBUG
    save()=vol;
    save.write("PPP2_no_bg.vol");
#endif

    // 3.-BSplines + diff = edges-------------------------------------------
    MultidimArray<double> vol_edge;
    computeEdges(vol,vol_edge);

#ifdef DEBUG
    save()=vol_edge;
    save.write("PPP3_edge.vol");
#endif
    // 4.-MEAN EDGE GRAY VALUE

    // 4.1.-Variable Neighbourhood
#define COMPUTE_STATISTICS(V,N,k,i,j,avg,stddev,cubeSize) \
     {\
   int k0,kF,i0,iF,j0,jF; \
     k0=XMIPP_MAX(k-N,STARTINGZ(V)); \
     kF=XMIPP_MIN(k+N,FINISHINGZ(V)); \
     i0=XMIPP_MAX(i-N,STARTINGY(V)); \
     iF=XMIPP_MIN(i+N,FINISHINGY(V)); \
     j0=XMIPP_MAX(j-N,STARTINGX(V)); \
         jF=XMIPP_MIN(j+N,FINISHINGX(V)); \
         double sum=0, sum2=0; \
         for (int kk=k0; kk<=kF; ++kk) \
    for (int ii=i0; ii<=iF; ++ii) \
     for (int jj=j0; jj<=jF; ++jj) \
     { \
      double v=A3D_ELEM(V,kk,ii,jj); \
      sum+=v; \
      sum2+=v*v; \
     } \
   cubeSize=(kF-k0+1)*(iF-i0+1)*(jF-j0+1); \
   avg=sum/cubeSize; \
   stddev=sum2/cubeSize-avg*avg; \
   stddev=sqrt(XMIPP_MAX(stddev,0)); \
  }

    // We use gaussian because we have a T-Student with more than 100
    float z=icdf_gauss(1-(0.01/2));  // degrees of freedom
    // We create a volume with the neighbor sizes
    MultidimArray<int> vol_tam;
    vol_tam.resize(vol);
    vol_tam.initConstant(4);

    FOR_ALL_ELEMENTS_IN_ARRAY3D(vol)
    {
        if (A3D_ELEM(mask,k,i,j)==1)
        {
            // only compute if the pixel is not background
            int cubeDim=1; // 3x3x3
            int newCubeDim=1; // 3x3x3
            int max_tam=4; // 9x9x9
            double avgSmall, stddevSmall, NSmall;
            COMPUTE_STATISTICS(vol,cubeDim,k,i,j,avgSmall,stddevSmall,NSmall);
            bool found=false;
            while (!found)
            {
                // Big Cube
                newCubeDim=newCubeDim+1;
                double avgBig, stddevBig, NBig;
                COMPUTE_STATISTICS(vol,newCubeDim,k,i,j,avgBig,stddevBig,NBig);
                // Computing the degrees of freedom
                //double num, den;
                //num=(((stddevSmall*stddevSmall)/NSmall)+((stddevBig*stddevBig)/NBig))*
                // (((stddevSmall*stddevSmall)/NSmall)+((stddevBig*stddevBig)/NBig));
                //den=((stddevSmall*stddevSmall)/((NSmall-1)*(NSmall*NSmall)))+
                // ((stddevBig*stddevBig)/((NBig-1)*(NBig*NBig)));
                //double df;
                //df=num/den; // Degrees of freedom
                // Confidence Intervals --> Comparison of two cubes
                double K, d, A, B;
                K=sqrt(((stddevSmall*stddevSmall)/NSmall)+((stddevBig*stddevBig)/NBig));
                d=avgSmall-avgBig;
                A=d-(K*z);
                B=d+(K*z);
                // If the interval [A,B] contains the cero there are the same
                if ((A<0) && (0<B))
                {  // equal continue or not
                    if (newCubeDim>=max_tam)
                    { // not continue
                        found=true;
                        A3D_ELEM(vol_tam,k,i,j)=newCubeDim;
                    }
                    else
                    { // continue searching
                        avgSmall=avgBig;
                        stddevSmall=stddevBig;
                        NSmall=NBig;
                    }
                }
                else // Not equal -> we stop searching
                {
                    found=true;  // Size found
                    A3D_ELEM(vol_tam,k,i,j)=newCubeDim-1; // it's the latest
                }

            } // end of while
        } // end of if
    } // end of FOR_ALL_ELEMENTS
#ifdef DEBUG
    save()=vol_tam;
    save.write("PPP4_vol_tam.vol");
#endif

    // 4.2.- Compute the mean edge gray value
    MultidimArray<double> VolxEdge;
    VolxEdge=vol;
    VolxEdge*=vol_edge;
    // Macro for compute de Mean Edge Gray Value
#define COMPUTE_MEGV(VxE,V_E,N,k,i,j,pixel) \
     {\
   int k0,kF,i0,iF,j0,jF; \
     k0=XMIPP_MAX(k-N,STARTINGZ(V_E)); \
     kF=XMIPP_MIN(k+N,FINISHINGZ(V_E)); \
     i0=XMIPP_MAX(i-N,STARTINGY(V_E)); \
     iF=XMIPP_MIN(i+N,FINISHINGY(V_E)); \
     j0=XMIPP_MAX(j-N,STARTINGX(V_E)); \
         jF=XMIPP_MIN(j+N,FINISHINGX(V_E)); \
         double sum1=0, sum2=0; \
         for (int kk=k0; kk<=kF; kk++) \
    for (int ii=i0; ii<=iF; ii++) \
     for (int jj=j0; jj<=jF; jj++) \
     { \
      double v1=A3D_ELEM(VxE,kk,ii,jj); \
      sum1=sum1+v1; \
      double v2=A3D_ELEM(V_E,kk,ii,jj); \
      sum2=sum2+v2; \
     } \
   pixel=sum1/sum2; \
  }

    MultidimArray<double> Vol_E;
    Vol_E=vol;
    int tam;
    double pixel;
    // We compute the MEGV for all the volume
    FOR_ALL_ELEMENTS_IN_ARRAY3D(vol)
    { // Only compute if no background
        if (A3D_ELEM(mask,k,i,j)==1)
        {
            tam=A3D_ELEM(vol_tam,k,i,j);
            COMPUTE_MEGV(VolxEdge,vol_edge,tam,k,i,j,pixel);
            A3D_ELEM(Vol_E,k,i,j)=pixel;
        }
    }
#ifdef DEBUG
    save()=Vol_E;
    save.write("PPP5_vol_megv.vol");
#endif

    //5.-Nonlinear function

    MultidimArray<double> vol_f=vol;
    double a, b, x, E;
    double threshold=lowerIntensity*255;

    FOR_ALL_ELEMENTS_IN_ARRAY3D(vol)
    {
        x=A3D_ELEM(vol,k,i,j);
        if (A3D_ELEM(mask,k,i,j)==1 && x>threshold)
        {
            E=A3D_ELEM(Vol_E,k,i,j);
            // We search the interval
            if ((x>=0)  && (x<=3))
            {
                a=0;
                b=3;
            }
            else if ((x>3)   && (x<=8))
            {
                a=3;
                b=8;
            }
            else if ((x>8)   && (x<=16))
            {
                a=8;
                b=16;
            }
            else if ((x>16)  && (x<=30))
            {
                a=16;
                b=30;
            }
            else if ((x>30)  && (x<=49))
            {
                a=30;
                b=49;
            }
            else if ((x>49)  && (x<=75))
            {
                a=49;
                b=75;
            }
            else if ((x>75)  && (x<=107))
            {
                a=75;
                b=107;
            }
            else if ((x>107) && (x<=147))
            {
                a=107;
                b=147;
            }
            else if ((x>147) && (x<=196))
            {
                a=147;
                b=196;
            }
            else
            {
                a=196;
                b=255;
            }
            // nonlinear function
            double x_new;
            if (x<=E)
            {
                x_new=E-sqrt(((E-a)*(E-a))-((x-a)*(x-a)));
            }
            if (x>E)
            {
                x_new=E+sqrt(((b-E)*(b-E))-((x-b)*(x-b)));
            }
            A3D_ELEM(vol_f,k,i,j)=x_new;
        }
    }
#ifdef DEBUG
    save()=vol_f;
    save.write("PPP6_vol_f.vol");
#endif

    // Unpadd volume with two new rows and cols
    if( FINISHINGZ(vol_f) - STARTINGZ(vol_f) > 0 )
    {
        vol_f.selfWindow( STARTINGZ(vol_f)+2,STARTINGY(vol_f)+2,STARTINGX(vol_f)+2,
                      FINISHINGZ(vol_f)-2, FINISHINGY(vol_f)-2,FINISHINGX(vol_f)-2 );
    }
    else
    {
        vol_f.selfWindow( STARTINGZ(vol_f),STARTINGY(vol_f)+2,STARTINGX(vol_f)+2,
                      FINISHINGZ(vol_f), FINISHINGY(vol_f)-2,FINISHINGX(vol_f)-2 );
    }

    //6.-Re-Scale (now we put the original values different from [0-255])
    double x_s,x_ns;
    FOR_ALL_ELEMENTS_IN_ARRAY3D(vol_f)
    {
        x_s=A3D_ELEM(vol_f,k,i,j);
        x_ns=((maxVal-minVal)/255)*(x_s+((255*minVal)/(maxVal-minVal)));
        A3D_ELEM(vol_f,k,i,j)=x_ns;
    }
    // OUTPUT
    vol=vol_f;
} // End of enhance

readParams()

void ProgVolumeEnhanceContrast::readParams ( )

protectedvirtual

Read parameters from command line

Reimplemented from XmippProgram.

Definition at line 51 of file volume_enhance_contrast.cpp.

{
    fnIn=getParam("-i");
    fnOut=getParam("-o");
    if (fnOut=="")
        fnOut=fnIn;
    alpha = getDoubleParam("--alpha");
    lowerIntensity = getDoubleParam("--lowerIntensity");
    removeBg = checkParam("--removeBackground");
    fnMask = getParam("--saveMask");
}

run()

void ProgVolumeEnhanceContrast::run ( )

protectedvirtual

This function will be start running the program. it also should be implemented by derived classes.

Reimplemented from XmippProgram.

Definition at line 63 of file volume_enhance_contrast.cpp.

{
  Image<double> img;
  img.read(fnIn);
  enhance(img());
  img.write(fnOut);
}

show()

void ProgVolumeEnhanceContrast::show

(

)

Show parameters. This function calls show_specific.

Definition at line 72 of file volume_enhance_contrast.cpp.

{
    if (verbose==0)
        return;
    std::cout
    << "Input: " << fnIn << std::endl
    << "Output: " << fnOut << std::endl
    << "Alpha: " << alpha << std::endl
    << "lowerIntensity: " << lowerIntensity << std::endl
    << "removeBackground: " << removeBg << std::endl
    << "saveMask: " << fnMask << std::endl;
}

Member Data Documentation

alpha

double ProgVolumeEnhanceContrast::alpha

Definition at line 54 of file volume_enhance_contrast.h.

fnIn

FileName ProgVolumeEnhanceContrast::fnIn

Definition at line 51 of file volume_enhance_contrast.h.

fnMask

FileName ProgVolumeEnhanceContrast::fnMask

Definition at line 63 of file volume_enhance_contrast.h.

fnOut

FileName ProgVolumeEnhanceContrast::fnOut

Definition at line 51 of file volume_enhance_contrast.h.

lowerIntensity

double ProgVolumeEnhanceContrast::lowerIntensity

Definition at line 57 of file volume_enhance_contrast.h.

removeBg

bool ProgVolumeEnhanceContrast::removeBg

Definition at line 60 of file volume_enhance_contrast.h.

---

The documentation for this class was generated from the following files:

• xmipp/legacy/libraries/reconstruction/volume_enhance_contrast.h
• xmipp/legacy/libraries/reconstruction/volume_enhance_contrast.cpp