volume_segment.cpp File Reference

#include <core/args.h>
#include <data/morphology.h>
#include <data/filters.h>
#include "volume_segment.h"

Include dependency graph for volume_segment.cpp:

Go to the source code of this file.

Functions
double	segment_threshold (const Image< double > *V_in, Image< double > *V_out, double threshold, bool do_prob)
void	wang_smoothing (const Image< double > *V_in, Image< double > *V_out, int radius)
void	probabilistic_solvent (Image< double > *V_in, Image< double > *V_out)

Function Documentation

probabilistic_solvent()

void probabilistic_solvent	(	Image< double > *	V_in,
		Image< double > *	V_out
	)

Definition at line 233 of file volume_segment.cpp.

{
    // Calculate mean and sigma for protein and solvent regions
    // according to the traditional segmentation
    double Np, sump, sum2p, Ns, sums, sum2s;
    double avgp, sigp, avgs, sigs, aux, solv_frac, prot_frac;
    double p_prot, p_solv;

    MultidimArray<double> &mVin=(*V_in)();
    MultidimArray<double> &mVout=(*V_out)();
    mVin.setXmippOrigin();
    mVout.setXmippOrigin();

    Np = sump = sum2p = Ns = sums = sum2s = 0.;
    FOR_ALL_ELEMENTS_IN_ARRAY3D(mVin)
    {
        aux = A3D_ELEM(mVin, k, i, j);
        if (A3D_ELEM(mVout, k, i, j) < 0.5)
        {
            sums += aux;
            sum2s += aux * aux;
            Ns += 1.;
        }
        else
        {
            sump += aux;
            sum2p += aux * aux;
            Np += 1.;
        }
    }
    if (Np > 0. && Ns > 0.)
    {
        avgs = sums / Ns;
        sigs = sum2s / Ns - avgs * avgs;
        avgp = sump / Np;
        sigp = sum2p / Np - avgp * avgp;
        prot_frac = Np / (Np + Ns);
        solv_frac = 1. - prot_frac;
    }
    else
    {
        REPORT_ERROR(ERR_NUMERICAL, "empty solvent or protein region");
    }

    // Terwilliger-like calculation of P(x|solv) & P(x|prot)
    // Bayes: P(prot|x)= P(x|prot)/{P(x|prot)+P(x|solv)}
    double isigs=-1/(2.0*sigs);
    double isigp=-1/(2.0*sigp);
    FOR_ALL_ELEMENTS_IN_ARRAY3D(mVin)
    {
        double voxelVal=A3D_ELEM(mVin, k, i, j);
        aux = voxelVal - avgs;
        p_solv = solv_frac * exp(aux * aux * isigs);
        aux = voxelVal - avgp;
        p_prot = prot_frac * exp(aux * aux * isigp);
        A3D_ELEM(mVout, k, i, j) = p_prot / (p_prot + p_solv);
    }
}

segment_threshold()

double segment_threshold	(	const Image< double > *	V_in,
		Image< double > *	V_out,
		double	threshold,
		bool	do_prob
	)

Definition at line 128 of file volume_segment.cpp.

{
    Image<double> aux;

    // Binarize input volume
    (*V_out)() = (*V_in)();
    (*V_out)().threshold("below", threshold, threshold);
    (*V_out)().binarize(threshold);

#ifdef DEBUG

    std::cout << threshold << std::endl;
    Image<double> save;
    save() = (*V_in)();
    save.write("PPP0.vol");
    save() = (*V_out)();
    save.write("PPP1.vol");
#endif

    if (!do_prob)
    {
        // Apply morphological opening to input volume
        aux().initZeros((*V_out)());
        opening3D((*V_out)(), aux(), 18, 0, 1);
        closing3D(aux(), (*V_out)(), 18, 0, 1);
#ifdef DEBUG

        save() = (*V_out)();
        save.write("PPP2.vol");
#endif
    }

    // Count the number of different objects
    int no_comp = labelImage3D((*V_out)(), aux());
    Matrix1D<double> count(no_comp + 1);
    const MultidimArray<double> &maux=aux();
    FOR_ALL_ELEMENTS_IN_ARRAY3D(maux)
    VEC_ELEM(count,(int)A3D_ELEM(maux,k, i, j))++;

#ifdef DEBUG

    std::cout << count << std::endl << std::endl;
    std::cout << "Press any key\n";
    char c;
    std::cin >> c;
#endif

    // Pick the maximum
    count(0) = 0; // We don't want to pick the background
    int imax=count.maxIndex();

    // Select the mask with only that piece
    const MultidimArray<double> &mVout=(*V_out)();
    FOR_ALL_ELEMENTS_IN_ARRAY3D(mVout)
    A3D_ELEM(mVout,k, i, j) = A3D_ELEM(maux,k, i, j) == imax;

    return count(imax);
}

wang_smoothing()

void wang_smoothing	(	const Image< double > *	V_in,
		Image< double > *	V_out,
		int	radius
	)

Definition at line 188 of file volume_segment.cpp.

{

    int r2;
    double radius2 = radius * radius;
    double sumw, weight;

    (*V_out)().initZeros((*V_in)());

    const MultidimArray<double> &mVin=(*V_in)();
    const MultidimArray<double> &mVout=(*V_out)();
    FOR_ALL_ELEMENTS_IN_ARRAY3D(mVin)
    {
        sumw = 0.;
        for (int kp = k - radius; kp < k + radius; kp++)
        {
            if (kp > STARTINGZ(mVin) && kp < FINISHINGZ(mVin))
            {
                for (int ip = i - radius; ip < i + radius; ip++)
                {
                    if (ip > STARTINGY(mVin) && ip < FINISHINGY(mVin))
                    {
                        for (int jp = j - radius; jp < j + radius; jp++)
                        {
                            if (jp > STARTINGX(mVin) && jp < FINISHINGX(mVin))
                            {
                                r2 = (kp - k) * (kp - k) + (ip - i) * (ip - i) + (jp - j) * (jp - j);
                                if ((r2 < radius2) && (A3D_ELEM(mVin, kp, ip, jp) > 0.))
                                {
                                    weight = 1. - sqrt(r2 / radius2);
                                    A3D_ELEM(mVout, k, i, j) += weight * XMIPP_MAX(0., A3D_ELEM(mVin, kp, ip, jp));
                                    sumw += weight;
                                }
                            }
                        }
                    }
                }
            }
        }
        if (sumw > 0.)
            A3D_ELEM(mVout, k, i, j) /= sumw;
    }
}