C++ VDOT函数代码示例

/**
 * R E C _ C U R V E
 *
 * Return the "curvature" of the cylinder.  If an endplate,
 * pick a principle direction orthogonal to the normal, and
 * indicate no curvature.  Otherwise, compute curvature.
 * Normal must have been computed before calling this routine.
 */
void
rt_rec_curve(struct curvature *cvp, struct hit *hitp, struct soltab *stp)
{
    struct rec_specific *rec =
	(struct rec_specific *)stp->st_specific;
    vect_t uu;
    fastf_t ax, bx, q;

    switch (hitp->hit_surfno) {
	case REC_NORM_BODY:
	    /* This could almost certainly be simpler if we used
	     * inverse A rather than inverse A squared, right Ed?
	     */
	    VMOVE(cvp->crv_pdir, rec->rec_Hunit);
	    VSUB2(uu, hitp->hit_point, rec->rec_V);
	    cvp->crv_c1 = 0;
	    ax = VDOT(uu, rec->rec_A) * rec->rec_iAsq;
	    bx = VDOT(uu, rec->rec_B) * rec->rec_iBsq;
	    q = sqrt(ax * ax * rec->rec_iAsq + bx * bx * rec->rec_iBsq);
	    cvp->crv_c2 = - rec->rec_iAsq * rec->rec_iBsq / (q*q*q);
	    break;
	case REC_NORM_TOP:
	case REC_NORM_BOT:
	    bn_vec_ortho(cvp->crv_pdir, hitp->hit_normal);
	    cvp->crv_c1 = cvp->crv_c2 = 0;
	    break;
	default:
	    bu_log("rt_rec_curve: bad surfno %d\n", hitp->hit_surfno);
	    break;
    }
}

/* planar image map */
color image_plane_texture(const vector * hit, const texture * tx, const ray * ry) {
  vector pnt;
  flt u, v, miprad, maxscale;
  standard_texture * tex = (standard_texture *) tx;
 
  pnt.x=hit->x - tex->ctr.x;
  pnt.y=hit->y - tex->ctr.y;
  pnt.z=hit->z - tex->ctr.z;

  VDOT(u, tex->uaxs, pnt);
  VDOT(v, tex->vaxs, pnt); 

  u = u * tex->scale.x;
  u = u + tex->rot.x;
  u = u - ((int) u);
  if (u < 0.0) u += 1.0;

  v = v * tex->scale.y;
  v = v + tex->rot.y;
  v = v - ((int) v);
  if (v < 0.0) v += 1.0;

  
  maxscale = (fabs(tex->scale.x) > fabs(tex->scale.y)) ? 
             tex->scale.x : tex->scale.y;

  miprad = 0.05 * ry->opticdist * fabs(maxscale);

  return MIPMap(tex->img, u, v, miprad);
}

static void light_intersect(point_light * l, ray * ry) {
  flt b, disc, t1, t2, temp;
  vector V;

  /* Lights do not cast shadows.. */
  if (ry->flags & RT_RAY_SHADOW)
    return;

  VSUB(l->ctr, ry->o, V);
  VDOT(b, V, ry->d); 
  VDOT(temp, V, V);  

  disc=b*b + l->rad*l->rad - temp;

  if (disc<=0.0) return;
  disc=sqrt(disc);

  t2=b+disc;
  if (t2 <= SPEPSILON) 
    return;
  add_intersection(t2, (object *) l, ry);  

  t1=b-disc;
  if (t1 > SPEPSILON) 
    add_intersection(t1, (object *) l, ry);  
}

/* planar image map */
color image_plane_texture(vector * hit, texture * tex, ray * ry) {
    vector pnt;
    flt u,v;

    pnt.x=hit->x - tex->ctr.x;
    pnt.y=hit->y - tex->ctr.y;
    pnt.z=hit->z - tex->ctr.z;

    VDOT(u, tex->uaxs, pnt);
    /*  VDOT(len, tex->uaxs, tex->uaxs);
      u = u / sqrt(len); */

    VDOT(v, tex->vaxs, pnt);
    /*  VDOT(len, tex->vaxs, tex->vaxs);
      v = v / sqrt(len); */


    u = u * tex->scale.x;
    u = u + tex->rot.x;
    u = fmod(u, 1.0);
    if (u < 0.0) u += 1.0;

    v = v * tex->scale.y;
    v = v + tex->rot.y;
    v = fmod(v, 1.0);
    if (v < 0.0) v += 1.0;

    return ImageMap((rawimage *)tex->img, u, v);
}

static void cylinder_intersect(const cylinder * cyl, ray * ry) {
  vector rc, n, D, O;  
  flt t, s, tin, tout, ln, d; 

  rc.x = ry->o.x - cyl->ctr.x;
  rc.y = ry->o.y - cyl->ctr.y;
  rc.z = ry->o.z - cyl->ctr.z; 

  VCross(&ry->d, &cyl->axis, &n);

  ln=SQRT(n.x*n.x + n.y*n.y + n.z*n.z);    /* finish length calculation */

  if (ln == 0.0) {  /* ray is parallel to the cylinder.. */
    VDOT(d, rc, cyl->axis);         
    D.x = rc.x - d * cyl->axis.x; 
    D.y = rc.y - d * cyl->axis.y;
    D.z = rc.z - d * cyl->axis.z;
    VDOT(d, D, D);
    d = SQRT(d);
    tin = -FHUGE;
    tout = FHUGE;
    /* if (d <= cyl->rad) then ray is inside cylinder.. else outside */
  }

  n.x /= ln;
  n.y /= ln;
  n.z /= ln;

  VDOT(d, rc, n);
  d = FABS(d); 

  if (d <= cyl->rad) {  /* ray intersects cylinder.. */
    VCross(&rc, &cyl->axis, &O);
    VDOT(t, O, n);
    t = - t / ln;
    VCross(&n, &cyl->axis, &O); 

    ln = SQRT(O.x*O.x + O.y*O.y + O.z*O.z);
    O.x /= ln;
    O.y /= ln;
    O.z /= ln;

    VDOT(s, ry->d, O);
    s = FABS(SQRT(cyl->rad*cyl->rad - d*d) / s);
    tin = t - s;
    ry->add_intersection(tin, (object *) cyl, ry); 
    tout = t + s;
    ry->add_intersection(tout, (object *) cyl, ry);
  }
}

static void fcylinder_intersect(cylinder * cyl, ray * ry) {
  vector rc, n, O, hit, tmp2, ctmp4;
  flt t, s, tin, tout, ln, d, tmp, tmp3;
 
  rc.x = ry->o.x - cyl->ctr.x;  
  rc.y = ry->o.y - cyl->ctr.y;
  rc.z = ry->o.z - cyl->ctr.z;
 
  VCross(&ry->d, &cyl->axis, &n);
 
  VDOT(ln, n, n);
  ln=sqrt(ln);    /* finish length calculation */
 
  if (ln == 0.0) {  /* ray is parallel to the cylinder.. */
    return;       /* in this case, we want to miss or go through the "hole" */
  }
 
  VNorm(&n);
  VDOT(d, rc, n);
  d = fabs(d);
 
  if (d <= cyl->rad) {  /* ray intersects cylinder.. */
    VCross(&rc, &cyl->axis, &O);
    VDOT(t, O, n);
    t = - t / ln;
    VCross(&n, &cyl->axis, &O);
    VNorm(&O);
    VDOT(s, ry->d, O);
    s = fabs(sqrt(cyl->rad*cyl->rad - d*d) / s);
    tin = t - s;

    RAYPNT(hit, (*ry), tin); 

    ctmp4=cyl->axis;
    VNorm(&ctmp4);

    tmp2.x = hit.x - cyl->ctr.x;   
    tmp2.y = hit.y - cyl->ctr.y;   
    tmp2.z = hit.z - cyl->ctr.z;   

    VDOT(tmp,  tmp2, ctmp4);
    VDOT(tmp3, cyl->axis, cyl->axis);

    if ((tmp > 0.0) && (tmp < sqrt(tmp3))) 
      add_intersection(tin, (object *) cyl, ry);
    tout = t + s;

    RAYPNT(hit, (*ry), tout); 

    tmp2.x = hit.x - cyl->ctr.x;   
    tmp2.y = hit.y - cyl->ctr.y;   
    tmp2.z = hit.z - cyl->ctr.z;   

    VDOT(tmp,  tmp2, ctmp4); 
    VDOT(tmp3, cyl->axis, cyl->axis);

    if ((tmp > 0.0) && (tmp < sqrt(tmp3))) 
      add_intersection(tout, (object *) cyl, ry);
  }
}

/*
 * R E F R A C T
 *
 * Compute the refracted ray 'v_2' from the incident ray 'v_1' with
 * the refractive indices 'ri_2' and 'ri_1' respectively.
 * Using Schnell's Law:
 *
 * theta_1 = angle of v_1 with surface normal
 * theta_2 = angle of v_2 with reversed surface normal
 * ri_1 * sin(theta_1) = ri_2 * sin(theta_2)
 *
 * sin(theta_2) = ri_1/ri_2 * sin(theta_1)
 *
 * The above condition is undefined for ri_1/ri_2 * sin(theta_1)
 * being greater than 1, and this represents the condition for total
 * reflection, the 'critical angle' is the angle theta_1 for which
 * ri_1/ri_2 * sin(theta_1) equals 1.
 *
 * Returns TRUE if refracted, FALSE if reflected.
 *
 * Note:  output (v_2) can be same storage as an input.
 */
HIDDEN int
rr_refract(vect_t v_1, vect_t norml, double ri_1, double ri_2, vect_t v_2)
{
    vect_t w, u;
    fastf_t beta;

    if (NEAR_ZERO(ri_1, 0.0001) || NEAR_ZERO(ri_2, 0.0001)) {
	bu_log("rr_refract:ri1=%g, ri2=%g\n", ri_1, ri_2);
	beta = 1;
    } else {
	beta = ri_1/ri_2;		/* temp */
	if (beta > 10000) {
	    bu_log("rr_refract:  beta=%g\n", beta);
	    beta = 1000;
	}
    }
    VSCALE(w, v_1, beta);
    VCROSS(u, w, norml);

    /*
     *	|w X norml| = |w||norml| * sin(theta_1)
     *	        |u| = ri_1/ri_2 * sin(theta_1) = sin(theta_2)
     */
    if ((beta = VDOT(u, u)) > 1.0) {
	/* Past critical angle, total reflection.
	 * Calculate reflected (bounced) incident ray.
	 */
	if (R_DEBUG&RDEBUG_REFRACT) bu_log("rr_refract: reflected.  ri1=%g ri2=%g beta=%g\n",
					   ri_1, ri_2, beta);
	VREVERSE(u, v_1);
	beta = 2 * VDOT(u, norml);
	VSCALE(w, norml, beta);
	VSUB2(v_2, w, u);
	return 0;		/* reflected */
    } else {
	/*
	 * 1 - beta = 1 - sin(theta_2)^^2
	 *	    = cos(theta_2)^^2.
	 * beta = -1.0 * cos(theta_2) - Dot(w, norml).
	 */
	if (R_DEBUG&RDEBUG_REFRACT) bu_log("rr_refract: refracted.  ri1=%g ri2=%g beta=%g\n",
					   ri_1, ri_2, beta);
	beta = -sqrt(1.0 - beta) - VDOT(w, norml);
	VSCALE(u, norml, beta);
	VADD2(v_2, w, u);
	return 1;		/* refracted */
    }
    /* NOTREACHED */
}

/**
 * Calculate a bounding RPP for an ARBN
 */
int
rt_arbn_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *UNUSED(tol)) {
    size_t i, j, k;
    struct rt_arbn_internal *aip;
    RT_CK_DB_INTERNAL(ip);
    aip = (struct rt_arbn_internal *)ip->idb_ptr;
    RT_ARBN_CK_MAGIC(aip);

    VSETALL((*min), INFINITY);
    VSETALL((*max), -INFINITY);

    /* Discover all vertices, use to calculate RPP */
    for (i = 0; i < aip->neqn-2; i++) {
	for (j = i+1; j < aip->neqn-1; j++) {
	    double dot;

	    /* If normals are parallel, no intersection */
	    dot = VDOT(aip->eqn[i], aip->eqn[j]);
	    if (((dot) <= -SMALL_FASTF) ? (NEAR_EQUAL((dot), -1.0,  RT_DOT_TOL)) : (NEAR_EQUAL((dot), 1.0, RT_DOT_TOL))) continue;

	    /* Have an edge line, isect with higher numbered planes */
	    for (k = j + 1; k < aip->neqn; k++) {
		size_t m;
		size_t next_k;
		point_t pt;

		next_k = 0;

		if (bn_mkpoint_3planes(pt, aip->eqn[i], aip->eqn[j], aip->eqn[k]) < 0) continue;

		/* See if point is outside arb */
		for (m = 0; m < aip->neqn; m++) {
		    if (i == m || j == m || k == m)
			continue;
		    if (VDOT(pt, aip->eqn[m])-aip->eqn[m][3] > RT_LEN_TOL) {
			next_k = 1;
			break;
		    }
		}
		if (next_k != 0) continue;

		VMINMAX((*min), (*max), pt);
	    }
	}
    }

    return 0;
}

/**
 * Make sure that the lu and fu orientation flags are consistent with
 * the geometric arrangement of the vertices and the faceuse normal.
 */
void
nmg_ck_lu_orientation(struct loopuse *lu, const struct bn_tol *tolp)
{
    struct faceuse *fu;
    plane_t fu_peqn;
    plane_t lu_peqn;
    fastf_t dot;

    NMG_CK_LOOPUSE(lu);
    fu = lu->up.fu_p;		/* parent had better be faceuse */
    NMG_CK_FACEUSE(fu);

    NMG_GET_FU_PLANE(fu_peqn, fu);
    nmg_loop_plane_newell(lu, lu_peqn);

    dot = VDOT(fu_peqn, lu_peqn);

    if (NEAR_ZERO(dot, tolp->perp))
	return;		/* can't determine geometric orientation */


    if (dot < 0.0) {
	bu_log("nmg_ck_lu_orientation() lu=%p, dot=%g, fu_orient=%s, lu_orient=%s\n",
	       (void *)lu, dot,
	       nmg_orientation(fu->orientation),
	       nmg_orientation(lu->orientation)
	    );
	bu_bomb("nmg_ck_lu_orientation() loop orientation flags do not match geometry\n");
    }
}

void
render_flos_work(render_t *render, struct tie_s *tie, struct tie_ray_s *ray, vect_t *pixel) {
    struct tie_id_s id, tid;
    vect_t vec;
    fastf_t angle;
    struct render_flos_s *rd;

    rd = (struct render_flos_s *)render->data;

    if (tie_work(tie, ray, &id, render_hit, NULL) != NULL) {
	VSET(*pixel, 0.0, 0.5, 0.0);
    } else
	return;

    VSUB2(vec, ray->pos, id.pos);
    VUNITIZE(vec);
    angle = VDOT(vec, id.norm);

    /* Determine if direct line of sight to fragment */
    VMOVE(ray->pos, rd->frag_pos);
    VSUB2(ray->dir, id.pos, rd->frag_pos);
    VUNITIZE(ray->dir);

    if (tie_work(tie, ray, &tid, render_hit, NULL)) {
	if (fabs (id.pos[0] - tid.pos[0]) < TIE_PREC
	    && fabs (id.pos[1] - tid.pos[1]) < TIE_PREC
	    && fabs (id.pos[2] - tid.pos[2]) < TIE_PREC)
	{
	    VSET(*pixel, 1.0, 0.0, 0.0);
	}
    }

    VSCALE(*pixel, *pixel, (0.5+angle*0.5));
}

void
rt_arbn_volume(fastf_t *volume, const struct rt_db_internal *ip)
{
    struct poly_face *faces;
    struct rt_arbn_internal *aip = (struct rt_arbn_internal *)ip->idb_ptr;
    size_t i;

    *volume = 0.0;
    /* allocate array of face structs */
    faces = (struct poly_face *)bu_calloc(aip->neqn, sizeof(struct poly_face), "rt_arbn_volume: faces");
    for (i = 0; i < aip->neqn; i++) {
	/* allocate array of pt structs, max number of verts per faces = (# of faces) - 1 */
	faces[i].pts = (point_t *)bu_calloc(aip->neqn - 1, sizeof(point_t), "rt_arbn_volume: pts");
    }
    rt_arbn_faces_area(faces, aip);
    for (i = 0; i < aip->neqn; i++) {
	vect_t tmp;

	/* calculate volume of pyramid */
	VSCALE(tmp, faces[i].plane_eqn, faces[i].area);
	*volume += VDOT(faces[i].pts[0], tmp)/3;
    }
    for (i = 0; i < aip->neqn; i++) {
	bu_free((char *)faces[i].pts, "rt_arbn_volume: pts");
    }
    bu_free((char *)faces, "rt_arbn_volume: faces");
}

static void cylinder_normal(cylinder * cyl, vector * pnt, ray * incident, vector * N) {
  vector a,b,c;
  flt t;

  VSub((vector *) pnt, &(cyl->ctr), &a);

  c=cyl->axis;

  VNorm(&c);
 
  VDOT(t, a, c);

  b.x = c.x * t + cyl->ctr.x; 
  b.y = c.y * t + cyl->ctr.y;
  b.z = c.z * t + cyl->ctr.z;

  VSub(pnt, &b, N); 
  VNorm(N);

  if (VDot(N, &(incident->d)) > 0.0)  { /* make cylinder double sided */
    N->x=-N->x;
    N->y=-N->y;
    N->z=-N->z;
  } 
}

/*
    computePlaneEq--takes as input a pointer 'polygon' to an
    arbitrary 3D polygon and returns in 'plane' the normalized
    (unit normal) plane equation of the polygon.
    Newell's method (see "Newell's Method for Computing the Plane
    Equation of a Polygon" in this volume) is used for the
    computation.
*/
static void computePlaneEq(Polygon* polygon, Plane plane)
{
    int i;
    Point refpt;
    Normal normal;
    float *u, *v, len;

    /* first, compute the normal of the input polygon */
    VZERO(normal);
    VZERO(refpt);
    for(i = 0; i < polygon->numVertices; i++) {
        u = polygon->vertices[i];
        v = polygon->vertices[(i + 1) % polygon->numVertices];
        normal[X] += (u[Y] - v[Y]) * (u[Z] + v[Z]);
        normal[Y] += (u[Z] - v[Z]) * (u[X] + v[X]);
        normal[Z] += (u[X] - v[X]) * (u[Y] + v[Y]);
        VINCR(refpt, u);
    }
    /* then, compute the normalized plane equation using the
       arithmetic average of the vertices of the input polygon to
       determine its last coefficient. Note that, for efficiency,
       'refpt' stores the sum of the vertices rather than the
       actual average; the division by 'polygon->numVertices' is
       carried out together with the normalization when computing
       'plane[D]'.
    */
    len = VNORM(normal);
    plane[X] = normal[X] / len;
    plane[Y] = normal[Y] / len;
    plane[Z] = normal[Z] / len;
    len *= polygon->numVertices;
    plane[D] = -VDOT(refpt, normal) / len;
}

static int
hit_headon(struct application *ap, struct partition *PartHeadp)
{
    register char diff_solid;
    vect_t	diff;
    register fastf_t len;

    if (PartHeadp->pt_forw->pt_forw != PartHeadp)
	Tcl_AppendResult(interp, "hit_headon: multiple partitions\n", (char *)NULL);

    VJOIN1(PartHeadp->pt_forw->pt_inhit->hit_point, ap->a_ray.r_pt,
	   PartHeadp->pt_forw->pt_inhit->hit_dist, ap->a_ray.r_dir);
    VSUB2(diff, PartHeadp->pt_forw->pt_inhit->hit_point, aim_point);

    diff_solid = (FIRST_SOLID(sp) !=
		  PartHeadp->pt_forw->pt_inseg->seg_stp->st_dp);
    len = MAGNITUDE(diff);

    if (	NEAR_ZERO(len, epsilon)
		||
		( diff_solid &&
		  VDOT(diff, ap->a_ray.r_dir) > 0 )
	)
	return(1);
    else
	return(0);
}

static void cylinder_normal(const cylinder * cyl, const vector * pnt, const ray * incident, vector * N) {
  vector a, b;
  flt t, invlen, invlen2;

  a.x = pnt->x - cyl->ctr.x;
  a.y = pnt->y - cyl->ctr.y;
  a.z = pnt->z - cyl->ctr.z;

  b=cyl->axis;

  invlen = 1.0 / SQRT(b.x*b.x + b.y*b.y + b.z*b.z);
  b.x *= invlen;
  b.y *= invlen;
  b.z *= invlen;
 
  VDOT(t, a, b);

  N->x = pnt->x - (b.x * t + cyl->ctr.x); 
  N->y = pnt->y - (b.y * t + cyl->ctr.y);
  N->z = pnt->z - (b.z * t + cyl->ctr.z);

  invlen2 = 1.0 / SQRT(N->x*N->x + N->y*N->y + N->z*N->z);
  N->x *= invlen2;
  N->y *= invlen2;
  N->z *= invlen2;

  /* Flip surface normal to point toward the viewer if necessary */
  if (VDot(N, &(incident->d)) > 0.0)  {
    N->x=-N->x;
    N->y=-N->y;
    N->z=-N->z;
  }
}