ProgResolutionIBW Class Reference

#include <resolution_ibw.h>

Inheritance diagram for ProgResolutionIBW:

Collaboration diagram for ProgResolutionIBW:

Public Member Functions
void	readParams ()
void	defineParams ()
void	show () const
void	run ()
void	edgeWidth (const MultidimArray< double > &volCoeffs, const MultidimArray< double > &edges, MultidimArray< double > &widths, const Matrix1D< double > &dir, double step=0.25) const
double	calculateIBW (MultidimArray< double > &widths) const
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	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	fnVol
	Input volume. More...
FileName	fnOut
	Output volume with local widths. More...
Image< double >	V
	Volume to evaluate. More...
Public Attributes inherited from XmippProgram
bool	doRun
bool	runWithoutArgs
int	verbose
	Verbosity level. More...
int	debug

Additional Inherited Members
Protected Member Functions inherited from XmippProgram
void	defineCommons ()
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

Parameter class for the resolution IBW program. The program is based on P. Marziliano, F. Dufaux, S. Winkler, T. Ebrahimi. Perceptual blur and ringing metrics: application to JPEG2000. Signal Processing: Image Communication, 19: 163-172 (2004)

Definition at line 42 of file resolution_ibw.h.

Member Function Documentation

calculateIBW()

double ProgResolutionIBW::calculateIBW

(

MultidimArray< double > &

widths

)

const

Calculate the inverse border width. The calculation is only performed at the border pixels.

Definition at line 307 of file resolution_ibw.cpp.

{
    double total, count;
    total=count=0;

    FOR_ALL_ELEMENTS_IN_ARRAY3D(widths)
    {
        double width=A3D_ELEM(widths,k,i,j);
        if (width<1e4)
        {
            total+=width;
            count++;
        }
    }
    double avg = total/count;
    return 1.0/avg;
}

defineParams()

void ProgResolutionIBW::defineParams ( )

virtual

Define parameters.

Reimplemented from XmippProgram.

Definition at line 39 of file resolution_ibw.cpp.

{
    addUsageLine("Evaluate the resolution of a volume through the inverse border widths");
    addParamsLine("   -i <file>              : Volume to evaluate");
    addParamsLine("  [-o <file=\"\">]        : Volume with the border widths of the edge voxels");
}

edgeWidth()

void ProgResolutionIBW::edgeWidth	(	const MultidimArray< double > &	volCoeffs,
		const MultidimArray< double > &	edges,
		MultidimArray< double > &	widths,
		const Matrix1D< double > &	dir,
		double	step = 0.25
	)		const

Compute the edge width at the edge pixels in a given direction. step controls how fine is the search of the border width.

Definition at line 205 of file resolution_ibw.cpp.

{
    double forward_count, backward_count, slope;
    Matrix1D<double> pos_aux_fw(3), pos_aux_bw(3), pos(3), pos_aux(3), next_pos(3), Kdir;

    Kdir=step*dir;

    //Visit all elements in volume
    FOR_ALL_ELEMENTS_IN_ARRAY3D(edges)
    {
        //Check for border pixels
        if (A3D_ELEM(edges,k,i,j)!=0)
        {
            //reset all counters
            forward_count=0;
            backward_count=0;
            VECTOR_R3(pos_aux_fw,j,i,k);
            pos_aux_bw=pos=pos_aux_fw;

            //find out if pixel magnitude grows or decreases
            pos_aux=pos;
            pos_aux+=dir;
            double value_plus_dir=volCoeffs.interpolatedElementBSpline3D(XX(pos_aux),YY(pos_aux),ZZ(pos_aux));

            pos_aux=pos;
            pos_aux-=dir;
            double value_minus_dir=volCoeffs.interpolatedElementBSpline3D(XX(pos_aux),YY(pos_aux),ZZ(pos_aux));

            slope=value_plus_dir-value_minus_dir;

            double sign;
            if (slope>0)
                sign=1;
            else
                sign=-1;

            //current_pixel is multiplied by the sign, so only one condition is enough to detect an
            //extremum no matter if the pixel values increase or decrease
            double current_pixel=sign*volCoeffs.interpolatedElementBSpline3D
                                 (XX(pos_aux_fw),YY(pos_aux_fw),ZZ(pos_aux_fw));

            double next_pixel;
            bool not_found;

            //Search for local extremum ahead of the edge in the given direction
            do
            {
                not_found=true;
                next_pos=pos_aux_fw+Kdir;
                next_pixel=sign*volCoeffs.interpolatedElementBSpline3D
                           (XX(next_pos),YY(next_pos),ZZ(next_pos));

                if(next_pixel>current_pixel)
                {
                    current_pixel=next_pixel;
                    pos_aux_fw=next_pos;
                    forward_count++;
                }
                else
                {
                    not_found=false;
                }
            }
            while(not_found);

            current_pixel=sign*volCoeffs.interpolatedElementBSpline3D
                          (XX(pos_aux_bw),YY(pos_aux_bw),ZZ(pos_aux_bw));

            //Search for local extremum behind of the edge in the given direction
            do
            {
                not_found=true;
                next_pos=pos_aux_bw-Kdir;
                next_pixel=sign*volCoeffs.interpolatedElementBSpline3D
                           (XX(next_pos),YY(next_pos),ZZ(next_pos));

                if(next_pixel<current_pixel)
                {
                    current_pixel=next_pixel;
                    pos_aux_bw=next_pos;
                    backward_count++;
                }
                else
                {
                    not_found=false;
                }
            }
            while(not_found);

            //If the width found for this position is smaller than the one stores in edges volume
            //before it is overwritten
            if ((forward_count+backward_count)<A3D_ELEM(widths,k,i,j))
            {
                A3D_ELEM(widths,k,i,j)=forward_count+backward_count;
            }
        }
    }
}

readParams()

void ProgResolutionIBW::readParams ( )

virtual

Read from a command line.

Reimplemented from XmippProgram.

Definition at line 32 of file resolution_ibw.cpp.

{
    fnVol = getParam("-i");
    fnOut = getParam("-o");
}

run()

void ProgResolutionIBW::run ( )

virtual

Compute the correlation for all micrographs

Reimplemented from XmippProgram.

Definition at line 58 of file resolution_ibw.cpp.

{
    V.read(fnVol);

    //Mask generation
    Image<double> aux;
    double bg_mean;
    MultidimArray<double> Vmask;
    detectBackground(V(),aux(),0.1,bg_mean);
#ifdef DEBUG

    aux.write("PPPmask_no_ero_03.vol");
#endif

    //Mask volume erosion to expand the mask boundaries
    Vmask.initZeros(V());
    erode3D(aux(),Vmask, 18,0,2);

    //Correction of some flaws produced in the edges of the mask volume
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(Vmask)
    if (k<=4 || i<=4 || j<=4 ||
        k>=ZSIZE(Vmask)-4 || i>=YSIZE(Vmask)-4 || j>=XSIZE(Vmask)-4)
        DIRECT_A3D_ELEM(Vmask,k,i,j)=1;

    aux()=Vmask;
#ifdef DEBUG

    aux.write("PPPmask_ero_03.vol");
#endif

    //Sobel edge detection applied to original volume
    Image<double> Vedge;
    computeEdges(V(),Vedge());
#ifdef DEBUG

    Vedge.write("PPPvolume_sobel_unmask_03.vol");
#endif

    //Masked volume generation
    const MultidimArray<double> &mVedge=Vedge();
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(mVedge)
    if (DIRECT_MULTIDIM_ELEM(Vmask,n)==1)
        DIRECT_MULTIDIM_ELEM(mVedge,n)=0;
#ifdef DEBUG

    Vedge.write("volume_sobel_mask_03.vol");
#endif

    double minval, maxval, avg, stddev;

    //Invert the mask to meet computeStats_within_binary_mask requirements
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(Vmask)
    if (DIRECT_MULTIDIM_ELEM(Vmask,n)==1)
        DIRECT_MULTIDIM_ELEM(Vmask,n)=0;
    else
        DIRECT_MULTIDIM_ELEM(Vmask,n)=1;

    //Threshold is 3 times the standard deviation of unmasked pixel values
    double thresh;
    computeStats_within_binary_mask(Vmask,mVedge,minval, maxval, avg, stddev);
    thresh=3*stddev;

    //Final edge volume generated by setting to 1 positions with values > threshold
    Image<double> Vaux;
    Vaux().initZeros(mVedge);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(mVedge)
    if (DIRECT_MULTIDIM_ELEM(mVedge,n)>=thresh)
        DIRECT_MULTIDIM_ELEM(Vaux(),n)=1;

#ifdef DEBUG

    Vaux.write("volumen_bordes_definitivo_03.vol");
#endif

    const MultidimArray<double> &mVaux=Vaux();

    //Spline coefficient volume from original volume, to allow <1 step sizes
    MultidimArray<double> Volcoeffs;
    Volcoeffs.initZeros(V());

    produceSplineCoefficients(3,Volcoeffs,V());

    //Width parameter volume initialization
    Image<double> widths;
    widths().resizeNoCopy(V());
    widths().initConstant(1e5);
    double step=0.25;

    Matrix1D<double> direction(3);

    //Calculation of edge width for 10 different directions, if a smaller value is found for a different
    //direction on a given position the former value is overwritten

    //Direction (1,0,0)
    VECTOR_R3(direction,1,0,0);
    edgeWidth(Volcoeffs, mVaux, widths(), direction, step);

    //Direction (0,1,0)
    VECTOR_R3(direction,0,1,0);
    edgeWidth(Volcoeffs, mVaux, widths(), direction, step);

    //Direction (0,0,1)
    VECTOR_R3(direction,0,0,1);
    edgeWidth(Volcoeffs, mVaux, widths(), direction, step);

    //Direction (1,1,0)
    VECTOR_R3(direction,(1/sqrt(2)),(1/sqrt(2)),0);
    edgeWidth(Volcoeffs, mVaux, widths(), direction, step);

    //Direction (1,0,1)
    VECTOR_R3(direction,(1/sqrt(2)),0,(1/sqrt(2)));
    edgeWidth(Volcoeffs, mVaux, widths(), direction, step);

    //Direction (0,1,1)
    VECTOR_R3(direction,0,(1/sqrt(2)),(1/sqrt(2)));
    edgeWidth(Volcoeffs, mVaux, widths(), direction, step);

    //Direction (1,1,1)
    VECTOR_R3(direction,(1/sqrt(3)),(1/sqrt(3)),(1/sqrt(3)));
    edgeWidth(Volcoeffs, mVaux, widths(), direction, step);

    //Direction (-1,1,1)
    VECTOR_R3(direction,-(1/sqrt(3)),(1/sqrt(3)),(1/sqrt(3)));
    edgeWidth(Volcoeffs, mVaux, widths(), direction, step);

    //Direction (1,1,-1)
    VECTOR_R3(direction,(1/sqrt(3)),(1/sqrt(3)),-(1/sqrt(3)));
    edgeWidth(Volcoeffs, mVaux, widths(), direction, step);

    //Direction (1,-1,1)
    VECTOR_R3(direction,(1/sqrt(3)),-(1/sqrt(3)),(1/sqrt(3)));
    edgeWidth(Volcoeffs, mVaux, widths(), direction, step);

#ifdef DEBUG

    std::cout << "width stats: ";
    widths().printStats();
    std::cout << std::endl;
    widths.write("PPPwidths.vol");
#endif

    double ibw=calculateIBW(widths());
    std::cout << "Resolution ibw= " << ibw << std::endl;
    if (fnOut!="")
        widths.write(fnOut);
}

show()

void ProgResolutionIBW::show ( ) const

virtual

Show parameters.

Reimplemented from XmippProgram.

Definition at line 47 of file resolution_ibw.cpp.

{
    if (verbose==0)
        return;
    std::cout << "Input volume:      " << fnVol << std::endl
              << "Output widths:     " << fnOut << std::endl
    ;
}

Member Data Documentation

fnOut

FileName ProgResolutionIBW::fnOut

Output volume with local widths.

Definition at line 49 of file resolution_ibw.h.

fnVol

FileName ProgResolutionIBW::fnVol

Input volume.

Definition at line 46 of file resolution_ibw.h.

V

Image<double> ProgResolutionIBW::V

Volume to evaluate.

Definition at line 53 of file resolution_ibw.h.

---

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

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