ARTReconsBase Class Reference

#include <base_art_recons.h>

Inheritance diagram for ARTReconsBase:

Collaboration diagram for ARTReconsBase:

Public Member Functions
virtual	~ARTReconsBase ()
virtual void	readParams (XmippProgram *program)
virtual void	print (std::ostream &o) const
virtual void	preProcess (GridVolume &vol_basis0, int level=FULL, int rank=-1)
	Produce Plain side information from the Class parameters. More...
virtual void	singleStep (GridVolume &vol_in, GridVolume *vol_out, Projection &theo_proj, Projection &read_proj, int sym_no, Projection &diff_proj, Projection &corr_proj, Projection &alig_proj, double &mean_error, int numIMG, double lambda, int act_proj, const FileName &fn_ctf, const MultidimArray< int > *maskPtr, bool refine)
virtual void	postProcess (GridVolume &vol_basis)
virtual void	applySymmetry (GridVolume &vol_in, GridVolume *vol_out, int grid_type)
void	initHistory (const GridVolume &vol_basis0)
void	iterations (GridVolume &vol_basis, int rank=-1)

Static Public Member Functions
static void	defineParams (XmippProgram *program, bool mpiMode=false)

Public Attributes
BasicARTParameters	artPrm

Friends
std::ostream &	operator<< (std::ostream &o, const ARTReconsBase &artRecons)

Detailed Description

ART Base Reconstruction. This class contains all the basic routines needed about the ART reconstruction process. This is the class used in order to reconstruct single particles. Any other type of reconstruction must inherit from this.

Definition at line 50 of file base_art_recons.h.

Constructor & Destructor Documentation

~ARTReconsBase()

virtual ARTReconsBase::~ARTReconsBase ( )

inlinevirtual

Definition at line 55 of file base_art_recons.h.

    {}

Member Function Documentation

applySymmetry()

void ARTReconsBase::applySymmetry ( GridVolume &  vol_in, GridVolume *  vol_out, int  grid_type  )

virtual

Force the trial volume to be symmetric. So far only implemented for crystals.

Reimplemented in CrystalARTRecons.

Definition at line 790 of file base_art_recons.cpp.

{
    REPORT_ERROR(ERR_NOT_IMPLEMENTED, "ARTReconsBase::applySymmetry: Function not implemented for single particles");
}

defineParams()

static void ARTReconsBase::defineParams ( XmippProgram *  program, bool  mpiMode = false  )

inlinestatic

Definition at line 58 of file base_art_recons.h.

    {
        BasicARTParameters::defineParams(program, mpiMode);
    }

initHistory()

void ARTReconsBase::initHistory

(

const GridVolume &

vol_basis0

)

Write first part of ART history. This function writes all ART parameters, projection angles, symmetry matrices, and grid structure in the history handler provided by the side information. At the end the routine makes a call to the operator << of the Extra_ART_Parameters to show the specific part of the History.

BasicARTParameters is not constant since things are written in BasicARTParameters::fh_hist.

Definition at line 607 of file base_art_recons.cpp.

{
    // Show general information ................................................
    *artPrm.fh_hist << " ============================================================\n";
    *artPrm.fh_hist << " ART RECONSTRUCTION HISTORY: \n";
    *artPrm.fh_hist << " ============================================================\n\n";
    *artPrm.fh_hist << " Parameters -------------------------------------------------\n";
    *artPrm.fh_hist << " Projections file     : " << artPrm.fn_sel << std::endl;
    *artPrm.fh_hist << " CTF file             : " << artPrm.fn_ctf << std::endl;
    *artPrm.fh_hist << " Goldmask             : " << artPrm.goldmask << std::endl;
    *artPrm.fh_hist << " Unmatched projectors : " << artPrm.unmatched << std::endl;
    *artPrm.fh_hist << " Sampling             : " << artPrm.sampling << std::endl;
    *artPrm.fh_hist << artPrm.basis << std::endl;
    switch (artPrm.grid_type)
    {
    case FCC:
        *artPrm.fh_hist << " Grid Type: FCC ";
        break;
    case BCC:
        *artPrm.fh_hist << " Grid Type: BCC ";
        break;
    case CC:
        *artPrm.fh_hist << " Grid Type: CC ";
        break;
    }
    *artPrm.fh_hist << " Shifted tomograms:" << artPrm.shiftedTomograms << std::endl;
    *artPrm.fh_hist << " Ray length: " << artPrm.ray_length << std::endl;
    *artPrm.fh_hist << "\n Radius of the interest sphere= " << artPrm.R
    << " pixels" << std::endl;
    *artPrm.fh_hist << " Grid unit=" << artPrm.grid_relative_size
    << " pixels" << std::endl;
    if (artPrm.proj_ext==0)
        *artPrm.fh_hist << " No Projection extension\n";
    else
        *artPrm.fh_hist << " Projection extension frame: " << artPrm.proj_ext
        << " pixels\n";
    if (artPrm.Zoutput_volume_size==0)
        *artPrm.fh_hist << " Output volume size as input images\n";
    else
        *artPrm.fh_hist << " Output volume size (ZxYxX): "
        << artPrm.Zoutput_volume_size << "x" << artPrm.Youtput_volume_size
        << "x" << artPrm.Xoutput_volume_size << std::endl;
    *artPrm.fh_hist << " Iterations:    lambda=" << artPrm.lambda_list << std::endl
    << "                No. Iter=" << artPrm.no_it << std::endl;
    if (artPrm.WLS)
    {
        *artPrm.fh_hist << " Perform weighted least-squares ART "<< std::endl;
        *artPrm.fh_hist << " Iterations:    kappa=" << artPrm.kappa_list << std::endl
        << "                No. Iter=" << artPrm.no_it << std::endl;
    }
    *artPrm.fh_hist << " Parallel mode: ";
    switch (artPrm.parallel_mode)
    {
    case BasicARTParameters::ART:
        *artPrm.fh_hist << "ART\n";
        break;
    case BasicARTParameters::SIRT:
        *artPrm.fh_hist << "SIRT\n";
        break;
    case BasicARTParameters::pSIRT:
        *artPrm.fh_hist << "Parallel SIRT\n";
        break;
    case BasicARTParameters::pfSIRT:
        *artPrm.fh_hist << "Parallel False SIRT\n";
        break;
    case BasicARTParameters::pSART:
        *artPrm.fh_hist << "SART, block size=" << artPrm.block_size << std::endl;
        break;
    case BasicARTParameters::pAVSP:
        *artPrm.fh_hist << "AVSP\n";
        break;
    case BasicARTParameters::pBiCAV:
        *artPrm.fh_hist << "BiCAV, block size=" << artPrm.block_size << std::endl;
        break;
    case BasicARTParameters::pCAV:
        *artPrm.fh_hist << "CAV (global algorithm)\n";
        break;
    }
    *artPrm.fh_hist << " Equation mode: ";
    if (artPrm.eq_mode==CAV)
        *artPrm.fh_hist << "CAV (Projection update)\n";
    else if (artPrm.eq_mode==CAVK)
        *artPrm.fh_hist << "CAVK\n";
    else if (artPrm.eq_mode==CAVARTK)
        *artPrm.fh_hist << "CAVARTK\n";
    else
        *artPrm.fh_hist << "ARTK\n";
    *artPrm.fh_hist << " Projections: Ydim x Xdim = " << artPrm.projYdim
    << " x " << artPrm.projXdim << std::endl;
    *artPrm.fh_hist << " Surface mask: " << artPrm.fn_surface_mask << std::endl;
    *artPrm.fh_hist << " POCS freq: " << artPrm.POCS_freq << std::endl;
    *artPrm.fh_hist << " Known volume: " << artPrm.known_volume << std::endl;
    *artPrm.fh_hist << " Positivity: " <<artPrm.positivity << std::endl;
    *artPrm.fh_hist << " Apply shifts in images headers: " << artPrm.apply_shifts << std::endl;
    *artPrm.fh_hist << " Symmetry file: " << artPrm.fn_sym << std::endl;
    *artPrm.fh_hist << " Force Symmetry each: " << artPrm.sym_each << " projections"
    <<std::endl;
    *artPrm.fh_hist << " Maximun absolute tilt angle: " << artPrm.max_tilt << " degrees"
    <<std::endl;
    *artPrm.fh_hist << " Generating symmetry group: " << !artPrm.do_not_generate_subgroup << std::endl;
    *artPrm.fh_hist << " Do not use symmetrized projections: "
    << !artPrm.do_not_use_symproj << std::endl;
    *artPrm.fh_hist << " Number of total projections (including symmetrized): "
    << artPrm.numIMG << std::endl;
    *artPrm.fh_hist << " Number of different projections: " << artPrm.trueIMG << std::endl;
    *artPrm.fh_hist << " Forcing symmetry: " << artPrm.force_sym << std::endl;
    *artPrm.fh_hist << " Stop at: " << artPrm.stop_at << std::endl;
    if (artPrm.random_sort)
        *artPrm.fh_hist << " Random sort" << std::endl;
    else
        *artPrm.fh_hist << " Sort with last " << artPrm.sort_last_N << " images\n";
    *artPrm.fh_hist << " Variability analysis: " << artPrm.variability_analysis << std::endl;
    *artPrm.fh_hist << " Refine projections: " << artPrm.refine << std::endl;
    *artPrm.fh_hist << " Sparsity epsilon: " << artPrm.sparseEps << std::endl;
    *artPrm.fh_hist << " Diffusion weight: " << artPrm.diffusionWeight << std::endl;
    if (artPrm.SL.symsNo()!=0)
    {
        Matrix2D<double> L(4,4),R(4,4); // A matrix from the list
        *artPrm.fh_hist << " Symmetry matrices -------\n";
        for (int j=0; j<artPrm.SL.symsNo(); j++)
        {
            artPrm.SL.getMatrices(j,L,R);
            *artPrm.fh_hist << " Left  Symmetry matrix " << j << std::endl << L;
            *artPrm.fh_hist << " Right Symmetry matrix " << j << std::endl << R << std::endl;
        }
    }
    *artPrm.fh_hist << " Saving intermidiate at every "
    << artPrm.save_intermidiate_every << " projections\n";
    if(artPrm.ref_trans_step > 0.)
    {
        *artPrm.fh_hist << " Refine translational alignement after "
        << artPrm.ref_trans_after << " projection presentations\n"
        << " Maximun allowed shift is: " << artPrm.ref_trans_step << "\n\n";
    }

    // Show angles .............................................................
    // Prepare info structure for showing
    MetaDataVec MD;
    double dfrot, dftilt, dfpsi;
    size_t id;
    for (int i=0; i<artPrm.numIMG; i++)
    {
        id=MD.addObject();
        MD.setValue(MDL_ANGLE_ROT, artPrm.IMG_Inf[i].rot,id);
        MD.setValue(MDL_ANGLE_TILT, artPrm.IMG_Inf[i].tilt,id);
        MD.setValue(MDL_ANGLE_PSI, artPrm.IMG_Inf[i].psi,id);
        MD.setValue(MDL_SYMNO, artPrm.IMG_Inf[i].sym,id);
    }

    // Now show
    *artPrm.fh_hist << " Projection angles -----------------------------------------\n";
    for (size_t objId : MD.ids())
    {
        MD.getValue(MDL_ANGLE_ROT, dfrot, objId);
        MD.getValue(MDL_ANGLE_TILT, dftilt, objId);
        MD.getValue(MDL_ANGLE_PSI, dfpsi, objId);

        *artPrm.fh_hist << "rot= "<<dfrot<<" tilt= "<<dftilt<<" psi= "<<dfpsi<<std::endl;
    }
    *artPrm.fh_hist << " -----------------------------------------------------------\n";

    // Show Initial volume and volume structure ................................
    if (artPrm.fn_start != "")
        *artPrm.fh_hist << " Starting from file: " << artPrm.fn_start << std::endl;
    else
        *artPrm.fh_hist << " Starting from a zero volume\n";

    *artPrm.fh_hist << " Grid structure ------\n" << vol_basis0.grid();

    // Show extra information ..................................................
    if (typeid(*this) != typeid(ARTReconsBase))
        *artPrm.fh_hist << " Extra: ";
    print(*artPrm.fh_hist);
}

iterations()

void ARTReconsBase::iterations	(	GridVolume &	vol_basis,
		int	rank = -1
	)

Perform all ART iterations. This function performs the iterations according to the ART parameters, it needs the side information to be fully computed. It throws a lot of information to the screen and to the history file (side.fh_hist), specially this one must exist.

The GridVolume must have as input an initial guess for the solution, and when this function finishes, it contains the final solution volume in basis.

The rank is the identification number of the process running this function. If it is -1, the function is run in sequential mode. If it is 0, then it is the root process.

See the BasicARTParameters for more information about how to generate the iterations.

This method is not virtual as it should be common for all ARTRecons classes.

Definition at line 45 of file base_art_recons.cpp.

{
    // Some variables .......................................................
    int        ART_numIMG;              // Normalizing factor
    GridVolume *ptr_vol_out;            // Pointer to output volume
    GridVolume vol_basis_out;           // Output volume (only useful in SIRT)
    Image<double> vol_voxels;             // This one is useful only in the
    // case of saving intermediate volumes
    int Xoutput_volume_size, Youtput_volume_size, Zoutput_volume_size;
    double aux_tilt;

    // Projection related ...................................................
    Projection      read_proj;          // Projection read from file
    Projection      theo_proj;          // Projection from the
    // reconstruction
    Projection      alig_proj;          // Projection with the correlation
    // maps needed for translation alignment
    Projection      corr_proj;          // Image with the correction
    // factors for unitary basis
    Projection      diff_proj;          // Difference between the
    // theoretical and real image

    //    preIterations(vol_basis);

    // Reconstruction results ...............................................
    double          mean_error = 0.0;
    double          global_mean_error,global_mean_error_1stblock;

    // Initialize residual image vector for wlsART ..........................
    if (artPrm.WLS)
    {
        artPrm.residual_imgs.clear();
        for (int iact_proj = 0; iact_proj < artPrm.numIMG ; iact_proj++)
        {
            read_proj.read(artPrm.IMG_Inf[iact_proj].fn_proj);
            read_proj().setXmippOrigin();
            read_proj().initZeros();
            artPrm.residual_imgs.push_back(read_proj);
        }
    }

    // Some initialization ..................................................
    Xoutput_volume_size = (artPrm.Xoutput_volume_size==0) ?
                          artPrm.projXdim : artPrm.Xoutput_volume_size;
    Youtput_volume_size = (artPrm.Youtput_volume_size==0) ?
                          artPrm.projYdim : artPrm.Youtput_volume_size;
    Zoutput_volume_size = (artPrm.Zoutput_volume_size==0) ?
                          artPrm.projXdim : artPrm.Zoutput_volume_size;

    // POCS constraints .....................................................
    POCSClass POCS(&artPrm, Zoutput_volume_size, Youtput_volume_size,
                   Xoutput_volume_size);

    // Variability analysis .................................................
    VariabilityClass VC(&artPrm, Zoutput_volume_size, Youtput_volume_size,
                        Xoutput_volume_size);

    // Noisy reconstruction .................................................
    GridVolume vol_basis_noisy; // Output volume (only for ART)
    Projection noisy_projection;
    MetaDataVec SF_noise, SF_signal;
    if (artPrm.noisy_reconstruction)
    {
        vol_basis_noisy = vol_basis;
        vol_basis_noisy.initZeros();
    }

    // If SIRT set normalizing factor and create output volume ..............
    if (artPrm.parallel_mode==BasicARTParameters::SIRT )
    {
        ART_numIMG = artPrm.numIMG;
        vol_basis_out = vol_basis;         // Copy the structure of vol_basis
        ptr_vol_out = &vol_basis_out;      // Pointer to output volume
        if (artPrm.variability_analysis)
            ptr_vol_out->initZeros();
    }
    else if ( artPrm.parallel_mode==BasicARTParameters::pfSIRT ||
              artPrm.parallel_mode==BasicARTParameters::pSIRT ||
              artPrm.parallel_mode==BasicARTParameters::pSART ||
              artPrm.parallel_mode==BasicARTParameters::pCAV )
    {
        // In those cases, normalization must be done at the top level program once we
        // have the reconstruction values. Have a look ar /Applications/Src/MPIArt/mpi_art.cc for
        // an example
        ART_numIMG = 1;
        ptr_vol_out = &vol_basis_out;      // Pointer to output volume
        vol_basis_out = vol_basis;         // Copy the structure of vol_basis
        // and possible initial values
    }
    else if (artPrm.eq_mode == CAV)
    {
        ART_numIMG = 1;                    // Normalizing factor = Total no. images
        ptr_vol_out = &vol_basis_out;      // Pointer to output volume
        vol_basis_out = vol_basis;         // Copy the structure of vol_basis
        // and possible initial values
    }
    else
    {
        ART_numIMG = 1;                    // No normalizing factor
        ptr_vol_out = &vol_basis;          // Output volume is the same as
        // input one
    }
    // Now iterate ..........................................................
    ProcessorTimeStamp time0;                    // For measuring the elapsed time
    annotate_processor_time(&time0);
    int images=0;
    double mean_error_2ndblock,pow_residual_imgs;
    bool iv_launched=false;
    for (int it = 0; it < artPrm.no_it; it++)
    {
        // Initialization of some variables
        global_mean_error = 0;
        global_mean_error_1stblock = 0;
        POCS.newIteration();
        VC.newIteration();
        if (rank==-1)
        {
            std::cout << "Running iteration " << it << " with lambda= " << artPrm.lambda(it)<< "\n" << std::endl;
            if (!(artPrm.tell&TELL_SHOW_ERROR))
                init_progress_bar(artPrm.numIMG);
        }

        // For each projection -----------------------------------------------
        for (int act_proj = 0; act_proj < artPrm.numIMG ; act_proj++)
        {
            POCS.newProjection();

            // Select next projection ........................................
            int iact_proj; // Index inside the sorting information for act_proj
            if (artPrm.tell&TELL_MANUAL_ORDER)
            {
                int proj_number;
                std::cout << "Introduce next projection to study: ";
                std::cin  >> proj_number;
                int sym_number=-1;
                if (artPrm.SL.symsNo()!=0)
                {
                    std::cout << "Introduce symmetry to study: ";
                    std::cin  >> sym_number;
                }
                iact_proj=0;
                while (iact_proj<artPrm.numIMG)
                {
                    if (artPrm.IMG_Inf[iact_proj].fn_proj.getNumber()==proj_number &&
                        artPrm.IMG_Inf[iact_proj].sym==sym_number)
                        break;
                    iact_proj++;
                }
            }
            else
            {
                iact_proj = artPrm.ordered_list(act_proj);
            }

            ReconsInfo &imgInfo = artPrm.IMG_Inf[iact_proj];

            read_proj.read(imgInfo.fn_proj, artPrm.apply_shifts, DATA, &imgInfo.row);
            read_proj().setXmippOrigin();
            read_proj.setEulerAngles(imgInfo.rot, imgInfo.tilt, imgInfo.psi);

            // If noisy reconstruction
            if (artPrm.noisy_reconstruction)
            {
                init_random_generator(imgInfo.seed);
                noisy_projection().resize(read_proj());
                noisy_projection().initRandom(0, 1, RND_GAUSSIAN);
                noisy_projection().setXmippOrigin();
                noisy_projection.setEulerAngles(imgInfo.rot,
                                                imgInfo.tilt,
                                                imgInfo.psi);
                if ( it == 0 && imgInfo.sym==-1 )
                {
                    FileName fn_noise;
                    MDRowVec row;
                    fn_noise.compose(read_proj.name().getPrefixNumber(),artPrm.fn_root+"_noise_proj.stk");

                    noisy_projection.write(fn_noise);
                    row.setValue(MDL_IMAGE, fn_noise);
                    row.setValue(MDL_ENABLED, 1);
                    row.setValue(MDL_ANGLE_PSI, read_proj.psi());
                    row.setValue(MDL_ANGLE_ROT, read_proj.rot());
                    row.setValue(MDL_ANGLE_TILT, read_proj.tilt());
                    SF_noise.addRow(row);

                    row.setValue(MDL_IMAGE, read_proj.name());
                    SF_signal.addRow(row);
                }
            }

            //skipping if  tilt greater than max_tilt
            //tilt is in between 0 and 360
            aux_tilt=read_proj.tilt();
            if((aux_tilt > artPrm.max_tilt && aux_tilt < 180.-artPrm.max_tilt) ||
               (aux_tilt > artPrm.max_tilt + 180 && aux_tilt < 360.-artPrm.max_tilt))
            {
                std::cout << "Skipping Proj no: " << iact_proj
                << "tilt=" << read_proj.tilt()  << std::endl;
                continue;
            }

            // Projection extension? .........................................
            if (artPrm.proj_ext!=0)
            {
                read_proj().selfWindow(
                    STARTINGY (read_proj())-artPrm.proj_ext,
                    STARTINGX (read_proj())-artPrm.proj_ext,
                    FINISHINGY(read_proj())+artPrm.proj_ext,
                    FINISHINGX(read_proj())+artPrm.proj_ext);
                noisy_projection().resize(read_proj());
            }

            //Skip if desired
            if (artPrm.tell&TELL_ONLY_SYM)
            {
                if( imgInfo.sym != -1)
                {
                    std::cout << "Skipping Proj no: " << iact_proj
                    << " with symmetry no: " << imgInfo.sym
                    <<  std::endl;
                    continue;
                }
                else
                    std::cout << "NO Skipping Proj no: " << iact_proj
                    << " with symmetry no: " << imgInfo.sym
                    <<  std::endl;
            }

            // For wlsART: use alig_proj for residual image!!
            if (artPrm.WLS)
                alig_proj=artPrm.residual_imgs[iact_proj];

            // Is there a mask ...............................................
            MultidimArray<int> mask;
            const MultidimArray<int> *maskPtr=nullptr;
            if (artPrm.goldmask<1e6 || artPrm.shiftedTomograms)
            {
                mask.resize(read_proj());
                FOR_ALL_ELEMENTS_IN_ARRAY2D(read_proj())
                {
                    mask(i,j)=1;
                    if ((read_proj(i,j)<artPrm.goldmask && artPrm.goldmask<1e6) ||
                        (ABS(read_proj(i,j))<1e-5 && artPrm.shiftedTomograms))
                        mask(i,j)=0;
                }
                maskPtr=&mask;
            }

            // Apply the reconstruction algorithm ............................
            // Notice that the following function is specific for each art type
            singleStep(vol_basis, ptr_vol_out,
                       theo_proj, read_proj, imgInfo.sym , diff_proj,
                       corr_proj, alig_proj,
                       mean_error, ART_numIMG, artPrm.lambda(it),
                       images, imgInfo.fn_ctf, maskPtr,
                       artPrm.refine);

            if (artPrm.WLS)
            {
                artPrm.residual_imgs[iact_proj]=alig_proj;
                global_mean_error_1stblock += diff_proj().sum2()/(XSIZE(diff_proj())*YSIZE(diff_proj()));
            }

            global_mean_error += mean_error;
            if (artPrm.noisy_reconstruction)
            {
                double noise_mean_error;
                singleStep(vol_basis_noisy, &vol_basis_noisy,
                           theo_proj, noisy_projection, imgInfo.sym,
                           diff_proj,  corr_proj, alig_proj,
                           noise_mean_error, ART_numIMG, artPrm.lambda(it),
                           images, imgInfo.fn_ctf,maskPtr,
                           false);
            }

            // Force symmetry in the volume ..................................
            // so far only crystallographic
            if (artPrm.sym_each &&
                act_proj%artPrm.sym_each==0 &&
                (act_proj!=0 || artPrm.sym_each==1))
            {
                //                apply_symmetry(vol_basis,ptr_vol_out, artPrm.grid_type);
                if (artPrm.noisy_reconstruction)
                    //                    apply_symmetry(vol_basis_noisy,&vol_basis_noisy, artPrm.grid_type);
                    ;
            }

            // Apply POCS ....................................................
            POCS.apply(vol_basis,it,images);
            if (artPrm.noisy_reconstruction)
                POCS.apply(vol_basis_noisy,it,images);

            // Variability analysis ..........................................
            if (artPrm.variability_analysis)
                VC.newUpdateVolume(ptr_vol_out,read_proj);

            // Apply anisotropic diffusion ...................................
            if (it>=1 && artPrm.diffusionWeight>-1)
            {
                artPrm.basis.changeToVoxels(vol_basis, &(vol_voxels()),
                                            Zoutput_volume_size, Youtput_volume_size, Xoutput_volume_size);
                Matrix1D<double> alpha;
                alpha=vectorR3(1.0,1.0,artPrm.diffusionWeight);
                double regError=tomographicDiffusion(vol_voxels(),alpha,
                                                     artPrm.lambda(it)/100);
                if (artPrm.tell&TELL_SHOW_ERROR)
                    std::cout << "Regularization error = " << regError << std::endl;
                *artPrm.fh_hist << "Regularization error = " << regError << std::endl;
                artPrm.basis.changeFromVoxels(vol_voxels(),vol_basis,artPrm.grid_type,
                                              artPrm.grid_relative_size, nullptr, nullptr, artPrm.R, artPrm.threads);
            }

            // Apply sparsity constraint .....................................
            if (it>=1 && artPrm.sparseEps>0 &&
                artPrm.parallel_mode!=BasicARTParameters::SIRT &&
                artPrm.parallel_mode!=BasicARTParameters::pSIRT &&
                artPrm.parallel_mode!=BasicARTParameters::pfSIRT)
            {
                artPrm.basis.changeToVoxels(vol_basis, &(vol_voxels()),
                                            Zoutput_volume_size, Youtput_volume_size, Xoutput_volume_size);
                forceDWTSparsity(vol_voxels(),artPrm.sparseEps);
                artPrm.basis.changeFromVoxels(vol_voxels(),vol_basis,artPrm.grid_type,
                                              artPrm.grid_relative_size, nullptr, nullptr, artPrm.R, artPrm.threads);
            }

            // Show results ..................................................
            *artPrm.fh_hist << imgInfo.fn_proj << ", sym="
            << imgInfo.sym << "\t\t" << mean_error;
            if (POCS.apply_POCS)
                *artPrm.fh_hist << "\tPOCS:" << POCS.POCS_mean_error;
            *artPrm.fh_hist << std::endl;
            if (artPrm.tell&TELL_SHOW_ERROR)
            {
                std::cout << imgInfo.fn_proj << ", sym="
                << imgInfo.sym << "\t\t" << mean_error;
                if (POCS.apply_POCS)
                    std::cout        << "\tPOCS:" << POCS.POCS_mean_error;
                std::cout << std::endl;
            }
            else if (act_proj%XMIPP_MAX(1,artPrm.numIMG/60)==0)
                progress_bar(act_proj);

            if (artPrm.tell&TELL_STATS)
            {
                std::cout << "   read      ";
                read_proj().printStats();
                std::cout << "\n   theo      ";
                theo_proj().printStats();
                std::cout << "\n   corr      ";
                corr_proj().printStats();
                std::cout << "\n   alig      ";
                alig_proj().printStats();
                std::cout << "\n   diff      ";
                diff_proj().printStats();
                std::cout << "\n   subvol(0) ";
                (*ptr_vol_out)(0)().printStats();
                std::cout << "\n";
            }

            if (artPrm.tell&TELL_SAVE_AT_EACH_STEP)
            {
                std::cout << "Stats PPPdiff.xmp: ";
                diff_proj().printStats();
                std::cout << std::endl;
                std::cout << "Stats PPPtheo.xmp: ";
                theo_proj().printStats();
                std::cout << std::endl;
                std::cout << "Stats PPPread.xmp: ";
                read_proj().printStats();
                std::cout << std::endl;
                std::cout << "Stats PPPcorr.xmp: ";
                corr_proj().printStats();
                std::cout << std::endl;
                diff_proj.write("PPPdiff.xmp");
                theo_proj.write("PPPtheo.xmp");
                read_proj.write("PPPread.xmp");
                corr_proj.write("PPPcorr.xmp");
                if(act_proj!=0)
                    alig_proj.write("PPPalign.xmp");
                vol_basis.write("PPPbasis.basis");
                artPrm.basis.changeToVoxels(vol_basis, &(vol_voxels()),
                                            Zoutput_volume_size, Youtput_volume_size, Xoutput_volume_size);
                std::cout << "Stats PPPvol.vol : ";
                vol_voxels().printStats();
                std::cout << std::endl;
                vol_voxels.write("PPPvol.vol");

                if (!iv_launched)
                {
                    if (system("xmipp_show -img PPPdiff.xmp PPPtheo.xmp PPPread.xmp PPPcorr.xmp -dont_apply_geo -poll &")==-1)
                        REPORT_ERROR(ERR_UNCLASSIFIED,"Cannot open shell");
                    if (system("xmipp_show -vol PPPvol.vol -poll &")==-1)
                        REPORT_ERROR(ERR_UNCLASSIFIED,"Cannot open shell");
                    iv_launched=true;
                }
                std::cout << "\nHit any key and enter\n";
                char c;
                std::cin >> c;
            }

            // Save intermediate
            if (((artPrm.tell&TELL_SAVE_INTERMIDIATE) | (artPrm.tell&TELL_IV)) &&
                artPrm.save_intermidiate_every!=0 &&
                act_proj%artPrm.save_intermidiate_every==0)
            {
                if (artPrm.tell&TELL_SAVE_INTERMIDIATE)
                    std::cout << "\nSaving intermediate ...\n"
                    << "Converting basis volume to voxels ...\n";
                // Save reconstructed volume
                artPrm.basis.changeToVoxels(vol_basis, &(vol_voxels()),
                                            Zoutput_volume_size, Youtput_volume_size, Xoutput_volume_size);
                if (artPrm.tell&TELL_SAVE_INTERMIDIATE)
                    vol_voxels.write(artPrm.fn_root+"it"+integerToString(it)+"proj"+
                                     integerToString(act_proj,5)+".vol");
                else
                    vol_voxels.write("PPPvol.vol");

                // Launch viewer
                if (!iv_launched)
                {
                    if (system("xmipp_show -i PPPvol.vol --poll &")==-1)
                        REPORT_ERROR(ERR_UNCLASSIFIED,"Cannot open shell");
                    iv_launched=true;
                }
            }

            // Check if algorithm must stop via stop_at
            if (++images==artPrm.stop_at)
                break;
        }

        if (!(artPrm.tell&TELL_SHOW_ERROR))
            progress_bar(artPrm.numIMG);

        // Update residual images for WLS
        if (artPrm.WLS)
        {
            double kappa=artPrm.kappa(it);
            updateResidualVector( artPrm, vol_basis, kappa,
                                  mean_error_2ndblock, pow_residual_imgs);
        }

        // Prepare for next global iteration ---------------------------------
        // Calculate norm and print
        if (rank==-1)
        {
            *artPrm.fh_hist << "Finished - Iteration " << it << std::endl;
            if (artPrm.WLS)
                std::cout        << "   Weighted error: " << global_mean_error
                << " 1st block: "        << global_mean_error_1stblock
                << " 2nd block: "        << mean_error_2ndblock
                << " residual: "         << pow_residual_imgs << std::endl;
            else
                std::cout << "   Global mean squared error: "
                <<  global_mean_error/artPrm.numIMG << std::endl;

            if (artPrm.WLS)
                *artPrm.fh_hist << "   Weighted error: " << global_mean_error
                << " 1st block: "        << global_mean_error_1stblock
                << " 2nd block: "        << mean_error_2ndblock
                << " residual: "         << pow_residual_imgs << std::endl;
            else
                *artPrm.fh_hist << "   Global mean squared error: "
                << global_mean_error/artPrm.numIMG << std::endl;

            if (POCS.apply_POCS)
            {
                std::cout        << "   POCS Global mean squared error: "
                << POCS.POCS_global_mean_error/POCS.POCS_N << std::endl;
                *artPrm.fh_hist << "   POCS Global mean squared error: "
                << POCS.POCS_global_mean_error/POCS.POCS_N << std::endl;
            }
        }

        // Convert volume and write if not last iteration
        // If in SIRT mode move the volume to the reference one
        if ((artPrm.parallel_mode==BasicARTParameters::SIRT
             && !artPrm.variability_analysis)||
            artPrm.parallel_mode==BasicARTParameters::pSIRT ||
            artPrm.parallel_mode==BasicARTParameters::pfSIRT ||
            artPrm.parallel_mode==BasicARTParameters::pSART ||
            artPrm.parallel_mode==BasicARTParameters::pCAV ||
            artPrm.eq_mode==CAV )
        {
            vol_basis=vol_basis_out;
            if (artPrm.sparseEps>0)
            {
                artPrm.basis.changeToVoxels(vol_basis, &(vol_voxels()),
                                            Zoutput_volume_size, Youtput_volume_size, Xoutput_volume_size);
                selfScaleToSize(xmipp_transformation::BSPLINE3,vol_voxels(),
                                (size_t)NEXT_POWER_OF_2(XSIZE(vol_voxels())),
                                (size_t)NEXT_POWER_OF_2(YSIZE(vol_voxels())),
                                (size_t)NEXT_POWER_OF_2(ZSIZE(vol_voxels())));
                Image<double> vol_wavelets, vol_wavelets_abs;
                set_DWT_type(DAUB12);
                DWT(vol_voxels(),vol_wavelets());
                vol_wavelets_abs()=vol_wavelets();
                vol_wavelets_abs().selfABS();
                double *begin=MULTIDIM_ARRAY(vol_wavelets_abs());
                double *end=MULTIDIM_ARRAY(vol_wavelets_abs())+
                            MULTIDIM_SIZE(vol_wavelets_abs());
                std::sort(begin,end);
                double threshold1=DIRECT_MULTIDIM_ELEM(vol_wavelets_abs(),
                                                       (long int)((1-artPrm.sparseEps)*MULTIDIM_SIZE(vol_wavelets_abs())));
                std::cout << "Threshold=" << threshold1 << std::endl;
                vol_wavelets().threshold("abs_below", threshold1, 0.0);
                IDWT(vol_wavelets(),vol_voxels());
                selfScaleToSize(xmipp_transformation::BSPLINE3,vol_voxels(),
                                Xoutput_volume_size,
                                Youtput_volume_size,
                                Zoutput_volume_size);
                artPrm.basis.changeFromVoxels(vol_voxels(),vol_basis,artPrm.grid_type,
                                              artPrm.grid_relative_size, nullptr, nullptr, artPrm.R, artPrm.threads);
            }
        }

        if ( (artPrm.tell & TELL_SAVE_INTERMIDIATE) && it != artPrm.no_it-1)
        {
            if (rank==-1)
                std::cout << "Converting basis volume to voxels ...\n";
            artPrm.basis.changeToVoxels(vol_basis, &(vol_voxels()),
                                        Zoutput_volume_size, Youtput_volume_size, Xoutput_volume_size);
            FileName fn_tmp = artPrm.fn_out.insertBeforeExtension(formatString("_it%06d", it));
            vol_voxels.write(fn_tmp);

            if (artPrm.tell & TELL_SAVE_BASIS)
                vol_basis.write(fn_tmp.removeLastExtension().addExtension("basis"));
        }

        // Check if algorithm must stop via stop_at
        if (images==artPrm.stop_at)
            break;
    }

    // Times on screen
    if (rank==-1)
    {
        std::cout << "\nTime of " << artPrm.no_it << " iterations: \n";
        print_elapsed_time(time0);
    }

    // Finish variability analysis
    VC.finishAnalysis();

    // Save the noisy reconstruction
    if (artPrm.noisy_reconstruction)
    {
        Image<double> vol_voxels_noisy;
        artPrm.basis.changeToVoxels(vol_basis_noisy, &(vol_voxels_noisy()),
                                    Zoutput_volume_size, Youtput_volume_size, Xoutput_volume_size);
        vol_voxels_noisy.write(artPrm.fn_out.insertBeforeExtension("_noise"));
        SF_noise.write(artPrm.fn_root+"_noise_proj.xmd");
        SF_signal.write(artPrm.fn_root+"_signal_proj.xmd");
    }

}

postProcess()

void ARTReconsBase::postProcess ( GridVolume &  vol_basis )

virtual

Finish iterations. For WLS: delete residual images Else: do nothing.

Reimplemented in SinPartARTRecons, and CrystalARTRecons.

Definition at line 604 of file base_art_recons.cpp.

{}

preProcess()

void ARTReconsBase::preProcess ( GridVolume &  vol_basis0, int  level = FULL, int  rank = -1  )

virtual

Produce Plain side information from the Class parameters.

Reimplemented in SinPartARTRecons, and CrystalARTRecons.

Definition at line 40 of file base_art_recons.cpp.

{
    artPrm.produceSideInfo(vol_basis0, level, rank);
}

print()

void ARTReconsBase::print ( std::ostream &  o ) const

virtual

Reimplemented in CrystalARTRecons.

Definition at line 782 of file base_art_recons.cpp.

{}

readParams()

void ARTReconsBase::readParams ( XmippProgram *  program )

virtual

Reimplemented in CrystalARTRecons.

Definition at line 35 of file base_art_recons.cpp.

{
    artPrm.readParams(program);
}

singleStep()

void ARTReconsBase::singleStep ( GridVolume &  vol_in, GridVolume *  vol_out, Projection &  theo_proj, Projection &  read_proj, int  sym_no, Projection &  diff_proj, Projection &  corr_proj, Projection &  alig_proj, double &  mean_error, int  numIMG, double  lambda, int  act_proj, const FileName &  fn_ctf, const MultidimArray< int > *  maskPtr, bool  refine  )

virtual

Run a single step of ART. An ART iteration is compound of as many steps as projections, this function runs a single step of the process. In ART the pointer to the output volume must point to the same vol_in, while in SIRT it should point to a second volume. The read projection must be provided to the algorithm but the rest of projections are output by the routine. The mean error is also an output. numIMG is a normalizing factor to be used in SIRT, if you are running pure ART then this factor should be 1.

The symmetry matrix from which the view is derived must be given in sym_no. In fact, it is not used in this version of ART, but it is needed for the crystal counterpart.

Reimplemented in SinPartARTRecons, and CrystalARTRecons.

Definition at line 601 of file base_art_recons.cpp.

{}

Friends And Related Function Documentation

operator<<

std::ostream& operator<< ( std::ostream &  o, const ARTReconsBase &  artRecons  )

friend

Definition at line 784 of file base_art_recons.cpp.

{
    artRecons.print(o);
    return o;
}

Member Data Documentation

artPrm

BasicARTParameters ARTReconsBase::artPrm

Definition at line 53 of file base_art_recons.h.

---

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

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