Phantom Class Reference

#include <phantom.h>

Collaboration diagram for Phantom:

Public Member Functions
	Phantom ()
	Phantom (const FileName &fn_phantom)
	Phantom (const Phantom &other)
	~Phantom ()
Phantom &	operator= (const Phantom &P)
	Assignment. More...
void	clear ()
int	FeatNo ()
Feature *	operator() (int i)
const Feature *	operator() (int i) const
void	add (Feature *f)
void	assign (const Phantom &P)
void	prepare ()
	Prepare for work. More...
double	max_distance () const
	Return the maximum distance of any feature to the volume center. More...
double	volume () const
	Return the volume of all the features. More...
void	read (const FileName &fn_phantom, bool apply_scale=true)
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
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
void	draw_in (MultidimArray< double > &V)
void	label (MultidimArray< double > &V)
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
void	surface (double z0, double radius, int direction, Image< double > *P) const

Public Attributes
FileName	fn
	Filename. More...
int	xdim
	Final volume X dimension. More...
int	ydim
	Final volume Y dimension. More...
int	zdim
	Final volume Z dimension. More...
double	Background_Density
	Final volume background density. More...
double	current_scale
	Has been the scale applied? More...
double	param_file_scale
	Param file scale. More...
double	phantom_scale
	Param file scale. More...
std::vector< Feature * >	VF
	List with the features. More...

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

Detailed Description

Phantom class. The phantom class is simply a list (STL vector) of features plus some information about the size of the final volume to generate and its background density. This is the class that will interact with the reconstruction programs as the features classes themselves haven't got enough information to generate the final volume. The file format to generate the phantom is described in the previous page (Phantoms).

This class is thought to be filled from a file, and doesn't give many facilities to update it from program. This is something to do.

Here goes an example of how to manage loops in the phantom class,

// Show all features
for (int i=1; i<=P.FeatNo(); i++) std::cout << P(i);

Definition at line 1352 of file phantom.h.

Constructor & Destructor Documentation

Phantom() [1/3]

Phantom::Phantom

(

)

Empty constructor. The empty phantom is 0x0x0, background density=0, no feature is inside and no name.

Definition at line 1967 of file phantom.cpp.

{
    xdim = ydim = zdim = 0;
    Background_Density = 0;
    fn = "";
    current_scale = 1;
    phantom_scale = 1.;
}

Phantom() [2/3]

Phantom::Phantom ( const FileName &  fn_phantom )

inline

Construct from a phantom file. Construct the phantom according to the specifications of the given file. The file must accomplish the structure given in Phantoms.

Definition at line 1390 of file phantom.h.

    {
        read(fn_phantom);
    }

Phantom() [3/3]

Phantom::Phantom

(

const Phantom &

other

)

Copy constructor. The empty phantom is 0x0x0, background density=0, no feature is inside and no name.

Definition at line 1976 of file phantom.cpp.

{
    *this = other;
}

~Phantom()

Phantom::~Phantom ( )

inline

Destructor. All features are freed.

Definition at line 1402 of file phantom.h.

    {
        clear();
    }

Member Function Documentation

add()

void Phantom::add ( Feature *  f )

inline

Add a feature to the phantom. This function allows you to add new features at the end of the phantom. \ Ex: Sphere S; S.radius(); S.prepare(); Phantom.add(&S);

Definition at line 1438 of file phantom.h.

    {
        VF.push_back(f);
    }

any_feature_intersects_sphere() [1/2]

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

Speeded up sphere intersecting any feature. This function returns the first feature in the list intersecting a sphere with center r in R3 and the given radius. In none, 0 is returned. This speeded up function needs two vectors with dimension 3 externally resized.

Definition at line 2416 of file phantom.cpp.

{
    bool intersects;
    for (size_t i = 0; i < VF.size(); i++)
    {
        intersects = VF[i]->intersects_sphere(r, radius, aux1, aux2, aux3);
        if (intersects)
            return i + 1;
    }
    return 0;
}

any_feature_intersects_sphere() [2/2]

int Phantom::any_feature_intersects_sphere ( const Matrix1D< double > &  r, double  radius  ) const

inline

Sphere intersecting feature. This function is based in the previous one. It makes the same but you needn't supply the auxiliar vectors.

Definition at line 1515 of file phantom.h.

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

assign()

void Phantom::assign

(

const Phantom &

P

)

Another function for assignment.

clear()

void Phantom::clear

(

)

Clear the phantom. Force the phantom to be empty. All features are freed.

Definition at line 1981 of file phantom.cpp.

{
    xdim = ydim = zdim = 0;
    Background_Density = 0;
    fn = "";
    for (size_t i = 0; i < VF.size(); i++)
        delete VF[i];
    VF.clear();
}

draw_in()

void Phantom::draw_in

(

MultidimArray< double > &

V

)

Draw the phantom in the volume. The volume is cleaned, resized to the phantom size and its origin is set at the center of the volume. Then every feature is drawn into the volume

Definition at line 2432 of file phantom.cpp.

{
    V.resize(zdim, ydim, xdim);
    V.setXmippOrigin();
    V.initConstant(Background_Density);
    for (size_t i = 0; i < VF.size(); i++)
        VF[i]->draw_in(V);
}

FeatNo()

int Phantom::FeatNo ( )

inline

Returns the number of features in the list.

Definition at line 1415 of file phantom.h.

    {
        return VF.size();
    }

label()

void Phantom::label

(

MultidimArray< double > &

V

)

Label a volume after the phantom. The volume is cleaned, resized to the phantom size and its origin is set at the center of the volume. Then every feature is drawn into the volume using different densities. Background has got density 0, the border for first feature density -1, the first feature has got density 1, the second 2 and its border -2, ...

Definition at line 2443 of file phantom.cpp.

{
    Matrix1D<double> r(3);
    Matrix1D<double> aux1(3);
    Matrix1D<double> aux2(3);
    V.resize(zdim, ydim, xdim);
    V.setXmippOrigin();
    FOR_ALL_ELEMENTS_IN_ARRAY3D(V)
    {
        ZZ(r) = k;
        YY(r) = i;
        XX(r) = j;
        int sel_feat = voxel_inside_any_feat(r, aux1, aux2);
        // If it is not in the background, check that it is completely
        // inside the feature, if not set it to border.
        if (sel_feat != 0)
            if (VF[sel_feat-1]->voxel_inside(r, aux1, aux2) != 8)
                sel_feat = -sel_feat;
        A3D_ELEM(V, k, i, j) = sel_feat;
    }
}

max_distance()

double Phantom::max_distance

(

)

const

Return the maximum distance of any feature to the volume center.

Definition at line 2070 of file phantom.cpp.

{
    double retval = 0;
    for (size_t i = 0; i < VF.size(); i++)
        retval = XMIPP_MAX(retval, VF[i]->max_distance + VF[i]->center.module());
    return retval;
}

operator()() [1/2]

Feature* Phantom::operator() ( int  i )

inline

Access to a feature pointer. You can address each one of the features using this operator, remember that features are numbered as 1, 2, ... FeatNo()

Definition at line 1423 of file phantom.h.

    {
        return VF[i-1];
    }

operator()() [2/2]

const Feature* Phantom::operator() ( int  i ) const

inline

Constant Access to a feature pointer. Same as the previous one but with constant access.

Definition at line 1430 of file phantom.h.

    {
        return VF[i-1];
    }

operator=()

Phantom & Phantom::operator=

(

const Phantom &

P

)

Assignment.

Definition at line 1991 of file phantom.cpp.

{
    if (&P == this)
        return *this;
    clear();
    fn = P.fn;
    xdim = P.xdim;
    ydim = P.ydim;
    zdim = P.zdim;
    phantom_scale = P.phantom_scale;
    Background_Density = P.Background_Density;
    Sphere     *sph;
    Blob       *blo;
    Gaussian   *gau;
    Cylinder   *cyl;
    DCylinder  *dcy;
    Cube       *cub;
    Ellipsoid  *ell;
    Cone       *con;
    for (size_t i = 0; i < P.VF.size(); i++)
        if (P.VF[i]->type == "sph")
        {
            sph = new Sphere;
            *sph = *((Sphere *)    P.VF[i]);
            add(sph);
        }
        else if (P.VF[i]->type == "blo")
        {
            blo = new Blob;
            *blo = *((Blob *)      P.VF[i]);
            add(blo);
        }
        else if (P.VF[i]->type == "gau")
        {
            gau = new Gaussian;
            *gau = *((Gaussian *)  P.VF[i]);
            add(gau);
        }
        else if (P.VF[i]->type == "cyl")
        {
            cyl = new Cylinder;
            *cyl = *((Cylinder *)  P.VF[i]);
            add(cyl);
        }
        else if (P.VF[i]->type == "dcy")
        {
            dcy = new DCylinder;
            *dcy = *((DCylinder *) P.VF[i]);
            add(dcy);
        }
        else if (P.VF[i]->type == "cub")
        {
            cub = new Cube;
            *cub = *((Cube *)      P.VF[i]);
            add(cub);
        }
        else if (P.VF[i]->type == "ell")
        {
            ell = new Ellipsoid;
            *ell = *((Ellipsoid *) P.VF[i]);
            add(ell);
        }
        else if (P.VF[i]->type == "con")
        {
            con = new Cone;
            *con = *((Cone      *) P.VF[i]);
            add(con);
        }
    return *this;
}

prepare()

void Phantom::prepare

(

)

Prepare for work.

Definition at line 2063 of file phantom.cpp.

{
    for (size_t i = 0; i < VF.size(); i++)
        VF[i]->prepare();
}

project_to() [1/3]

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

Project phantom from a direction. The direction is specified by the 3 Euler angles (as usual, rot=1st, tilt=2nd, psi=3rd), the projection size is (Ydim x Xdim) given in the function call, the projection logical origin is set at the physical center of the image. The projection is cleaned and then the projection of the phantom is computed. A matrix A (3x3) can be supplied in order to apply a deformation in the projection plane. A must be such that

deformed position=A*undeformed position

Definition at line 2511 of file phantom.cpp.

{
#ifdef DEBUG
    std::cout << "Ydim=" << Ydim << " Xdim=" << Xdim << std::endl
    << "rot=" << rot << " tilt=" << tilt << " psi=" << psi << std::endl
    << "A\n" << A;
#endif
    // Initialise projection
    P().initZeros(Ydim, Xdim);
    P().setXmippOrigin();
    P.setAngles(rot, tilt, psi);

    // Compute volume to Projection matrix
    Matrix2D<double> VP = P.euler;
    if (A != nullptr)
        VP = (*A) * VP;
    Matrix2D<double> PV = VP.inv();
    // Project all features
    for (size_t i = 0; i < VF.size(); i++)
        VF[i]->project_to(P, VP, PV);
}

project_to() [2/3]

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

Project phantom from a direction. The same as before but this time the projection is supposed to be already resized and with the right center. The phantom projection is added to the already drawn projection.

Definition at line 2535 of file phantom.cpp.

{
    P.setAngles(rot, tilt, psi);

    // Compute volume to Projection matrix
    Matrix2D<double> VP = P.euler;
    if (A != nullptr)
        VP = (*A) * VP;
    Matrix2D<double> PV = VP.inv();

    // Project all features
    for (size_t i = 0; i < VF.size(); i++)
        VF[i]->project_to(P, VP, PV);
}

project_to() [3/3]

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

Project phantom using a conversion matrix. The same as before but this time the projection is supposed to be already resized and with the right center. The phantom projection is added to the already drawn projeciton. Euler angles of the projection are not set and the phantom volume is projected to the projection plane making use of the matrix provided (3x3) which is the needed transformation from the volume to the projection plane. Inside the projection plane only the first two coordinates are valid. If disappearing_th < 1 then the ith phantom feature is skiped with a provability = 1-disappearing_th

Definition at line 2551 of file phantom.cpp.

{
    Matrix2D<double> PV = VP.inv();

    // Project all features
    for (size_t i = 0; i < VF.size(); i++)
    {
        if (rnd_unif(0, 1) < disappearing_th)
            VF[i]->project_to(P, VP, PV);
    }
}

read()

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

Read a phantom file. The file must accomplish the structure given in Phantoms .

If you don't apply the scale then all spatial coordinates are expressed in the given scale units. I.e., if the scale is 0.25 that means that every voxel in the voxel phantom represent 4 length units (usually Angstroms). If the scale is applied, then coordinates are expressed in voxel edge units. Otherwise, coordinates are expressed in length units (usually Angstroms).

We recommend apply_scale=false only if you plan to modify the description file without produstd::cing projections, voxel phantoms, ...

Definition at line 2088 of file phantom.cpp.

{

    FILE *fh_phantom;
    char line[256];
    int Global_Feature_Read = 0; // Indicates if the line with volume dimensions
    // has been already read
    int        stat;
    Sphere     *sph;
    Blob       *blo;
    Gaussian   *gau;
    Cylinder   *cyl;
    DCylinder  *dcy;
    Cube       *cub;
    Ellipsoid  *ell;
    Cone       *con;
    Feature    *feat; 
    Feature    *scaled_feat;
    std::string feat_type;
    double     scale = 1.;          // The scale factor is not stored
    char       straux[6];

    // Clear actual phantom
    clear();

    if (fn_phantom.isMetaData())
    {
        MetaDataVec MD1;  //MetaData for the first block (phantom parameters)
        MetaDataVec MD2; //MetaData for the second block (phantom parameters)
        std::vector <double> TempVec; // A temporary vector for reading vector data
        size_t objId;

        // Assign different blocks to different MetaDatas
        MD1.read((std::string)"block1@"+fn_phantom.c_str());
        MD2.read((std::string)"block2@"+fn_phantom.c_str());

        // Read the first block containing parameters of phantom
        objId = MD1.firstRowId();
        MD1.getValue(MDL_DIMENSIONS_3D, TempVec, objId);
        if (TempVec.size()<3)
            REPORT_ERROR(ERR_ARG_MISSING, MDL::label2Str(MDL_DIMENSIONS_3D) + " problems with project dimensions");
        xdim = (int)TempVec[0];
        ydim = (int)TempVec[1];
        zdim = (int)TempVec[2];
        if (!MD1.getValue(MDL_PHANTOM_BGDENSITY, Background_Density, objId))
            Background_Density = 0;
        if (!MD1.getValue(MDL_SCALE, scale, objId))
            scale = 1;
        if (apply_scale)
        {
            xdim = (int) CEIL(scale * xdim);
            ydim = (int) CEIL(scale * ydim);
            zdim = (int) CEIL(scale * zdim);
            current_scale = 1;
        }
        else
            current_scale = scale;

        // Read the second block
        for (auto& FeatureRow: MD2)
        {
            if(!FeatureRow.getValue(MDL_PHANTOM_FEATURE_TYPE, feat_type))
                REPORT_ERROR(ERR_ARG_MISSING, MDL::label2Str(MDL_PHANTOM_FEATURE_TYPE) + " feature type not present");
            if (feat_type == "sph")
            {
                sph = new Sphere;
                feat = sph;
                sph->read(FeatureRow);
            }
            else if (feat_type == "blo")
            {
                blo = new Blob;
                feat = blo;
                blo->read(FeatureRow);
            }
            else if (feat_type == "gau")
            {
                gau = new Gaussian;
                feat = gau;
                gau->read(FeatureRow);
            }
            else if (feat_type == "cyl")
            {
                cyl = new Cylinder;
                feat = cyl;
                cyl->read(FeatureRow);
            }
            else if (feat_type == "dcy")
            {
                dcy = new DCylinder;
                feat = dcy;
                dcy->read(FeatureRow);
            }
            else if (feat_type == "cub")
            {
                cub = new Cube;
                feat = cub;
                cub->read(FeatureRow);
            }
            else if (feat_type == "ell")
            {
                ell = new Ellipsoid;
                feat = ell;
                ell->read(FeatureRow);
            }
            else if (feat_type == "con")
            {
                con = new Cone;
                feat = con;
                con->read(FeatureRow);
            }
            else
                REPORT_ERROR(ERR_ARG_INCORRECT, MDL::label2Str(MDL_PHANTOM_FEATURE_TYPE) + "Unknown feature type");
            if (apply_scale)
            {
                scaled_feat = feat->scale(scale);
                scaled_feat->center = scaled_feat->center * scale;
                delete feat;

                // Store feature
                VF.push_back(scaled_feat);
            }
            else
                VF.push_back(feat);
        }
    }
    else
    {
        // Open Volume Description File
        if ((fh_phantom = fopen(fn_phantom.c_str(), "r")) == nullptr)
            REPORT_ERROR(ERR_IO_NOTOPEN, (std::string)"Phantom::read: Cannot open the phantom file: "
                         + fn_phantom);
        fn = fn_phantom;

        size_t lineNumber = 0;
        // Read the file
        while (fgets(line, 256, fh_phantom) != nullptr)
        {
            ++lineNumber;
            if (line[0] == 0)
                continue;
            if (line[0] == '#')
                continue;
            if (line[0] == '\n')
                continue;

            // Read volume dimensions and global density .........................
            if (Global_Feature_Read == 0)
            {
                Global_Feature_Read = 1;
                stat = sscanf(line, "%d %d %d %lf %lf", &xdim, &ydim, &zdim,
                              &Background_Density, &scale);
                if (stat < 3)
                    REPORT_ERROR(ERR_IO_NOREAD, "Phantom::read: check the volume"
                                 " dimensions and global density in volume description file");
                if (stat <= 3)
                    Background_Density = 0;
                if (stat <= 4)
                    scale = 1;
                if (apply_scale)
                {
                    xdim = (int) CEIL(scale * xdim);
                    ydim = (int) CEIL(scale * ydim);
                    zdim = (int) CEIL(scale * zdim);
                    current_scale = 1;
                }
                else
                    current_scale = scale;
                continue;
            }

            // Read feature description ..........................................
            stat = sscanf(line, "%s", straux);
            feat_type = straux;
            if (stat != 1)
                REPORT_ERROR(ERR_IO_NOREAD, formatString("Phantom::read: Not correct feature type in line number %ld : \n%s",lineNumber, line));

            if (feat_type == "sph")
            {
                sph = new Sphere;
                feat = sph;
                sph->readCommon(line);
                sph->read_specific(line);
            }
            else if (feat_type == "blo")
            {
                blo = new Blob;
                feat = blo;
                blo->readCommon(line);
                blo->read_specific(line);
            }
            else if (feat_type == "gau")
            {
                gau = new Gaussian;
                feat = gau;
                gau->readCommon(line);
                gau->read_specific(line);
            }
            else if (feat_type == "cyl")
            {
                cyl = new Cylinder;
                feat = cyl;
                cyl->readCommon(line);
                cyl->read_specific(line);
            }
            else if (feat_type == "dcy")
            {
                dcy = new DCylinder;
                feat = dcy;
                dcy->readCommon(line);
                dcy->read_specific(line);
            }
            else if (feat_type == "cub")
            {
                cub = new Cube;
                feat = cub;
                cub->readCommon(line);
                cub->read_specific(line);
            }
            else if (feat_type == "ell")
            {
                ell = new Ellipsoid;
                feat = ell;
                ell->readCommon(line);
                ell->read_specific(line);
            }
            else if (feat_type == "con")
            {
                con = new Cone;
                feat = con;
                con->readCommon(line);
                con->read_specific(line);
            }
            else
                REPORT_ERROR(ERR_IO_NOREAD, (std::string)"Phantom::read: Unknown feature type: " + line);

            // Scale and Store feature
            if (apply_scale)
            {
                scaled_feat = feat->scale(scale);
                scaled_feat->center = scaled_feat->center * scale;
                delete feat;

                // Store feature
                VF.push_back(scaled_feat);
            }
            else
                VF.push_back(feat);
        }
        fclose(fh_phantom);
        phantom_scale = scale;
    }
}

rotate()

void Phantom::rotate

(

const Matrix2D< double > &

E

)

Rotate. Rotate a phantom after a given 3D rotation.

Definition at line 2486 of file phantom.cpp.

{
    for (size_t i = 0; i < VF.size(); i++)
        VF[i]->rotate(E);
}

selfApplyGeometry()

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

Apply a general geometric transformation. The transformation must be defined by a 4x4 matrix that can be generated using the geometric functions in xmippGeometry or xmippMatrix2D. The inv argument is either IS_INV or IS_NOT_INV. Matrices coming out from xmippGeometry for rotations, shifts, ... are not INV.

An exception is thrown if the matrix A is not valid.

Definition at line 2493 of file phantom.cpp.

{
    if ((MAT_XSIZE(A) != 4) || (MAT_YSIZE(A) != 4))
        REPORT_ERROR(ERR_MATRIX_SIZE, "Apply_geom3D: geometrical transformation is not 4x4");
    if (A.isIdentity())
        return;
    Matrix2D<double> T;
    if (inv == xmipp_transformation::IS_INV)
        T = A.inv();
    else
        T = A;

    for (size_t i = 0; i < VF.size(); i++)
        VF[i]->selfApplyGeometry(T);
}

shift()

void Phantom::shift	(	double	shiftX,
		double	shiftY,
		double	shiftZ
	)

Shift. Shift all features in the phantom a given amount of voxels. See Feature::shift.

Definition at line 2479 of file phantom.cpp.

{
    for (size_t i = 0; i < VF.size(); i++)
        VF[i]->shift(shiftX, shiftY, shiftZ);
}

sketch_in()

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

Sketch the surface of the phantom in the volume. This function allows you to draw only the surface of every feature (see Feature::sketch_in to see when a voxel is said to belong to the feature surface). The input volume might be cleaned (resized and the logical origin set at the physical center) or not according to the labels CLEAN or DONT_CLEAN (by default). The grey level for those voxels can be defined (by default, 2) and the colour is assigned

Definition at line 2467 of file phantom.cpp.

{
    if (clean)
    {
        V.resize(zdim, ydim, xdim);
        V.setXmippOrigin();
    }
    for (size_t i = 0; i < VF.size(); i++)
        VF[i]->sketch_in(V, colour);
}

surface()

void Phantom::surface	(	double	z0,
		double	radius,
		int	direction,
		Image< double > *	P
	)		const

Surface of the phantom as in AFM. This function produces the surface of a phantom in a given image when measuring rays go POS_NEG or NEG_POS (from Z positive to Z negative, and viceversa). If rays reach point z0 and the phantom is not reached the surface for that point is z0, otherwise it is the z of the last voxel which is still outside the phantom.

This function throws an exception if z0 is outside the volume definition. If z0 is equal to zdim then it is not used.

The radius is the probe radius for the surface generation.

If the output image is not resized, it is resized to the Y and X dimensions of the phantom.

Definition at line 2565 of file phantom.cpp.

{
    if (z0 != zdim)
        if (z0 < FIRST_XMIPP_INDEX(zdim) || z0 > LAST_XMIPP_INDEX(zdim))
            REPORT_ERROR(ERR_INDEX_OUTOFBOUNDS, "Phantom::surface: z0 outside phantom");
#ifdef DEBUG

    std::cout << "Direction: " << direction << std::endl;
    std::cout << "z0:        " << z0        << std::endl;
    std::cout << "zdim:      " << zdim      << std::endl;
#endif

    Matrix1D<double> aux1(3);
    Matrix1D<double> aux2(3);
    Matrix1D<double> aux3(3);
    Matrix1D<double> r(3);
    if (XSIZE((*P)()) == 0)
    {
        (*P)().resize(ydim, xdim);
        (*P)().setXmippOrigin();
    }
    FOR_ALL_ELEMENTS_IN_ARRAY2D(IMGMATRIX(*P))
    {
#ifdef DEBUG
        std::cout << "Processing (" << i << "," << j << ")" << std::endl;
#endif
        // Init ray
        int k;
        if (direction == POS_NEG)
            k = LAST_XMIPP_INDEX(zdim) + 1;
        else
            k = FIRST_XMIPP_INDEX(zdim) - 1;
        bool finished;
        finished = false;

#ifdef DEBUG

        std::cout << "Initial k=" << k << std::endl;
#endif
        // Check that it is not inside and move ray
        // at the end k takes the right value for the image
        while (!finished)
        {
            VECTOR_R3(r, j, i, k);
#ifdef DEBUG

            std::cout << "Checking " << r.transpose() << std::endl;
#endif
            // If it is inside move a step backward and finish
            if (any_feature_intersects_sphere(r, radius, aux1, aux2, aux3))
            {
                finished = true;
                if (direction == POS_NEG)
                    k++;
                else
                    k--;
            }
            else
            {
                // Else, move until you find z0
                if (z0 != zdim)
                    if (direction == POS_NEG)
                    {
                        k--;
                        if (k < z0)
                        {
                            finished = true;
                            k = CEIL(z0);
                        }
                    }
                    else
                    {
                        k++;
                        if (k > z0)
                        {
                            finished = true;
                            k = FLOOR(z0);
                        }
                    }
                // or you reach the end of the volume
                else
                    if (direction == POS_NEG)
                    {
                        k--;
                        if (k < FIRST_XMIPP_INDEX(zdim))
                        {
                            finished = true;
                        }
                    }
                    else
                    {
                        k++;
                        if (k > LAST_XMIPP_INDEX(zdim))
                        {
                            finished = true;
                        }
                    }
            }
        }

        IMGPIXEL(*P, i, j) = k;
    }
}

volume()

double Phantom::volume

(

)

const

Return the volume of all the features.

Definition at line 2079 of file phantom.cpp.

{
    double retval = 0;
    for (size_t i = 0; i < VF.size(); i++)
        retval += VF[i]->volume();
    return retval;
}

voxel_inside_any_feat() [1/2]

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

Speeded up voxel inside any feature. This function tells you if the voxel of size 1x1x1 whose center is at position r is inside any of the features. A voxel is said to be inside a feature if it shares at least 1 corner with the feature. This function returns the first feature of the list containing the voxel. Features are numbered as 1, 2, ... If the voxel is not in any feature the function returns 0.

aux1 and aux2, are two vectors of dimension 3. They must be supplied in order to gain speed in the calculations. This is very useful when checking if many voxels are inside any feature. See also Feature::voxel_inside to know more.

Definition at line 2395 of file phantom.cpp.

{
    int inside;
    int current_i;
    double current_density;
    current_i = 0;
    current_density = Background_Density;
    for (size_t i = 0; i < VF.size(); i++)
    {
        inside = VF[i]->voxel_inside(r, aux1, aux2);
        if (inside != 0 && VF[i]->density > current_density)
        {
            current_i = i + 1;
            current_density = VF[i]->density;
        }
    }
    return current_i;
}

voxel_inside_any_feat() [2/2]

int Phantom::voxel_inside_any_feat ( const Matrix1D< double > &  r ) const

inline

Voxel inside any feature. The same as the previous one but you needn't supply the extra auxiliar vectors.

Definition at line 1495 of file phantom.h.

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

write()

void Phantom::write

(

const FileName &

fn_phantom = ""

)

Write a phantom file in the standard feature format. You may rename the file or not giving a different name in the write call.

Definition at line 2355 of file phantom.cpp.

{
    MetaDataVec MD1;  //MetaData for phanto global parameters
    MetaDataVec MD2;  //MetaData for Feature parameters
    std::vector<double> FCVect(3);  //For the center of feature
    size_t id;
    // Write global parameters to the first block
    std::vector<double> PCVector;  //For the center of Phantom
    MD1.setColumnFormat(false);
    id = MD1.addObject();
    PCVector.push_back(xdim);
    PCVector.push_back(ydim);
    PCVector.push_back(zdim);
    MD1.setValue(MDL_DIMENSIONS_3D, PCVector, id);
    MD1.setValue(MDL_PHANTOM_BGDENSITY, Background_Density, id);
    if (current_scale != 1)
        MD1.setValue(MDL_SCALE, current_scale, id);
    else
        MD1.setValue(MDL_SCALE, 1.0, id);
    MD1.write((std::string)"block1@"+fn_phantom.c_str(), MD_OVERWRITE);

    // Write specific parameters
    std::string SAddAssign;  // string variab for feature operation (+/=)
    for (size_t i = 0; i < VF.size(); i++)
    {
        id = MD2.addObject();
        SAddAssign = VF[i]->add_assign;
        MD2.setValue(MDL_PHANTOM_FEATURE_TYPE,VF[i]->type, id);
        MD2.setValue(MDL_PHANTOM_FEATURE_OPERATION, SAddAssign, id);
        MD2.setValue(MDL_PHANTOM_FEATURE_DENSITY, VF[i]->density, id);
        FCVect[0] = XX(VF[i]->center);
        FCVect[1] = YY(VF[i]->center);
        FCVect[2] = ZZ(VF[i]->center);
        MD2.setValue(MDL_PHANTOM_FEATURE_CENTER, FCVect, id);
        VF[i]->feat_printm(MD2, id);
    }
    MD2.write((std::string)"block2@"+fn_phantom.c_str(), MD_APPEND);
}

Friends And Related Function Documentation

operator<<

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

friend

Show a phantom file. The more descriptive format is used instead of the standard Feature format

Definition at line 2343 of file phantom.cpp.

{
    std::cout << "Phantom ---------------------------------------\n";
    std::cout << "Dimensions: " << P.xdim << " x " << P.ydim << " x " << P.zdim << std::endl;
    std::cout << "Background density: " << P.Background_Density << std::endl;
    std::cout << "phantom_scale : " << P.phantom_scale << std::endl;
    for (size_t i = 0; i < P.VF.size(); i++)
        o << P.VF[i];
    return o;
}

Member Data Documentation

Background_Density

double Phantom::Background_Density

Final volume background density.

Definition at line 1368 of file phantom.h.

current_scale

double Phantom::current_scale

Has been the scale applied?

Definition at line 1371 of file phantom.h.

fn

FileName Phantom::fn

Filename.

Definition at line 1356 of file phantom.h.

param_file_scale

double Phantom::param_file_scale

Param file scale.

Definition at line 1374 of file phantom.h.

phantom_scale

double Phantom::phantom_scale

Param file scale.

Definition at line 1377 of file phantom.h.

VF

std::vector<Feature*> Phantom::VF

List with the features.

Definition at line 1380 of file phantom.h.

xdim

int Phantom::xdim

Final volume X dimension.

Definition at line 1359 of file phantom.h.

ydim

int Phantom::ydim

Final volume Y dimension.

Definition at line 1362 of file phantom.h.

zdim

int Phantom::zdim

Final volume Z dimension.

Definition at line 1365 of file phantom.h.

---

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

• xmipp/libraries/data/phantom.h
• xmipp/libraries/data/phantom.cpp