phantom.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'
 ***************************************************************************/
/* ------------------------------------------------------------------------- */
/* PHANTOMS                                                                  */
/* ------------------------------------------------------------------------- */

#ifndef _PHANTOM_HH
#define _PHANTOM_HH

#include "core/matrix1d.h"
#include "core/matrix2d.h"
#include "core/xmipp_filename.h"

class Metadata;
template<typename T>
class MultidimArray;
template<typename T>
class Image;
class Projection;
class MetaData;
class MDRow;

/* FEATURE ================================================================= */
class Feature
{
    /* Structure --------------------------------------------------------------- */
public:
    std::string       type;

    char              add_assign;

    double             density;

    Matrix1D<double>   center;

    double             max_distance;

public:
    virtual void prepare() = 0;

    void assign(const Feature &F);

    virtual void rotate(const Matrix2D<double> &E);

    virtual void rotate_center(const Matrix2D<double> &E);

    /* Inside ------------------------------------------------------------------ */
    virtual int point_inside(const Matrix1D<double> &r,
                             Matrix1D<double> &aux) const = 0;

    int point_inside(const Matrix1D<double> &r) const
    {
        Matrix1D<double> aux(3);
        return point_inside(r, aux);
    }

    virtual double density_inside(const Matrix1D<double> &r,
                                  Matrix1D<double> &aux) const = 0;


    int voxel_inside(const Matrix1D<double> &r, Matrix1D<double> &aux1,
                     Matrix1D<double> &aux2) const;

    int voxel_inside(const Matrix1D<double> &r) const
    {
        Matrix1D<double> aux1(3);
        Matrix1D<double> aux2(3);
        return voxel_inside(r, aux1, aux2);
    }

    double voxel_inside_by_normalized_density(const Matrix1D<double> &r, Matrix1D<double> &aux1,
            Matrix1D<double> &aux2) const;


    int intersects_sphere(const Matrix1D<double> &r, double radius,
                          Matrix1D<double> &aux1, Matrix1D<double> &aux2, Matrix1D<double> &aux3)
    const;

    int intersects_sphere(const Matrix1D<double> &r, double radius) const
    {
        Matrix1D<double> aux1(3);
        Matrix1D<double> aux2(3);
        Matrix1D<double> aux3(3);
        return intersects_sphere(r, radius, aux1, aux2, aux3);
    }

    Feature *encircle(double radius = 0) const;

    virtual Feature *scale(double factor) const = 0;

#define ENLARGE_MODE 1
#define SPHERE_MODE  2

    Feature *background(int back_mode, double back_param) const;

    virtual double intersection(const Matrix1D<double> &direction,
                                const Matrix1D<double> &passing_point, Matrix1D<double> &r,
                                Matrix1D<double> &u) const = 0;

    double intersection(const Matrix1D<double> &direction,
                        const Matrix1D<double> &passing_point) const
    {
        Matrix1D<double> r(3);
        Matrix1D<double> u(3);
        return intersection(direction, passing_point, r, u);
    }

    virtual double volume() const = 0;

    void mean_variance_in_plane(Image<double> *V, double z, double &mean, double &var);

    void project_to(Projection &P, const Matrix2D<double> &VP,
                    const Matrix2D<double> &PV) const;

    void corners(const MultidimArray<double> &V, Matrix1D<double> &corner1,
                 Matrix1D<double> &corner2);

#define INTERNAL 0
#define EXTERNAL 1

    void draw_in(MultidimArray<double> &V, int color_mode = INTERNAL, double colour = -1);

    void sketch_in(MultidimArray<double> &V, double colour = 2);

    void shift(double shiftX, double shiftY, double shiftZ);

    void selfApplyGeometry(const Matrix2D<double> &A);

    virtual void feat_printf(FILE *fh) const = 0;
    virtual void feat_printm(MetaData &MD, size_t id) = 0;
    void readCommon(char *line);
    void readCommon(MDRow &row);

 //   void read(char *line);
      void read(MDRow &row);

//    virtual void read_specific(char *line) = 0;
    virtual void read_specific(const std::vector<double> &vector) = 0;

    friend std::ostream& operator << (std::ostream &o, const Feature *F);
};

/* ORIENTED FEATURE ======================================================== */
class Oriented_Feature: public Feature
{
public:
    double             rot;

    double             tilt;

    double             psi;

    Matrix2D<double>   euler;

    Matrix2D<double>   eulert;
public:
    void prepare_Euler();

    void assign(const Oriented_Feature & OF);

    virtual void rotate(const Matrix2D<double> &E);
};

/* SPHERES ================================================================= */
class Sphere: public Feature
{
public:
    double             radius;
public:
    void prepare();

    void assign(const Sphere &F);

    int point_inside(const Matrix1D<double> &r, Matrix1D<double> &aux) const;

    double density_inside(const Matrix1D<double> &r, Matrix1D<double> &aux) const
    {
        return 1.;
    }

    Feature *scale(double factor) const;

    __attribute__ ((deprecated)) void scale(double factor, Feature **_f) const;

    double intersection(const Matrix1D<double> &direction,
                        const Matrix1D<double> &passing_point,
                        Matrix1D<double> &r, Matrix1D<double> &u) const;

    double volume() const
    {
        return 4 / 3*PI*radius*radius*radius;
    }

    void read_specific(char *line);
    void read_specific(const std::vector<double> &vect);

    void feat_printf(FILE *fh) const;
    void feat_printm(MetaData &MD, size_t id);
    friend std::ostream& operator << (std::ostream &o, const Sphere &f);

    void init_rnd(
        double minradius,   double maxradius,
        double minden = 1.0,  double maxden = 1.0,
        double minx0 = 0,     double maxx0 = 0,
        double miny0 = 0,     double maxy0 = 0,
        double minz0 = 0,     double maxz0 = 0);
};

/* BLOB ================================================================= */
class Blob: public Feature
{
public:

    /* I do not use the structure "blobtype" because common code with sphere becomes
    more difficult*/
    double             radius;
    double             alpha;
    int         m;
public:
    void prepare();

    void assign(const Blob &F);

    int point_inside(const Matrix1D<double> &r, Matrix1D<double> &aux) const;

    double density_inside(const Matrix1D<double> &r, Matrix1D<double> &aux) const;

    Feature *scale(double factor) const;

    __attribute__ ((deprecated)) void scale(double factor, Feature **_f) const;

    double intersection(const Matrix1D<double> &direction,
                        const Matrix1D<double> &passing_point,
                        Matrix1D<double> &r, Matrix1D<double> &u) const;

    double volume() const;

    void read_specific(char *line);
    void read_specific(const std::vector<double> &vect);

    void feat_printf(FILE *fh) const;
    void feat_printm(MetaData &MD, size_t id);
    friend std::ostream& operator << (std::ostream &o, const Blob &f);

    void init_rnd(
        double minradius,   double maxradius,
        double minalpha,    double maxalpha,
        double minorder,    double maxorder,
        double minden = 1.0,  double maxden = 1.0,
        double minx0 = 0,     double maxx0 = 0,
        double miny0 = 0,     double maxy0 = 0,
        double minz0 = 0,     double maxz0 = 0);
};

/* GAUSSIAN ================================================================ */
class Gaussian: public Feature
{
public:

    double             sigma;
public:
    void prepare();

    void assign(const Gaussian &F);

    int point_inside(const Matrix1D<double> &r, Matrix1D<double> &aux) const;

    double density_inside(const Matrix1D<double> &r, Matrix1D<double> &aux) const;

    Feature *scale(double factor) const;

    __attribute__ ((deprecated)) void scale(double factor, Feature **_f) const;

    double intersection(const Matrix1D<double> &direction,
                        const Matrix1D<double> &passing_point,
                        Matrix1D<double> &r, Matrix1D<double> &u) const;

    double volume() const
    {
        return 1;
    }

    void read_specific(char *line);
    void read_specific(const std::vector<double> &vect);
    void feat_printf(FILE *fh) const;
    void feat_printm(MetaData &MD, size_t id);
    friend std::ostream& operator << (std::ostream &o, const Gaussian &f);

    void init_rnd(
        double minsigma,   double maxsigma,
        double minden = 1.0,  double maxden = 1.0,
        double minx0 = 0,     double maxx0 = 0,
        double miny0 = 0,     double maxy0 = 0,
        double minz0 = 0,     double maxz0 = 0);
};

/* CYLINDERS =============================================================== */
class Cylinder: public Oriented_Feature
{
public:
    double             xradius;

    double             yradius;

    double             height;
public:
    void prepare();

    void assign(const Cylinder &F);

    int point_inside(const Matrix1D<double> &r, Matrix1D<double> &aux) const;

    double density_inside(const Matrix1D<double> &r, Matrix1D<double> &aux) const
    {
        return 1.;
    }

    Feature *scale(double factor) const;

    __attribute__ ((deprecated)) void scale(double factor, Feature **_f) const;

    double intersection(const Matrix1D<double> &direction,
                        const Matrix1D<double> &passing_point,
                        Matrix1D<double> &r, Matrix1D<double> &u) const;

    double volume() const
    {
        return 4 / 3*PI*xradius*yradius*height;
    }

    void read_specific(char *line);
    void read_specific(const std::vector<double> &vect);
    void  feat_printf(FILE *fh) const;
    void feat_printm(MetaData &MD, size_t id);
    friend std::ostream& operator << (std::ostream &o, const Cylinder &f);

    void  init_rnd(
        double minxradius,  double maxxradius,
        double minyradius,  double maxyradius,
        double minheight,   double maxheight,
        double minden = 1.0,  double maxden = 1.0,
        double minx0 = 0,     double maxx0 = 0,
        double miny0 = 0,     double maxy0 = 0,
        double minz0 = 0,     double maxz0 = 0,
        double minrot = 0,    double maxrot = 360,
        double mintilt = 0,   double maxtilt = 180,
        double minpsi = 0,    double maxpsi = 360);
};

/* DOUBLE CYLINDERS ======================================================== */
class DCylinder: public Oriented_Feature
{
public:
    double             radius;

    double             height;

    double             separation;
public:
    void prepare();

    void assign(const DCylinder &F);

    int point_inside(const Matrix1D<double> &r, Matrix1D<double> &aux) const;

    double density_inside(const Matrix1D<double> &r, Matrix1D<double> &aux) const
    {
        return 1.;
    }

    Feature *scale(double factor) const;

    __attribute__ ((deprecated)) void scale(double factor, Feature **_f) const;

    double intersection(const Matrix1D<double> &direction,
                        const Matrix1D<double> &passing_point,
                        Matrix1D<double> &r, Matrix1D<double> &u) const;

    double volume() const
    {
        return 2* 4 / 3*PI*radius*radius*height;
    }

    void read_specific(char *line);
    void read_specific(const std::vector<double> &vect);
    void feat_printf(FILE *fh) const;
    void feat_printm(MetaData &MD, size_t id);
    friend std::ostream& operator << (std::ostream &o, const DCylinder &f);

    void init_rnd(
        double minradius,   double maxradius,
        double minheight,   double maxheight,
        double minsep,      double maxsep,
        double minden = 1.0,  double maxden = 1.0,
        double minx0 = 0,     double maxx0 = 0,
        double miny0 = 0,     double maxy0 = 0,
        double minz0 = 0,     double maxz0 = 0,
        double minrot = 0,    double maxrot = 360,
        double mintilt = 0,   double maxtilt = 180,
        double minpsi = 0,    double maxpsi = 360);
};

/* CUBE ==================================================================== */
class Cube: public Oriented_Feature
{
public:
    double             xdim;

    double             ydim;

    double             zdim;
public:
    void prepare();

    void assign(const Cube &F);

    int point_inside(const Matrix1D<double> &r, Matrix1D<double> &aux) const;

    double density_inside(const Matrix1D<double> &r, Matrix1D<double> &aux) const
    {
        return 1.;
    }

    Feature *scale(double factor) const;

    __attribute__ ((deprecated)) void scale(double factor, Feature **_f) const;

    double intersection(const Matrix1D<double> &direction,
                        const Matrix1D<double> &passing_point,
                        Matrix1D<double> &r, Matrix1D<double> &u) const;

    double volume() const
    {
        return xdim*ydim*zdim;
    }

    void read_specific(char *line);
    void read_specific(const std::vector<double> &vect);
    void feat_printf(FILE *fh) const;
    void feat_printm(MetaData &MD, size_t id);
    friend std::ostream& operator << (std::ostream &o, const Cube &f);

    void  init_rnd(
        double minXdim,   double maxXdim,
        double minYdim = 0, double maxYdim = 0,
        double minZdim = 0, double maxZdim = 0,
        double minden = 1.0,  double maxden = 1.0,
        double minx0 = 0,   double maxx0 = 0,
        double miny0 = 0,   double maxy0 = 0,
        double minz0 = 0,   double maxz0 = 0,
        double minrot = 0,  double maxrot = 360,
        double mintilt = 0, double maxtilt = 180,
        double minpsi = 0,  double maxpsi = 360);
};

/* ELLIPSOID =============================================================== */
class Ellipsoid: public Oriented_Feature
{
public:
    double             xradius;

    double             yradius;

    double             zradius;
public:
    void prepare();

    void assign(const Ellipsoid &F);

    int point_inside(const Matrix1D<double> &r, Matrix1D<double> &aux) const;

    double density_inside(const Matrix1D<double> &r, Matrix1D<double> &aux) const
    {
        return 1.;
    }

    Feature *scale(double factor) const;

    __attribute__ ((deprecated)) void scale(double factor, Feature **_f) const;

    double intersection(const Matrix1D<double> &direction,
                        const Matrix1D<double> &passing_point,
                        Matrix1D<double> &r, Matrix1D<double> &u) const;

    double volume() const
    {
        return 4 / 3*PI*xradius*yradius*zradius;
    }

    void read_specific(char *line);
    void read_specific(const std::vector<double> &vect);
    void feat_printf(FILE *fh) const;
    void feat_printm(MetaData &MD, size_t id);
    friend std::ostream& operator << (std::ostream &o, const Ellipsoid &f);

    void init_rnd(
        double minXradius, double maxXradius,
        double minYradius, double maxYradius,
        double minZradius, double maxZradius,
        double minden = 1.0, double maxden = 1.0,
        double minx0 = 0,    double maxx0 = 0,
        double miny0 = 0,    double maxy0 = 0,
        double minz0 = 0,    double maxz0 = 0,
        double minrot = 0,   double maxrot = 360,
        double mintilt = 0,  double maxtilt = 180,
        double minpsi = 0,   double maxpsi = 360);
};

/* CONE ==================================================================== */
class Cone: public Oriented_Feature
{
public:
    double             radius;

    double             height;

public:
    void prepare();

    void assign(const Cone &F);

    int point_inside(const Matrix1D<double> &r, Matrix1D<double> &aux) const;

    double density_inside(const Matrix1D<double> &r, Matrix1D<double> &aux) const
    {
        return 1.;
    }

    Feature *scale(double factor) const;

    __attribute__ ((deprecated)) void scale(double factor, Feature **_f) const;

    double intersection(const Matrix1D<double> &direction,
                        const Matrix1D<double> &passing_point,
                        Matrix1D<double> &r, Matrix1D<double> &u) const;

    double volume() const
    {
        return 1 / 3*PI*radius*radius*height;
    }

    void read_specific(char *line);
    void read_specific(const std::vector<double> &vect);
    void feat_printf(FILE *fh) const;
    void feat_printm(MetaData &MD, size_t id);
    friend std::ostream& operator << (std::ostream &o, const Cone &f);

    void init_rnd(
        double minradius,   double maxradius,
        double minheight,   double maxheight,
        double minden = 1.0,  double maxden = 1.0,
        double minx0 = 0,     double maxx0 = 0,
        double miny0 = 0,     double maxy0 = 0,
        double minz0 = 0,     double maxz0 = 0,
        double minrot = 0,    double maxrot = 360,
        double mintilt = 0,   double maxtilt = 180,
        double minpsi = 0,    double maxpsi = 360);
};

/* PHANTOM ================================================================= */
class Phantom
{
public:
    FileName       fn;

    int            xdim;

    int            ydim;

    int            zdim;

    double          Background_Density;

    double          current_scale;

    double          param_file_scale;

    double          phantom_scale;

    std::vector<Feature*> VF;
public:
    Phantom();

    Phantom(const FileName &fn_phantom)
    {
        read(fn_phantom);
    }

    Phantom(const Phantom &other);

    ~Phantom()
    {
        clear();
    }

    Phantom & operator = (const Phantom &P);

    void clear();

    int FeatNo()
    {
        return VF.size();
    }

    Feature * operator()(int i)
    {
        return VF[i-1];
    }

    const Feature * operator()(int i) const
    {
        return VF[i-1];
    }

    void add(Feature *f)
    {
        VF.push_back(f);
    }

    void assign(const Phantom &P);

    void prepare();

    double max_distance() const;

    double volume() const;

    void read(const FileName &fn_phantom, bool apply_scale = true);

    friend std::ostream& operator << (std::ostream &o, const Phantom &f);

    void write(const FileName &fn_phantom = "");

    int voxel_inside_any_feat(const Matrix1D<double> &r,
                              Matrix1D<double> &aux1, Matrix1D<double> &aux2) const;

    int voxel_inside_any_feat(const Matrix1D<double> &r) const
    {
        Matrix1D<double> aux1(3);
        Matrix1D<double> aux2(3);
        return voxel_inside_any_feat(r, aux1, aux2);
    }

    int any_feature_intersects_sphere(const Matrix1D<double> &r, double radius,
                                      Matrix1D<double> &aux1, Matrix1D<double> &aux2, Matrix1D<double> &aux3)
    const;

    int any_feature_intersects_sphere(const Matrix1D<double> &r,
                                      double radius) const
    {
        Matrix1D<double> aux1(3);
        Matrix1D<double> aux2(3);
        Matrix1D<double> aux3(3);
        return any_feature_intersects_sphere(r, radius, aux1, aux2, aux3);
    }

    void draw_in(MultidimArray<double> &V);

    void label(MultidimArray<double> &V);

#define DONT_CLEAN 0
#define CLEAN      1

    void sketch_in(MultidimArray<double> &V, int clean = DONT_CLEAN, double colour = 2);

    void shift(double shiftX, double shiftY, double shiftZ);

    void rotate(const Matrix2D<double> &E);

    void selfApplyGeometry(const Matrix2D<double> &A, int inv);

    void project_to(Projection &P, int Ydim, int Xdim,
                    double rot, double tilt, double psi, const Matrix2D<double> *A = nullptr) const;

    void project_to(Projection &P,
                    double rot, double tilt, double psi, const Matrix2D<double> *A = nullptr) const;

    void project_to(Projection &P, const Matrix2D<double> &VP, double    disappearing_th = 1.0) const;

#define POS_NEG 1
#define NEG_POS 2

    void surface(double z0, double radius, int direction, Image<double> *P)
    const;
};
#endif