projection.h

Go to the documentation of this file.

/***************************************************************************
 *
 * Authors:     Carlos Oscar S. Sorzano (coss@cnb.csic.es)
 *
 * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
 * 02111-1307  USA
 *
 *  All comments concerning this program package may be sent to the
 *  e-mail address 'xmipp@cnb.csic.es'
 ***************************************************************************/

#ifndef CORE_PROJECTION_H
#define CORE_PROJECTION_H

#include "core/xmipp_filename.h"
#include "core/xmipp_threads.h"
#include "core/matrix1d.h"
#include "grids.h"

class Projection;
template<typename T>
class Image;
template<typename T>
class Matrix2D;
template<typename T>
class MultidimArray;
class SimpleGrid;
class Basis;
class XmippProgram;


// These two structures are needed when projecting and backprojecting using
// threads. They make mutual exclusion and synchronization possible.
extern barrier_t project_barrier;
extern pthread_mutex_t project_mutex;

/* Projection parameters for tomography --------------------------- */
class ParametersProjectionTomography
{
public:
    FileName fnPhantom;
    FileName  fnRoot;
    FileName  fnOut;
    bool singleProjection;
    int      starting;
    std::string   fn_projection_extension;

    int      proj_Xdim;
    int      proj_Ydim;

    bool only_create_angles;
    // Show angles calculation in std::out
    bool show_angles;

    double axisRot;
    double axisTilt;
    Matrix1D<double> raxis;
    double tilt0;
    double tiltF;
    double tiltStep;

    double    Npixel_avg;
    double    Npixel_dev;

    double    Ncenter_avg;
    double    Ncenter_dev;

    double    Nangle_avg;
    double    Nangle_dev;

public:

    ParametersProjectionTomography();

    void read(const FileName &fn_proj_param);
    void defineParams(XmippProgram * program);
    void readParams(XmippProgram * program);

    void calculateProjectionAngles(Projection &P, double angle, double inplaneRot,
                                   const Matrix1D<double> &sinplane);
};

typedef struct
{
    int thread_id;
    int threads_count;
    Image<double> * vol;
    const SimpleGrid * grid;
    const Basis * basis;
    Projection * global_proj;
    Projection * global_norm_proj;
    int FORW;
    int eq_mode;
    const Image<int> *VNeq;
    Matrix2D<double> *M;
    const MultidimArray<int> *mask;
    double ray_length;
    double rot;
    double tilt;
    double psi;
    bool destroy;
}
project_thread_params;

extern project_thread_params * project_threads;

template <class T>
void project_SimpleGrid(Image<T> *vol, const SimpleGrid *grid,
                        const Basis *basis,
                        Projection *proj, Projection *norm_proj, int FORW, int eq_mode,
                        const Image<int> *VNeq, Matrix2D<double> *M,
                        const MultidimArray<int> *mask=nullptr,
                        double ray_length = -1.0,
                        int thread_id = -1, int num_threads = 1);

/*---------------------------------------------------------------------------*/
/* PROJECTION GENERATION                                                     */
/*---------------------------------------------------------------------------*/
// Projecting functions ====================================================
constexpr int FORWARD =  1;
constexpr int BACKWARD = 0;

constexpr int ARTK =     1;
constexpr int CAVK =     2;
constexpr int COUNT_EQ = 3;
constexpr int CAV =      4;
constexpr int CAVARTK =  5;

void projectVolume(MultidimArray<double> &V, Projection &P, int Ydim, int Xdim,
                   double rot, double tilt, double psi,
                   const Matrix1D<double> *roffset=nullptr);

void projectVolumeOffCentered(MultidimArray<double> &V, Projection &P,
                              int Ydim, int Xdim);

void singleWBP(MultidimArray<double> &V, Projection &P);

void count_eqs_in_projection(GridVolumeT<int> &GVNeq,
                             const Basis &basis, Projection &read_proj);

void project_Crystal_Volume(GridVolume &vol, const Basis &basis,
                            Projection &proj, Projection &norm_proj,
                            int Ydim, int Xdim,
                            double rot, double tilt, double psi, const Matrix1D<double> &shift,
                            const Matrix1D<double> &aint, const Matrix1D<double> &bint,
                            const Matrix2D<double> &D, const Matrix2D<double> &Dinv,
                            const MultidimArray<int> &mask, int FORW, int eq_mode = ARTK);

// Implementations =========================================================
// Some aliases
#define x0   STARTINGX(IMGMATRIX(*proj))
#define xF   FINISHINGX(IMGMATRIX(*proj))
#define y0   STARTINGY(IMGMATRIX(*proj))
#define yF   FINISHINGY(IMGMATRIX(*proj))
#define xDim XSIZE(IMGMATRIX(*proj))
#define yDim YSIZE(IMGMATRIX(*proj))

// Projections from single particles #######################################
// Projections from basis volumes ==========================================

/* Projection of a simple volume ------------------------------------------- */
// The projection is not precleaned (set to 0) before projecting and its
// angles are supposed to be already written (and all Euler matrices
// precalculated
// The projection plane is supposed to pass through the Universal coordinate
// origin

/* Time measures ...........................................................
   This function has been optimized in time, several approaches have been
   tried and here are the improvements in time
   Case: Base (unit time measure)
         The points within the footprint are calculated inside the
         inner loop instead of computing the foot coordinate for the
         first corner and knowing the sampling rate in the oversampled
         image go on the next points. The image access was done directly
         (physical access).

   Notation simplification:
   Case: VOLVOXEL(V,k,i,j) ----> V(k,i,j):   + 38% (Inacceptable)
   Case: DIRECT_IMGPIXEL   ----> IMGPIXEL:   +  5% (Acceptable)

   Algorithmic changes:
   Case: project each basis position ---->   +325% (Inacceptable)
         get rid of all beginZ, beginY, prjX, ... and project each
         basis position onto the projection plane
   Case: footprint coordinates outside ----> - 33% (Great!!)
         the footprint coordinates are computed outside the inner
         loop, but you have to use the sampling rate instead to
         move over the oversampled image.
*/

const int ART_PIXEL_SUBSAMPLING = 2;

template <class T>
void *project_SimpleGridThread( void * params );

template <class T>
void project_SimpleGrid(Image<T> *vol, const SimpleGrid *grid,
                        const Basis *basis,
                        Projection *proj, Projection *norm_proj, int FORW, int eq_mode,
                        const Image<int> *VNeq, Matrix2D<double> *M,
                        const MultidimArray<int> *mask,
                        double ray_length,
                        int thread_id, int numthreads);

/* Project a Grid Volume --------------------------------------------------- */
template <class T>
void project_GridVolume(
    GridVolumeT<T> &vol,                  // Volume
    const Basis &basis,                   // Basis
    Projection       &proj,               // Projection
    Projection       &norm_proj,          // Projection of a unitary volume
    int              Ydim,                // Dimensions of the projection
    int              Xdim,
    double rot, double tilt, double psi,  // Euler angles
    int              FORW,                // 1 if we are projecting a volume
    //   norm_proj is calculated
    // 0 if we are backprojecting
    //   norm_proj must be valid
    int              eq_mode,             // ARTK, CAVARTK, CAVK or CAV
    GridVolumeT<int> *GVNeq,              // Number of equations per blob
    Matrix2D<double> *M,                  // System matrix
    const MultidimArray<int> *mask,            // mask(i,j)=0 => do not update this pixel
    double            ray_length,   // Ray length of the projection
    int               threads);

#undef x0
#undef xF
#undef y0
#undef yF
#undef xDim
#undef yDim

#endif