POCSClass Class Reference

#include <recons_misc.h>

Collaboration diagram for POCSClass:

Public Types
enum	t_POCS_status {   POCS_measuring, POCS_use, POCS_lowering, POCS_N_measure,   POCS_N_use }

Public Member Functions
	POCSClass (BasicARTParameters *_prm, int _Zoutput_volume_size, int _Youtput_volume_size, int _Xoutput_volume_size)
	Constructor. More...
void	newIteration ()
	Start New ART iteration. More...
void	newProjection ()
	Start new Projection. More...
void	apply (GridVolume &vol_basis, int it, int images)
	Apply. More...

Public Attributes
t_POCS_status	POCS_state
double	POCS_avg
double	POCS_stddev
double	POCS_min
double	POCS_max
double	POCS_mean_error
double	POCS_max_error
double	POCS_global_mean_error
int	POCS_freq
int	POCS_i
int	POCS_vec_i
int	POCS_used
int	POCS_N
int	Zoutput_volume_size
int	Youtput_volume_size
int	Xoutput_volume_size
bool	apply_POCS
MultidimArray< double >	POCS_errors
BasicARTParameters *	prm

Detailed Description

POCS structure

Definition at line 111 of file recons_misc.h.

Member Enumeration Documentation

t_POCS_status

enum POCSClass::t_POCS_status

Enumerator
POCS_measuring
POCS_use
POCS_lowering
POCS_N_measure
POCS_N_use

Definition at line 114 of file recons_misc.h.

                 {POCS_measuring, POCS_use, POCS_lowering, POCS_N_measure,
                  POCS_N_use} t_POCS_status;

Constructor & Destructor Documentation

POCSClass()

POCSClass::POCSClass	(	BasicARTParameters *	_prm,
		int	_Zoutput_volume_size,
		int	_Youtput_volume_size,
		int	_Xoutput_volume_size
	)

Constructor.

Definition at line 699 of file recons_misc.cpp.

{
    prm = _prm;
    POCS_state = POCS_measuring;
    POCS_freq = prm->POCS_freq;
    POCS_i = 0;
    POCS_vec_i = 0;
    POCS_used = 0;
    POCS_N = 0;
    POCS_errors.initZeros(POCS_N_measure);
    Zoutput_volume_size = _Zoutput_volume_size;
    Youtput_volume_size = _Youtput_volume_size;
    Xoutput_volume_size = _Xoutput_volume_size;
    apply_POCS = (prm->surface_mask != nullptr ||
                  prm->positivity || (prm->force_sym != 0 && !prm->is_crystal) ||
                  prm->known_volume != -1);
}

Member Function Documentation

apply()

void POCSClass::apply	(	GridVolume &	vol_basis,
		int	it,
		int	images
	)

Apply.

Definition at line 730 of file recons_misc.cpp.

{
    Image<double> vol_POCS, theo_POCS_vol, corr_POCS_vol, vol_voxels;

    if (apply_POCS && POCS_i % POCS_freq == 0)
    {
        Image<double> vol_aux;
        Image<double> *desired_volume = nullptr;

        // Compute the corresponding voxel volume
        prm->basis.changeToVoxels(vol_basis, &(vol_voxels()),
                                  Zoutput_volume_size, Youtput_volume_size, Xoutput_volume_size);
        if (prm->tell&TELL_SAVE_AT_EACH_STEP)
        {
            vol_voxels.write("PPPvolPOCS0.vol");
            std::cout << "Stats PPPvolPOCS0.vol: ";
            vol_voxels().printStats();
            std::cout << std::endl;
        }
        // Apply surface restriction
        if (prm->surface_mask != nullptr)
        {
            vol_POCS() = (*(prm->surface_mask))();
        }
        else
        {
            vol_POCS().resize(vol_voxels());
            vol_POCS().initZeros();
        }
        if (prm->tell&TELL_SAVE_AT_EACH_STEP)
        {
            vol_POCS.write("PPPvolPOCS1.vol");
            std::cout << "Stats PPPvolPOCS1.vol: ";
            vol_POCS().printStats();
            std::cout << std::endl;
        }
        // Force symmetry
        if (prm->force_sym != 0)
        {
            symmetrizeVolume(prm->SL, vol_voxels(), vol_aux());
            desired_volume = &vol_aux;
            if (prm->tell&TELL_SAVE_AT_EACH_STEP)
            {
                vol_aux.write("PPPvolPOCS2.vol");
                std::cout << "Stats PPPvolPOCS2.vol: ";
                vol_aux().printStats();
                std::cout << std::endl;
            }
        }
        // Apply volume constraint
        if (prm->known_volume != -1)
        {
            Histogram1D hist;
            MultidimArray<int> aux_mask;
            aux_mask.resize(vol_POCS());
            FOR_ALL_ELEMENTS_IN_ARRAY3D(aux_mask)
            aux_mask(k, i, j) = 1 - (int)vol_POCS(k, i, j);
            long mask_voxels = vol_POCS().countThreshold("below", 0.5, 0);
            compute_hist_within_binary_mask(
                aux_mask, vol_voxels(), hist, 300);
            double known_percentage;
            known_percentage = XMIPP_MIN(100, 100 * prm->known_volume / mask_voxels);
            double threshold;
            threshold = hist.percentil(100 - known_percentage);
            FOR_ALL_ELEMENTS_IN_ARRAY3D(vol_voxels())
            if (vol_voxels(k, i, j) < threshold)
                vol_POCS(k, i, j) = 1;
        }
        if (prm->tell&TELL_SAVE_AT_EACH_STEP)
        {
            vol_POCS.write("PPPvolPOCS3.vol");
            std::cout << "Stats PPPvolPOCS3.vol: ";
            vol_POCS().printStats();
            std::cout << std::endl;
        }

        // Do not allow positivity outside interest region
        // and do not allow negativity inside the interest region
        // if positivity restrictions are to be applied
        // int bg = (int) vol_POCS().sum();
        // int fg = MULTIDIM_SIZE(vol_POCS()) - bg;
        int relax = 0, posi = 0;
        const MultidimArray<double> &mVolVoxels=vol_voxels();
        const MultidimArray<double> &mVolPOCS=vol_POCS();
        FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(mVolVoxels)
        if (DIRECT_A3D_ELEM(mVolPOCS,k, i, j) == 1 && DIRECT_A3D_ELEM(mVolVoxels,k, i, j) < 0)
        {
            DIRECT_A3D_ELEM(mVolPOCS,k, i, j) = 0;
            relax++;
        }
        else if (DIRECT_A3D_ELEM(mVolPOCS,k, i, j) == 0 && DIRECT_A3D_ELEM(mVolVoxels,k, i, j) < 0 && prm->positivity)
        {
            DIRECT_A3D_ELEM(mVolPOCS,k, i, j) = 1;
            posi++;
        }
        // Debugging messages
        //std::cout << "Relaxation/Positivity " << (double)relax/(double)bg << " "
        //     << (double)posi/(double)fg << " " << std::endl;

        // Solve volumetric equations
        switch (prm->basis.type)
        {
        case Basis::blobs:
            if (desired_volume == nullptr)
                ART_voxels2blobs_single_step(vol_basis, &vol_basis,
                                             prm->basis.blob, prm->D, prm->lambda(it),
                                             &(theo_POCS_vol()), nullptr,
                                             &(corr_POCS_vol()),
                                             &(vol_POCS()),
                                             POCS_mean_error, POCS_max_error, VARTK);
            else
            {
                FOR_ALL_ELEMENTS_IN_ARRAY3D(vol_POCS())
                if (vol_POCS(k, i, j) == 1)
                    (*desired_volume)(k, i, j) = 0;
                for (int i = 0; i < prm->force_sym; i++)
                {
                    ART_voxels2blobs_single_step(vol_basis, &vol_basis,
                                                 prm->basis.blob, prm->D, prm->lambda(it),
                                                 &(theo_POCS_vol()), &((*desired_volume)()),
                                                 &(corr_POCS_vol()),
                                                 nullptr,
                                                 POCS_mean_error, POCS_max_error, VARTK);
                    if (prm->tell&TELL_SAVE_AT_EACH_STEP)
                        std::cout << "    POCS Iteration " << i
                        << " POCS Error=" <<  POCS_mean_error << std::endl;
                }
            }
            break;
        case Basis::voxels:
            if (desired_volume == nullptr)
            {
                FOR_ALL_ELEMENTS_IN_ARRAY3D(vol_POCS())
                if (vol_POCS(k, i, j))
                    vol_basis(0)(k, i, j) = 0;
            }
            else
            {
                vol_basis(0)().initZeros();
                FOR_ALL_ELEMENTS_IN_ARRAY3D((*desired_volume)())
                vol_basis(0)(k, i, j) = (*desired_volume)(k, i, j);
            }
            POCS_mean_error = -1;
            break;
        default:
            REPORT_ERROR(ERR_ARG_INCORRECT,"This function cannot work with this basis");
        }
        POCS_i = 1;
        POCS_global_mean_error += POCS_mean_error;
        POCS_N++;

        // Now some control logic
        if (prm->numIMG - images < 100 || images % 100 == 0 ||
            desired_volume != nullptr)
        {
            POCS_freq = 1;
            POCS_state = POCS_measuring;
            POCS_vec_i = 0;
        }
        else
        {
            double dummy;
            switch (POCS_state)
            {
            case POCS_measuring:
#ifdef DEBUG_POCS

                std::cout << "M:" << POCS_vec_i << " " << POCS_mean_error << std::endl;
#endif

                POCS_errors(POCS_vec_i++) = POCS_mean_error;
                if (POCS_vec_i == POCS_N_measure)
                {
                    POCS_vec_i = 0;
                    // Change to use state
                    POCS_used = 0;
                    POCS_freq++;
                    POCS_state = POCS_use;
#ifdef DEBUG_POCS

                    std::cout << "1: Changing to " << POCS_freq << std::endl;
#endif

                }
                break;
            case POCS_use:
                POCS_used++;
                POCS_errors.computeStats(POCS_avg,
                                         POCS_stddev, dummy, POCS_min);
#ifdef DEBUG_POCS

                std::cout << "Reference errors: " << POCS_errors.transpose() << std::endl;
                std::cout << "Checking " << ABS(POCS_mean_error - POCS_avg) << " " << 1.2*1.96*POCS_stddev << std::endl;
#endif

                if (ABS(POCS_mean_error - POCS_avg) < 1.2*1.96*POCS_stddev)
                {
                    if (POCS_mean_error < POCS_avg)
                    {
                        double max_error = POCS_errors(0);
                        POCS_vec_i = 0;
                        for (int i = 1; i < POCS_N_measure; i++)
                            if (POCS_errors(i) > max_error)
                            {
                                max_error = POCS_errors(i);
                                POCS_vec_i = i;
                            }
                        POCS_errors(POCS_vec_i) = POCS_mean_error;
                    }
                    if (POCS_used < POCS_N_use)
                    { // While not enough uses
                    }
                    else if (POCS_freq < 3)
                    { // increase frequency
                        POCS_freq++;
#ifdef DEBUG_POCS

                        std::cout << "2: Changing to " << POCS_freq << std::endl;
#endif

                        POCS_used = 0;
                    }
                }
                else
                {
                    // It is behaving worse
                    if (POCS_freq > prm->POCS_freq + 1)
                    {
                        POCS_freq = prm->POCS_freq + 1;
                        POCS_used = 0;
#ifdef DEBUG_POCS

                        std::cout << "3: Changing to " << POCS_freq << std::endl;
#endif

                    }
                    else if (POCS_used > 2)
                    {
                        POCS_freq = prm->POCS_freq;
                        // Change status
                        POCS_used = 0;
                        POCS_state = POCS_lowering;
#ifdef DEBUG_POCS

                        std::cout << "Lowering\n";
#endif

                    }
                }
                break;
            case POCS_lowering:
                // Lower the POCS error before measuring again
                POCS_errors.computeStats(POCS_avg,
                                         POCS_stddev, POCS_max, dummy);
                POCS_used++;
                if (POCS_mean_error < POCS_max || POCS_used > 2*POCS_N_measure)
                {
                    // Change status
                    POCS_vec_i = 0;
                    POCS_state = POCS_measuring;
                }
                break;
            default:
                REPORT_ERROR(ERR_ARG_INCORRECT,"Unknown equation mode");
            }
        }
    }
    else
    {
        POCS_i++;
        POCS_mean_error = -1;
    }
}

newIteration()

void POCSClass::newIteration

(

)

Start New ART iteration.

Definition at line 719 of file recons_misc.cpp.

{
    POCS_global_mean_error = 0;
}

newProjection()

void POCSClass::newProjection

(

)

Start new Projection.

Definition at line 724 of file recons_misc.cpp.

{
    POCS_N = 0;
}

Member Data Documentation

apply_POCS

bool POCSClass::apply_POCS

Definition at line 132 of file recons_misc.h.

POCS_avg

double POCSClass::POCS_avg

Definition at line 117 of file recons_misc.h.

POCS_errors

MultidimArray<double> POCSClass::POCS_errors

Definition at line 133 of file recons_misc.h.

POCS_freq

int POCSClass::POCS_freq

Definition at line 124 of file recons_misc.h.

POCS_global_mean_error

double POCSClass::POCS_global_mean_error

Definition at line 123 of file recons_misc.h.

POCS_i

int POCSClass::POCS_i

Definition at line 125 of file recons_misc.h.

POCS_max

double POCSClass::POCS_max

Definition at line 120 of file recons_misc.h.

POCS_max_error

double POCSClass::POCS_max_error

Definition at line 122 of file recons_misc.h.

POCS_mean_error

double POCSClass::POCS_mean_error

Definition at line 121 of file recons_misc.h.

POCS_min

double POCSClass::POCS_min

Definition at line 119 of file recons_misc.h.

POCS_N

int POCSClass::POCS_N

Definition at line 128 of file recons_misc.h.

POCS_state

t_POCS_status POCSClass::POCS_state

Definition at line 116 of file recons_misc.h.

POCS_stddev

double POCSClass::POCS_stddev

Definition at line 118 of file recons_misc.h.

POCS_used

int POCSClass::POCS_used

Definition at line 127 of file recons_misc.h.

POCS_vec_i

int POCSClass::POCS_vec_i

Definition at line 126 of file recons_misc.h.

prm

BasicARTParameters* POCSClass::prm

Definition at line 134 of file recons_misc.h.

Xoutput_volume_size

int POCSClass::Xoutput_volume_size

Definition at line 131 of file recons_misc.h.

Youtput_volume_size

int POCSClass::Youtput_volume_size

Definition at line 130 of file recons_misc.h.

Zoutput_volume_size

int POCSClass::Zoutput_volume_size

Definition at line 129 of file recons_misc.h.

---

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

• xmipp/libraries/reconstruction/recons_misc.h
• xmipp/libraries/reconstruction/recons_misc.cpp