Feature Class Reference

#include <phantom.h>

Inheritance diagram for Feature:

Collaboration diagram for Feature:

Public Member Functions
virtual void	prepare ()=0
void	assign (const Feature &F)
virtual void	rotate (const Matrix2D< double > &E)
virtual void	rotate_center (const Matrix2D< double > &E)
virtual int	point_inside (const Matrix1D< double > &r, Matrix1D< double > &aux) const =0
int	point_inside (const Matrix1D< double > &r) const
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
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
Feature *	encircle (double radius=0) const
virtual Feature *	scale (double factor) const =0
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
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)
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 (MDRow &row)
virtual void	read_specific (const std::vector< double > &vector)=0

Public Attributes
std::string	type
char	add_assign
double	density
Matrix1D< double >	center
double	max_distance

Friends
std::ostream &	operator<< (std::ostream &o, const Feature *F)

Detailed Description

Feature superclass. This is a superclass from which all features (cones, cylinders, ...) inherit. It contains all general information common to all feature classes, then the subclasses give the specific parameters for the feature.

Definition at line 93 of file phantom.h.

Member Function Documentation

assign()

void Feature::assign

(

const Feature &

F

)

Another function for assigmnet.

Definition at line 92 of file phantom.cpp.

{
    *this = F;
}

background()

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

Return a pointer to a feature which is the background of the actual one. You can specify two ways of computing the background, either as a sphere surrounding the feature or as an enlarged version of the feature. This is chosen giving the modes ENLARGE_MODE or SPHERE_MODE. Depending on the mode used the background parameter is understood as the sphere radius or as the scaling factor.

Definition at line 1662 of file phantom.cpp.

{
    switch (back_mode)
    {
    case ENLARGE_MODE:
        return scale(back_param);
        break;
    case SPHERE_MODE:
        return encircle(back_param);
        break;
    default:
        REPORT_ERROR(ERR_VALUE_INCORRECT, "Feature::background: mode not supported");
        break;
    }
}

corners()

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

Define 3D corners for a feature. This function returns two Z3 points where the feature is confined. The volume borders are taken into account and you might make a for using these two values like this: \Ex:

F.corners(V,corner1,corner2);
for (int k=ZZ(corner1); k<=ZZ(corner2); k++)
    for (int i=YY(corner1); i<=YY(corner2); i++)
        for (int j=XX(corner1); j<=XX(corner2); j++) {
            ...
}

Definition at line 967 of file phantom.cpp.

{
    corner1.resize(3);
    corner2.resize(3);
    XX(corner1) = XMIPP_MAX(FLOOR(XX(center) - max_distance), STARTINGX(V));
    YY(corner1) = XMIPP_MAX(FLOOR(YY(center) - max_distance), STARTINGY(V));
    ZZ(corner1) = XMIPP_MAX(FLOOR(ZZ(center) - max_distance), STARTINGZ(V));
    XX(corner2) = XMIPP_MIN(CEIL(XX(center) + max_distance), FINISHINGX(V));
    YY(corner2) = XMIPP_MIN(CEIL(YY(center) + max_distance), FINISHINGY(V));
    ZZ(corner2) = XMIPP_MIN(CEIL(ZZ(center) + max_distance), FINISHINGZ(V));
}

density_inside()

virtual double Feature::density_inside ( const Matrix1D< double > &  r, Matrix1D< double > &  aux  ) const

pure virtual

Speeded up density inside a feature, VIRTUAL!!. This function MUST be implemented for each subclass with NON constant density as blobs and tells you if a point is inside the feature or not plus give you information about the density of the feature at that point. If the point is inside returns 1 multiplied by the density of the NORMALIZED feature and if not, returns 0. The point is given in the vector r, and aux is an auxiliar vector with dimension 3 (must be externally resized) used to do some computations. This auxiliar vector must be supplied in order to gain speed as no memory allocating and freeing is needed.

Implemented in Cone, Ellipsoid, Cube, DCylinder, Cylinder, Gaussian, Blob, and Sphere.

draw_in()

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

Draw a feature in a Volume. This function draws a feature in a volume (*V). The word "colour" is a little misleading as the colour really refers to a grey-level global density. Ie, the default mode (INTERNAL) draws the feature with its own density and discards the "colour" given in the function call. However, you may change this colour, set a new global density and draw the volume with this given density. This option is useful to generate easily a labelled volume.

The Add-Assign behaviour is the one of the feature if the colour is internally defined, or Assign if the colour is externally given. In the assign behaviour a voxel value is assigned to the volume if the voxel there is smaller than the value we pretend to assign.

A voxel is drawn with a colour proportional to the number of Vertex inside the features. The volume is not cleaned at the beginning. \ Ex:

f.draw_in(V)              --> Internal density
f.draw_in(V,INTERNAL,0.5) --> Internal density, 0.5 is discarded
f.draw_in(V,EXTERNAL,0.5) --> External density=0.5

Definition at line 983 of file phantom.cpp.

{
    Matrix1D<double>   aux1(3);
    Matrix1D<double>   aux2(3);
    Matrix1D<double>   corner1(3);
    Matrix1D<double>   corner2(3);
    Matrix1D<double>   r(3);
    int                add;
    double             inside;
    double             final_colour;

    if (colour_mode == INTERNAL)
    {
        final_colour = density;
        add = add_assign == '+';
    }
    else
    {
        final_colour = colour;
        add = 0;
    }

    corners(V, corner1, corner2);
#ifdef DEBUG

    std::cout << "Drawing \n";
    std::cout << this;
    std::cout << "colour_mode=" << colour_mode << std::endl;
    std::cout << "add_assign= " << add_assign  << std::endl;
    std::cout << "add=        " << add         << std::endl;
    std::cout << "   Corner 1" << corner1.transpose() << std::endl;
    std::cout << "   Corner 2" << corner2.transpose() << std::endl;
#endif

    FOR_ALL_ELEMENTS_IN_ARRAY3D_BETWEEN(corner1, corner2)
    {
        inside = voxel_inside_by_normalized_density(r, aux1, aux2);
#ifdef DEBUG
        //int condition=(ZZ(r)==-12) && (YY(r)==1);
        int condition = 1;
        if (condition)
            std::cout << "   r=" << r.transpose() << " inside= " << inside;
#endif

        if (inside != 0)
        {
            double drawing_colour = final_colour * inside / 8;
            if (add)
                Vr += drawing_colour;
            else
                Vr  = drawing_colour; // It does not select the maximum between Vr and drawing_colour anymore -> it fails when adding less dense features
#ifdef DEBUG

            if (condition)
                std::cout << "   V(r)=" << VOLVOXEL(V, (int)ZZ(r), (int)YY(r), (int)XX(r));
#endif

        }
#ifdef DEBUG
        if (condition)
            std::cout << std::endl;
#endif

    }
}

encircle()

Feature * Feature::encircle

(

double

radius = 0

)

const

Produce a sphere envolving the feature. This function returns a pointer to a feature (a sphere) with the same center, density, and feature behaviour as the given feature. The radius could be given in voxel units or if you leave it 0 then the radius is computed as 1.5 times the max_distance of this feature

Definition at line 1644 of file phantom.cpp.

{
    Sphere *f;
    f = new Sphere;

    if (radius == 0)
        radius = 1.5 * max_distance;

    f->type         = "sph";
    f->add_assign   = add_assign;
    f->density      = density;
    f->center       = center;
    f->max_distance = radius;
    f->radius       = radius;

    return (Feature *)f;
}

feat_printf()

virtual void Feature::feat_printf ( FILE *  fh ) const

pure virtual

Print the feature in the Feature format, VIRTUAL!!!. This function prints the feature in the Standard Feature format (readable by this library). Notice that the standard format is different for each specific feature. See Phantoms for more information about the supported file format.

Implemented in Cone, Ellipsoid, Cube, DCylinder, Cylinder, Gaussian, Blob, and Sphere.

feat_printm()

virtual void Feature::feat_printm ( MetaData &  MD, size_t  id  )

pure virtual

Implemented in Cone, Ellipsoid, Cube, DCylinder, Cylinder, Gaussian, Blob, and Sphere.

intersection() [1/2]

virtual double Feature::intersection ( const Matrix1D< double > &  direction, const Matrix1D< double > &  passing_point, Matrix1D< double > &  r, Matrix1D< double > &  u  ) const

pure virtual

Speeded Up intersection of a feature with a ray, VIRTUAL!!!. This function returns the length of the intersection between the ray defined by its direction and a passing point and the actual feature (whose center and dimensions are known by itself).

r and u are auxiliar vectors of dimension 3, they must be supplied in order to gain speed. r is the passing point expressed in the feature coordinate system, and u is the direction in the same coordinate system.

Implemented in Cone, Ellipsoid, Cube, DCylinder, Cylinder, Gaussian, Blob, and Sphere.

intersection() [2/2]

double Feature::intersection ( const Matrix1D< double > &  direction, const Matrix1D< double > &  passing_point  ) const

inline

Intersection of a feature with a ray. This function does the same as the previous one but you needn't provide extra auxiliar vectors.

Definition at line 274 of file phantom.h.

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

intersects_sphere() [1/2]

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

Speeded up sphere intersecting feature. This function returns TRUE if a sphere of a given radius has any voxel inside this feature. r is the center of the sphere in R3. This speeded up function needs two vectors with dimension 3 externally resized.

Definition at line 936 of file phantom.cpp.

{
    double radius2 = radius * radius;
    bool intersects = false;
    for (int k = FLOOR(ZZ(r) - radius); k <= CEIL(ZZ(r) + radius) && !intersects; k++)
    {
        auto dk=(double) k;
        double distk2=(dk - ZZ(r))*(dk - ZZ(r));
        for (int i = FLOOR(YY(r) - radius); i <= CEIL(YY(r) + radius) && !intersects; i++)
        {
            auto di=(double) i;
            double distki2=distk2+(di - YY(r))*(di - YY(r));
            for (int j = FLOOR(XX(r) - radius); j <= CEIL(XX(r) + radius) && !intersects; j++)
            {
                auto dj=(double) j;
                if (distki2+(dj - XX(r))*(dj - XX(r))>radius2)
                    continue;
                VECTOR_R3(aux3, j, i, k);
                intersects = voxel_inside(aux3, aux1, aux2);
            }
        }
    }
    return intersects;
}

intersects_sphere() [2/2]

int 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 222 of file phantom.h.

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

mean_variance_in_plane()

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

Mean and variance in a given plane. Given a plane z=z0, this function returns the mean and variance values of the volume (*V) in the voxels inside this feature. A voxel is considered to belong to the feature if it is totally inside the feature. If the plane is outside the volume scope the result is mean=variance=0

Definition at line 1929 of file phantom.cpp.

{
    double sum1 = 0;
    double sum2 = 0;
    double no_points = 0;
    Matrix1D<double> r(3);
    Matrix1D<double> aux1(3);
    Matrix1D<double> aux2(3);

    mean = 0;
    var = 0;
    if (z < STARTINGZ(VOLMATRIX(*V)) || z > FINISHINGZ(VOLMATRIX(*V)))
        return;

    ZZ(r) = z;
    for (YY(r) = STARTINGY(VOLMATRIX(*V)); YY(r) <= FINISHINGY(VOLMATRIX(*V)); YY(r)++)
        for (XX(r) = STARTINGX(VOLMATRIX(*V)); XX(r) <= FINISHINGX(VOLMATRIX(*V)); XX(r)++)
        {
            if (voxel_inside(r, aux1, aux2) == 8)
            {
                double voxel = VOLVOXEL(*V, (int)ZZ(r), (int)YY(r), (int)XX(r));
                sum1 += voxel;
                sum2 += voxel * voxel;
                no_points++;
            }
        }
    if (no_points != 0)
    {
        mean = sum1 / no_points;
        var  = sum2 / no_points - mean * mean;
    }
}

point_inside() [1/2]

virtual int Feature::point_inside ( const Matrix1D< double > &  r, Matrix1D< double > &  aux  ) const

pure virtual

Speeded up point inside a feature, VIRTUAL!!. This function MUST be implemented for each subclass and tells you if a point is inside the feature or not. If the point is inside returns 1 and if not, returns 0. The point is given in the vector r, and aux is an auxiliar vector with dimension 3 (must be externally resized) used to do some computations. This auxiliar vector must be supplied in order to gain speed as no memory allocating and freeing is needed.

Implemented in Cone, Ellipsoid, Cube, DCylinder, Cylinder, Gaussian, Blob, and Sphere.

point_inside() [2/2]

int Feature::point_inside ( const Matrix1D< double > &  r ) const

inline

Point inside a feature. This function is based in the previous one. It makes the same but you needn't supply the auxiliar vector.

Definition at line 160 of file phantom.h.

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

prepare()

virtual void Feature::prepare ( )

pure virtual

Prepare feature for work. This function computes the maximum distance and possibly the Euler and inverse Euler matrices.

Implemented in Cone, Ellipsoid, Cube, DCylinder, Cylinder, Gaussian, Blob, and Sphere.

project_to()

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

Project feature onto a projection plane. Projection is a class itself which has got inside the direction of projection. The projection plane is not cleaned (set all values to 0) when entering this function, but the projection of this feature is added to the already stored one. The projection plane is supposed to have its logical origin in the right place (normally, the middle of the image) before entering the function. An oversampling technique is used in order to produce better projections, the projection in one pixel is computed from the value of 4 points near the pixel, you can modify the subsampling value (actually 2x2=4) by changing a constant at the beginning of the function.

The matrix VP (3x3) is used to define how to relate a point in the 3D universal coordinate to the projection plane. So a point ru in the universal coordinate system, projects to VP*r. PV must be the inverse of the forwarding projection matrix.

Definition at line 1293 of file phantom.cpp.

{
constexpr float SUBSAMPLING = 2;                  // for every measure 2x2 line
    // integrals will be taken to
    // avoid numerical errors
constexpr float SUBSTEP = 1/(SUBSAMPLING*2.0);

    Matrix1D<double> origin(3);
    Matrix1D<double> direction;
    VP.getRow(2, direction);
    direction.selfTranspose();
    Matrix1D<double> corner1(3);
    Matrix1D<double> corner2(3);
    Matrix1D<double> act(3);
    SPEED_UP_temps012;

    // Find center of the feature in the projection plane ...................
    // Step 1). Project the center to the plane, the result is in the
    //          universal coord system
    M3x3_BY_V3x1(origin, VP, center);

    //   Matrix1D<double> origin_debug(3);
    //   Uproject_to_plane(center,P.direction,0,origin_debug);

    //#define DEBUG_LITTLE
#ifdef DEBUG_LITTLE

    std::cout << "Actual feature\n"     << this << std::endl;
    std::cout << "center              " << center.transpose() << std::endl;
    std::cout << "VP matrix\n"          << VP << std::endl;
    std::cout << "P.direction         " << P.direction.transpose() << std::endl;
    std::cout << "direction           " << direction.transpose() << std::endl;
    std::cout << "P.euler matrix      " << P.euler << std::endl;
    std::cout << "max_distance        " << max_distance << std::endl;
    std::cout << "origin              " << origin.transpose() << std::endl;
    //      std::cout << "origin_debug (Univ.coord) " << origin_debug.transpose() << std::endl;
#endif
    /*
       // Step 2). Express this projected center in the projection coord system
       M3x3_BY_V3x1(origin_debug,P.euler,origin_debug);
    //   if (A!=NULL) M2x2_BY_V2x1(origin,*A,origin_);
       #ifdef DEBUG_LITTLE
          std::cout << "origin (Proj.coord) " << origin_debug.transpose() << std::endl;
       #endif
    */

    // Find limits for projection ...........................................
    // Choose corners for the projection of this feature. It is supposed
    // to have at the worst case a projection of size max_distance
    VECTOR_R3(corner1, max_distance, max_distance, max_distance);
    VECTOR_R3(corner2, -max_distance, -max_distance, -max_distance);
#ifdef DEBUG_LITTLE

    std::cout << "Corner1 : " << corner1.transpose() << std::endl
    << "Corner2 : " << corner2.transpose() << std::endl;
#endif

    box_enclosing(corner1, corner2, VP, corner1, corner2);
    //   if (A!=NULL) {
    //      rectangle_enclosing(corner1,corner2,*A,corner1,corner2);
#ifdef DEBUG_LITTLE

    std::cout << "Corner1 moves to : " << corner1.transpose() << std::endl
    << "Corner2 moves to : " << corner2.transpose() << std::endl;
#endif
    //   }

    V3_PLUS_V3(corner1, origin, corner1);
    V3_PLUS_V3(corner2, origin, corner2);
#ifdef DEBUG_LITTLE

    std::cout << "Corner1 finally is : " << corner1.transpose() << std::endl
    << "Corner2 finally is : " << corner2.transpose() << std::endl;
#endif
    /*
       Matrix1D<double> corner1_debug(2),corner2_debug(2);
       VECTOR_R2(corner1_debug, max_distance, max_distance);
       VECTOR_R2(corner2_debug,-max_distance,-max_distance);
       #ifdef DEBUG_LITTLE
          std::cout << "Corner1_debug : " << corner1_debug.transpose() << std::endl
               << "Corner2_debug : " << corner2_debug.transpose() << std::endl;
       #endif
       V2_PLUS_V2(corner1_debug,origin_debug,corner1_debug);
       V2_PLUS_V2(corner2_debug,origin_debug,corner2_debug);
       #ifdef DEBUG_LITTLE
          std::cout << "Corner1_debug finally is : " << corner1_debug.transpose() << std::endl
               << "Corner2_debug finally is : " << corner2_debug.transpose() << std::endl;
       #endif
    */
    // Discard not necessary components
    corner1.resize(2);
    corner2.resize(2);

    // Clip to image size
    sortTwoVectors(corner1, corner2);
    XX(corner1) = CLIP(ROUND(XX(corner1)), STARTINGX(P()), FINISHINGX(P()));
    YY(corner1) = CLIP(ROUND(YY(corner1)), STARTINGY(P()), FINISHINGY(P()));
    XX(corner2) = CLIP(ROUND(XX(corner2)), STARTINGX(P()), FINISHINGX(P()));
    YY(corner2) = CLIP(ROUND(YY(corner2)), STARTINGY(P()), FINISHINGY(P()));

#ifdef DEBUG_LITTLE

    std::cout << "corner1      " << corner1.transpose() << std::endl;
    std::cout << "corner2      " << corner2.transpose() << std::endl;
    std::cout.flush();
#endif

    // Check if there is something to project
    if (XX(corner1) == XX(corner2))
        return;
    if (YY(corner1) == YY(corner2))
        return;

    // Study the projection for each point in the projection plane ..........
    // (u,v) are in the deformed projection plane (if any deformation)
    for (auto v = (int)YY(corner1); v <= (int)YY(corner2); v++)
        for (auto u = (int)XX(corner1); u <= (int)XX(corner2); u++)
        {
            double length = 0;
#ifdef DEBUG_EVEN_MORE

            std::cout << "Studying point (" << u << "," << v << ")\n";
            std::cout.flush();
#endif

            // Perform subsampling ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
            double u0 = u - (int)(SUBSAMPLING / 2.0) * SUBSTEP;
            double v0 = v - (int)(SUBSAMPLING / 2.0) * SUBSTEP;
            double actv = v0;
            for (int subv = 0; subv < SUBSAMPLING; subv++)
            {
                double actu = u0;
                for (int subu = 0; subu < SUBSAMPLING; subu++)
                {
                    // Compute the coordinates of point (subu,subv) which is
                    // within the plane in the universal coordinate system
                    XX(act) = actu;
                    YY(act) = actv;
                    ZZ(act) = 0;
                    //               if (Ainv!=NULL) M2x2_BY_V2x1(act,*Ainv,act);
                    //               M3x3_BY_V3x1(act,P.eulert,act);
                    M3x3_BY_V3x1(act, PV, act);

                    // Compute the intersection of a ray which passes through
                    // this point and its direction is perpendicular to the
                    // projection plane
                    double possible_length = intersection(direction, act);
                    if (possible_length > 0)
                        length += possible_length;

#ifdef DEBUG_EVEN_MORE

                    std::cout << "Averaging at (" << actu << "," << actv << ")\n";
                    std::cout << "   which in univ. coords is " << act.transpose() << std::endl;
                    std::cout << "   intersection there " << possible_length << std::endl;
#endif
                    // Prepare for next iteration
                    actu += SUBSTEP * 2.0;
                }
                actv += SUBSTEP * 2.0;
            }
            length /= (SUBSAMPLING * SUBSAMPLING);
#ifdef DEBUG

            std::cout << "Final value added at position (" << u << "," << v << ")="
            << length << std::endl;
#endif

            // Add at the correspondent pixel the found intersection ,,,,,,,,,,
            IMGPIXEL(P, v, u) += length * density;
        }
}

read()

void Feature::read

(

MDRow &

row

)

Definition at line 213 of file phantom.cpp.

{
    readCommon(row);
    std::vector<double> VectSpecific;  // Vector for specific parameters of feature
    row.getValue(MDL_PHANTOM_FEATURE_SPECIFIC, VectSpecific);
    read_specific(VectSpecific);
}

read_specific()

virtual void Feature::read_specific ( const std::vector< double > &  vector )

pure virtual

Read a feature from a file, VIRTUAL!!!. The format must be the one given in Phantoms, and each subclass must implement its own I/O routines. These routines must fill only the non common part of the feature description, but they receive the whole line with the description.

Implemented in Cone, Ellipsoid, Cube, DCylinder, Cylinder, Gaussian, Blob, and Sphere.

readCommon() [1/2]

void Feature::readCommon

(

char *

line

)

Read common part of the feature description. The common part is the feature type, the behaviour, density and center. The description is passed as a line. Exceptions are thrown if the description doesn't conform the standard specification.

Definition at line 177 of file phantom.cpp.

{
    int        stat;
    char       straux[6];
    center.resize(3);
    stat = sscanf(line, "%s %c %lf %lf %lf %lf",
                  straux,
                  &add_assign,
                  &density,
                  &(XX(center)),
                  &(YY(center)),
                  &(ZZ(center)));
    if (stat != 6)
        REPORT_ERROR(ERR_IO_NOREAD,
                     (std::string)"Error when reading common part of feature: " + line);
    type = straux;
}

readCommon() [2/2]

void Feature::readCommon

(

MDRow &

row

)

Definition at line 196 of file phantom.cpp.

{
    center.resize(3);
    std::vector <double> VecFeatureCenter;  // Keep the center of the feature
    std::string s_op;  // As no label for char in MD_TYPE  (for add/assign)
    if (!row.getValue(MDL_PHANTOM_FEATURE_TYPE,type) ||
        !row.getValue(MDL_PHANTOM_FEATURE_OPERATION,s_op) ||
        !row.getValue(MDL_PHANTOM_FEATURE_DENSITY,density) ||
        !row.getValue(MDL_PHANTOM_FEATURE_CENTER,VecFeatureCenter))
        REPORT_ERROR(ERR_ARG_MISSING, (std::string)"Error when reading common part of feature");
    XX(center) = VecFeatureCenter[0];
    YY(center) = VecFeatureCenter[1];
    ZZ(center) = VecFeatureCenter[2];
    add_assign = s_op[0];
}

rotate()

void Feature::rotate ( const Matrix2D< double > &  E )

virtual

Rotate the whole feature. Rotate this feature using a rotation matrix. The center as well as the feature itself is rotated.

Reimplemented in Oriented_Feature.

Definition at line 159 of file phantom.cpp.

{
    rotate_center(E);
}

rotate_center()

void Feature::rotate_center ( const Matrix2D< double > &  E )

virtual

Rotate only the center. Rotate the center of this feature only around the phantom center

Definition at line 151 of file phantom.cpp.

{
    Matrix2D<double> inverse_angles_matrix;
    inverse_angles_matrix = E.inv();
    center = inverse_angles_matrix * center;
    //center=E*center;
}

scale()

virtual Feature* Feature::scale ( double  factor ) const

pure virtual

Return a scaled version of this feature, VIRTUAL!!!. This function returns a pointer to a feature (of the same type as the feature for which the function was called, ie, if you scale a sphere the result is a sphere, if you scale a cone, the result is a cone, ...) that is a scaled version of the actual feature (see each specific implementation to see which is the relatioship between the two features). This function is useful to define backgrounds with the same shape of the given feature.

Implemented in Cone, Ellipsoid, Cube, DCylinder, Cylinder, Gaussian, Blob, and Sphere.

selfApplyGeometry()

void Feature::selfApplyGeometry

(

const Matrix2D< double > &

A

)

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 matrix must be the desired transformation (i.e., new coordinate=A*old_coordinate. Don't worry because the selfApplyGeometry of Phantom take cares of passing to this function the appropriate matrix. No check is done about the size of A.

Only the center is transformed, the feature will keep the same size.

Definition at line 1078 of file phantom.cpp.

{
    Matrix1D<double> r(4);
    XX(r) = XX(center);
    YY(r) = YY(center);
    ZZ(r) = ZZ(center);
    r(3) = 1;
    r = A * r;
    XX(center) = XX(r);
    YY(center) = YY(r);
    ZZ(center) = ZZ(r);
}

shift()

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

Shift. Shift the feature a given amount of voxels. The new center is the old center plus the given shift, and that's all

Definition at line 1070 of file phantom.cpp.

{
    XX(center) += shiftX;
    YY(center) += shiftY;
    ZZ(center) += shiftZ;
}

sketch_in()

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

Draw the surface of the feature. This function draws the surface of the feature at the given volume. A voxel is said to belong to the surface if the number of corners inside the feature meets 1<=n<=7. The default gray_level with which voxels will be drawn is 2 (normally higher than the usual volume grey levels, and the behaviour of the voxels drawn is Assign.

Definition at line 1051 of file phantom.cpp.

{
    Matrix1D<double>   aux1(3);
    Matrix1D<double>   aux2(3);
    Matrix1D<double>   corner1(3);
    Matrix1D<double>   corner2(3);
    Matrix1D<double>   r(3);
    int                inside;

    corners(V, corner1, corner2);
    FOR_ALL_ELEMENTS_IN_ARRAY3D_BETWEEN(corner1, corner2)
    {
        inside = voxel_inside(r, aux1, aux2);
        if (inside != 0 && inside != 8)
            A3D_ELEM(V, (int)ZZ(r), (int)YY(r), (int)XX(r)) = colour;
    }
}

volume()

virtual double Feature::volume ( ) const

pure virtual

Volume of the feature in (voxel units)� , VIRTUAL!!!. This function returns the volume of each feature supposing that a voxel is of size 1x1x1.

Implemented in Cone, Ellipsoid, Cube, DCylinder, Cylinder, Gaussian, Blob, and Sphere.

voxel_inside() [1/2]

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

Speeded up voxel inside a feature. A voxel is compound of 8 subvoxels. This function returns the number of subvoxel centers falling inside the feature. The voxel size is supposed to be 1, and r is the center of the voxel in R3. This speeded up function needs two vectors with dimension 3 externally resized.

Definition at line 845 of file phantom.cpp.

{

    // The subvoxels are visited following a Gray code, so the number
    // of operations is minimized
    XX(aux1) = XX(r) + 0.25;
    YY(aux1) = YY(r) + 0.25;
    ZZ(aux1) = ZZ(r) + 0.25; // 000
    int inside = point_inside(aux1, aux2);
    DEBUG_SHOW;
    ZZ(aux1) -= 0.5;
    inside += point_inside(aux1, aux2);             // 001
    DEBUG_SHOW;
    YY(aux1) -= 0.5;
    inside += point_inside(aux1, aux2);             // 011
    DEBUG_SHOW;
    ZZ(aux1) += 0.5;
    inside += point_inside(aux1, aux2);             // 010
    DEBUG_SHOW;
    XX(aux1) -= 0.5;
    inside += point_inside(aux1, aux2);             // 110
    DEBUG_SHOW;
    ZZ(aux1) -= 0.5;
    inside += point_inside(aux1, aux2);             // 111
    DEBUG_SHOW;
    YY(aux1) += 0.5;
    inside += point_inside(aux1, aux2);             // 101
    DEBUG_SHOW;
    ZZ(aux1) += 0.5;
    inside += point_inside(aux1, aux2);             // 100
    DEBUG_SHOW;
    return inside;
}

voxel_inside() [2/2]

int Feature::voxel_inside ( const Matrix1D< double > &  r ) const

inline

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

Definition at line 192 of file phantom.h.

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

voxel_inside_by_normalized_density()

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

A voxel is compound of 8 subvoxels. This function returns the number of subvoxel centers falling inside the feature multiplied by the normalized feature density. That is, the value of the density is always 1 at the origin The voxel size is supposed to be 1, and r is the center of the voxel in R3. This speeded up function needs two vectors with dimension 3 externally resized.

Definition at line 880 of file phantom.cpp.

{
#ifdef NEVER
    if (type == "blo")
    {
        std::cout << "den=" <<   density_inside(r, aux2) << std::endl;
        return(density_inside(r, aux2));
    }
    else
        return((double)voxel_inside(r, aux1, aux2));
#endif
    // The subvoxels are visited following a Gray code, so the number
    // of operations is minimized
    XX(aux1) = XX(r) + 0.25;
    YY(aux1) = YY(r) + 0.25;
    ZZ(aux1) = ZZ(r) + 0.25; // 000
    double inside = (double)point_inside(aux1, aux2)
                    * density_inside(r, aux2);
    DEBUG_SHOW;
    ZZ(aux1) -= 0.5;
    inside += (double)point_inside(aux1, aux2)
              * density_inside(r, aux2); // 001
    DEBUG_SHOW;
    YY(aux1) -= 0.5;
    inside += (double)point_inside(aux1, aux2)
              * density_inside(r, aux2);  // 011
    DEBUG_SHOW;
    ZZ(aux1) += 0.5;
    inside += (double)point_inside(aux1, aux2)
              * density_inside(r, aux2); // 010
    DEBUG_SHOW;
    XX(aux1) -= 0.5;
    inside += (double)point_inside(aux1, aux2)
              * density_inside(r, aux2); // 110
    DEBUG_SHOW;
    ZZ(aux1) -= 0.5;
    inside += (double)point_inside(aux1, aux2)
              * density_inside(r, aux2); // 111
    DEBUG_SHOW;
    YY(aux1) += 0.5;
    inside += (double)point_inside(aux1, aux2)
              * density_inside(r, aux2); // 101
    DEBUG_SHOW;
    ZZ(aux1) += 0.5;
    inside += (double)point_inside(aux1, aux2)
              * density_inside(r, aux2); // 100
    DEBUG_SHOW;
    return inside;

}

Friends And Related Function Documentation

operator<<

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

friend

Show feature not in the standard format but more informatively. This function is based on the std::cout << ... of each subclass. First shows the common part of the feature and then its specific part. Be careful that you must show a pointer to the feature!! \ Ex: Sphere sphere; std::cout << (Feature *) &sphere;

Definition at line 569 of file phantom.cpp.

{
    if (F != nullptr)
    {
        o << "Feature --------" << std::endl;
        o << "   type:        " << F->type << std::endl;
        o << "   add_assign:  " << F->add_assign << std::endl;
        o << "   density:     " << F->density << std::endl;
        o << "   center:      " << F->center.transpose() << std::endl;
        if (F->type == "sph")
            o << *((Sphere *) F);
        else if (F->type == "blo")
            o << *((Blob *) F);
        else if (F->type == "gau")
            o << *((Gaussian *) F);
        else if (F->type == "cyl")
            o << *((Cylinder *) F);
        else if (F->type == "dcy")
            o << *((DCylinder *) F);
        else if (F->type == "cub")
            o << *((Cube *) F);
        else if (F->type == "ell")
            o << *((Ellipsoid *) F);
        else if (F->type == "con")
            o << *((Cone *) F);
    }
    return o;
}

Member Data Documentation

add_assign

char Feature::add_assign

Feature behaviour. This flag indicates how the feature behaves inside the voxel volume. If this flag is set to '+' then the voxels occupied by this feature are incremented with the feature value. If the flag is set to '=' then the voxels are set to to the same value of the feature.

Definition at line 108 of file phantom.h.

center

Matrix1D<double> Feature::center

Center of the feature. The center of the feature is understood differently according to the specific class, see them to know exactly how this value is interpreted.

Definition at line 119 of file phantom.h.

density

double Feature::density

Feature Density. Density of the feature in grey level values. It needn't be between 0 and 1, or 0 and 255. What is more, it can be even negative, and so you can build holes inside volumes.

Definition at line 114 of file phantom.h.

max_distance

double Feature::max_distance

Maximum distance from the center. This value is a precalculated and tells the maximum distance from any point belonging to the feature to its center. This is used to speed up some functions not considering voxels which we know are beyond the scope of this feature.

Definition at line 126 of file phantom.h.

type

std::string Feature::type

Feature type. A three letters string telling what kind of feature this object is. For example, "cyl", "con", "cub", ... See the specific classes to know exactly which is each label.

Definition at line 101 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