BRL-CAD
#include "common.h"
#include <stddef.h>
#include <string.h>
#include <math.h>
#include "bio.h"
#include "bu/cv.h"
#include "vmath.h"
#include "db.h"
#include "nmg.h"
#include "rtgeom.h"
#include "raytrace.h"
#include "../../librt_private.h"
Include dependency graph for epa.c:

Go to the source code of this file.

Data Structures

struct  epa_specific
 

Macros

#define EPA_NORM_BODY   (1) /* compute normal */
 
#define EPA_NORM_TOP   (2) /* copy epa_N */
 
#define ELLOUT(n)   ov+(n-1)*3
 

Functions

int rt_epa_bbox (struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *tol)
 
int rt_epa_prep (struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
 
void rt_epa_print (const struct soltab *stp)
 
int rt_epa_shot (struct soltab *stp, struct xray *rp, struct application *ap, struct seg *seghead)
 
void rt_epa_norm (struct hit *hitp, struct soltab *stp, struct xray *rp)
 
void rt_epa_curve (struct curvature *cvp, struct hit *hitp, struct soltab *stp)
 
void rt_epa_uv (struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp)
 
void rt_epa_free (struct soltab *stp)
 
int rt_epa_adaptive_plot (struct rt_db_internal *ip, const struct rt_view_info *info)
 
int rt_epa_plot (struct bu_list *vhead, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol, const struct rt_view_info *info)
 
void rt_ell_norms (fastf_t *ov, fastf_t *A, fastf_t *B, fastf_t *h_vec, fastf_t t, int sides)
 
void rt_ell (fastf_t *ov, const fastf_t *V, const fastf_t *A, const fastf_t *B, int sides)
 
int rt_epa_tess (struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
 
int rt_epa_import4 (struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
 
int rt_epa_export4 (struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
 
int rt_epa_import5 (struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
 
int rt_epa_export5 (struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
 
int rt_epa_describe (struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
 
void rt_epa_ifree (struct rt_db_internal *ip)
 
int rt_epa_params (struct pc_pc_set *ps, const struct rt_db_internal *ip)
 
void rt_epa_volume (fastf_t *vol, const struct rt_db_internal *ip)
 
void rt_epa_centroid (point_t *cent, const struct rt_db_internal *ip)
 
void rt_epa_surf_area (fastf_t *area, const struct rt_db_internal *ip)
 

Variables

const struct bu_structparse rt_epa_parse []
 

Detailed Description

Intersect a ray with an Elliptical Paraboloid.

Algorithm -

Given V, H, R, and B, there is a set of points on this epa

{ (x, y, z) | (x, y, z) is on epa }

Through a series of Affine Transformations, this set of points will be transformed into a set of points on an epa located at the origin with a semi-major axis R1 along the +Y axis, a semi-minor axis R2 along the -X axis, a height H along the -Z axis, and a vertex V at the origin.

{ (x', y', z') | (x', y', z') is on epa at origin }

The transformation from X to X' is accomplished by:

X' = S(R(X - V))

where R(X) = (R1/(-|R1|)) (R2/(|R2|)) . X (H /(-|H |))

and S(X) = (1/|R1| 0 0) (0 1/|R2| 0) . X (0 0 1/|H |)

To find the intersection of a line with the surface of the epa, consider the parametric line L:

L : { P(n) | P + t(n) . D }

Call W the actual point of intersection between L and the epa. Let W' be the point of intersection between L' and the unit epa.

L' : { P'(n) | P' + t(n) . D' }

W = invR(invS(W')) + V

Where W' = k D' + P'.

If Dy' and Dz' are both 0, then there is no hit on the epa; but the end plates need checking. If there is now only 1 hit point, the top plate needs to be checked as well.

Line L' hits the infinitely long epa at W' when

A * k**2 + B * k + C = 0

where

A = Dx'**2 + Dy'**2 B = 2 * (Dx' * Px' + Dy' * Py') - Dz' C = Px'**2 + Py'**2 - Pz' - 1 b = |Breadth| = 1.0 h = |Height| = 1.0 r = 1.0

The quadratic formula yields k (which is constant):

k = [ -B +/- sqrt(B**2 - 4 * A * C)] / (2.0 * A)

Now, D' = S(R(D)) and P' = S(R(P - V))

Substituting,

W = V + invR(invS[ k *(S(R(D))) + S(R(P - V)) ]) = V + invR((k * R(D)) + R(P - V)) = V + k * D + P - V = k * D + P

Note that ``k'' is constant, and is the same in the formulations for both W and W'.

The hit at ``k'' is a hit on the canonical epa IFF Wz' <= 0.

NORMALS. Given the point W on the surface of the epa, what is the vector normal to the tangent plane at that point?

Map W onto the unit epa, i.e.: W' = S(R(W - V)).

Plane on unit epa at W' has a normal vector N' where

N' = <Wx', Wy', -.5>.

The plane transforms back to the tangent plane at W, and this new plane (on the original epa) has a normal vector of N, viz:

N = inverse[ transpose(inverse[ S o R ]) ] (N')

because if H is perpendicular to plane Q, and matrix M maps from Q to Q', then inverse[ transpose(M) ] (H) is perpendicular to Q'. Here, H and Q are in "prime space" with the unit sphere. [Somehow, the notation here is backwards]. So, the mapping matrix M = inverse(S o R), because S o R maps from normal space to the unit sphere.

N = inverse[ transpose(inverse[ S o R ]) ] (N') = inverse[ transpose(invR o invS) ] (N') = inverse[ transpose(invS) o transpose(invR) ] (N') = inverse[ inverse(S) o R ] (N') = invR o S (N')

because inverse(R) = transpose(R), so R = transpose(invR), and S = transpose(S).

Note that the normal vector produced above will not have unit length.

THE TOP PLATE.

If Dz' == 0, line L' is parallel to the top plate, so there is no hit on the top plate. Otherwise, rays intersect the top plate with k = (0 - Pz')/Dz'. The solution is within the top plate IFF Wx'**2

  • Wy'**2 <= 1.

The normal for a hit on the top plate is -Hunit.

Definition in file epa.c.