Index: src/librt/mkbundle.c
===================================================================
--- src/librt/mkbundle.c	(revision 63934)
+++ src/librt/mkbundle.c	(working copy)
@@ -27,6 +27,7 @@
 #include "common.h"
 
 #include <math.h>
+#include <string.h>
 #include "bio.h"
 
 #include "vmath.h"
@@ -51,9 +52,9 @@
  */
 
 int
-rt_raybundle_maker(struct xray *rp, double radius, const fastf_t *avec, const fastf_t *bvec, int rays_per_ring, int nring)
+rt_raybundle_maker(struct xrays *rays, const struct xray *center_ray, double radius, const fastf_t *avec, const fastf_t *bvec, int rays_per_ring, int nring)
 {
-    register struct xray *rayp = rp+1;
+    register struct xrays *xrayp;
     int ring;
     double fraction = 1.0;
     double theta;
@@ -61,8 +62,12 @@
     double radial_scale;
     int count = 0;
 
-    rp[0].index = count++;
-    rp[0].magic =RT_RAY_MAGIC;
+    BU_ALLOC(xrayp, struct xrays);
+    VMOVE(xrayp->ray.r_pt, center_ray->r_pt);
+    VMOVE(xrayp->ray.r_dir, center_ray->r_dir);
+    xrayp->ray.index = count++;
+    xrayp->ray.magic = RT_RAY_MAGIC;
+    BU_LIST_APPEND(&rays->l, &xrayp->l);
 
     for (ring=0; ring < nring; ring++) {
 	register int i;
@@ -77,12 +82,13 @@
 	    /* pt = V + cos(theta) * A + sin(theta) * B */
 	    ct = cos(theta) * radial_scale;
 	    st = sin(theta) * radial_scale;
-	    VJOIN2(rayp->r_pt, rp[0].r_pt, ct, avec, st, bvec);
-	    VMOVE(rayp->r_dir, rp[0].r_dir);
-	    rayp->index = count++;
-	    rayp->magic = RT_RAY_MAGIC;
+	    BU_ALLOC(xrayp, struct xrays);
+	    VJOIN2(xrayp->ray.r_pt, center_ray->r_pt, ct, avec, st, bvec);
+	    VMOVE(xrayp->ray.r_dir, center_ray->r_dir);
+	    xrayp->ray.index = count++;
+	    xrayp->ray.magic = RT_RAY_MAGIC;
+	    BU_LIST_APPEND(&rays->l, &xrayp->l);
 	    theta += delta;
-	    rayp++;
 	}
     }
     return count;
@@ -190,13 +196,12 @@
     VCROSS(a_dir, orig_dir, up_vector);
     VUNITIZE(a_dir);
 
-    for (y = -radius; y <= radius; y += gridsize) {
-	vect_t tmp;
-	printf("y:%f\n", y);
-	VSCALE(tmp, b_dir, y);
-	printf("y_partofit: %f,%f,%f\n", V3ARGS(tmp));
-	for (x = -radius; x <= radius; x += gridsize) {
-	    if (((x*x)/rsq + (y*y)/rsq) <= 1) {
+    /* This adds BN_TOL_DIST to various conditions because in some
+     * cases, floating-point problems can make extremely close numbers
+     * compare incorrectly. */
+    for (y = -radius; y <= radius + BN_TOL_DIST; y += gridsize) {
+	for (x = -radius; x <= radius + BN_TOL_DIST; x += gridsize) {
+	    if (((x*x)/rsq + (y*y)/rsq) <= 1 + BN_TOL_DIST) {
 		BU_ALLOC(xrayp, struct xrays);
 		VMOVE(xrayp->ray.r_pt, start);
 		VJOIN2(xrayp->ray.r_dir, orig_dir, x, a_dir, y, b_dir);
@@ -217,7 +222,7 @@
 rt_gen_frustum(struct xrays *rays, const struct xray *center_ray,
 	       const vect_t a_vec, const vect_t b_vec,
 	       const fastf_t a_theta, const fastf_t b_theta,
-	       const fastf_t a_num, const fastf_t UNUSED(b_num))
+	       const fastf_t a_num, const fastf_t b_num)
 {
     int count;
 
@@ -226,9 +231,10 @@
 
     fastf_t x, y;
 
-    fastf_t a_length = tan(a_theta);
-    fastf_t b_length = tan(b_theta);
+    fastf_t a_length = tan(a_theta) * MAGNITUDE(center_ray->r_dir);
+    fastf_t b_length = tan(b_theta) * MAGNITUDE(center_ray->r_dir);
     fastf_t a_inc = 2 * a_length / (a_num - 1);
+    fastf_t b_inc = 2 * b_length / (b_num - 1);
 
     vect_t a_dir, b_dir;
 
@@ -245,7 +251,7 @@
     /* This adds BN_TOL_DIST to the *_length variables in the
      * condition because in some cases, floating-point problems can
      * make extremely close numbers compare incorrectly. */
-    for (y = -b_length; y <= b_length + BN_TOL_DIST;) {
+    for (y = -b_length; y <= b_length + BN_TOL_DIST; y += b_inc) {
 	for (x = -a_length; x <= a_length + BN_TOL_DIST; x += a_inc) {
 	    BU_ALLOC(xrayp, struct xrays);
 	    VMOVE(xrayp->ray.r_pt, start);
@@ -288,8 +294,8 @@
     VMOVE(b_dir, b_vec);
     VUNITIZE(b_dir);
 
-    for (y = -b_length; y <= b_length; y += db) {
-	for (x = -a_length; x <= a_length; x += da) {
+    for (y = -b_length; y <= b_length + BN_TOL_DIST; y += db) {
+	for (x = -a_length; x <= a_length +  BN_TOL_DIST; x += da) {
 	    BU_ALLOC(xrayp, struct xrays);
 	    VJOIN2(xrayp->ray.r_pt, orig_start, x, a_dir, y, b_dir);
 	    VMOVE(xrayp->ray.r_dir, dir);
@@ -301,6 +307,135 @@
     return count;
 }
 
+/* Used to carry data between rt_gen_stencil() and
+ * gen_stencil_ray_hit() */
+struct gen_stencil_data {
+    struct xrays *rays;
+    int count;
+};
+
+static int
+gen_stencil_ray_hit(register struct application *ap, struct partition *PartHeadp, struct seg *segp)
+{
+    register struct partition *pp;
+    struct partition_bundle *bundle = (struct partition_bundle *)ap->a_uptr;
+    struct partition_list *new_shotline;
+
+    struct gen_stencil_data *data;
+    struct xrays *xrayp;
+
+    if ((pp=PartHeadp->pt_forw) == PartHeadp) {
+	bundle->misses++;
+	return 0;
+    }
+
+    bundle->hits++;
+
+    if (bundle->list == NULL) {
+	/*
+	 * setup partition collection
+	 */
+	BU_ALLOC(bundle->list, struct partition_list);
+	BU_LIST_INIT(&(bundle->list->l));
+    }
+
+    /* add a new partition to list */
+    BU_ALLOC(new_shotline, struct partition_list);
+
+    /* steal partition list */
+    BU_LIST_INIT_MAGIC((struct bu_list *)&new_shotline->PartHeadp, PT_HD_MAGIC);
+    BU_LIST_APPEND_LIST((struct bu_list *)&new_shotline->PartHeadp, (struct bu_list *)PartHeadp);
+
+    BU_LIST_INIT_MAGIC(&new_shotline->segHeadp.l, RT_SEG_MAGIC);
+    BU_LIST_APPEND_LIST(&new_shotline->segHeadp.l, &segp->l);
+
+    new_shotline->ap = ap;
+    BU_LIST_PUSH(&(bundle->list->l), &(new_shotline->l));
+
+    /* Since this ray hit, append it to the xrays list. */
+    data = (struct gen_stencil_data *)bundle->ap->a_uptr;
+    BU_ALLOC(xrayp, struct xrays);
+    VMOVE(xrayp->ray.r_pt, ap->a_ray.r_pt);
+    VMOVE(xrayp->ray.r_dir, ap->a_ray.r_dir);
+    xrayp->ray.index = data->count++;
+    xrayp->ray.magic = RT_RAY_MAGIC;
+    BU_LIST_APPEND(&data->rays->l, &xrayp->l);
+
+    return 1;
+}
+
+
+int rt_gen_stencil(struct xrays *rays, const struct xray *start, struct rt_i *rtip,
+		   const vect_t up_vector, const fastf_t grid_size, const int backout)
+{
+    fastf_t backout_amt = 0;
+
+    /* For generating the initial sampling of rays */
+    vect_t avec, bvec;
+    struct xray grid_center;
+    struct application ap;
+    struct application_bundle b;
+    struct xrays *xr;
+
+    struct gen_stencil_data data;
+
+    VMOVE(grid_center.r_pt, start->r_pt);
+    VMOVE(grid_center.r_dir, start->r_dir);
+    VUNITIZE(grid_center.r_dir);
+
+    if (backout) {
+	/* ``Back out'' to make sure the point rays are fired from is
+	 * outside the model */
+	fastf_t center_dist, delta;
+	vect_t center, dvec;
+
+	VADD2SCALE(center, rtip->mdl_max, rtip->mdl_min, 0.5);
+	center_dist = DIST_PT_PT(center, grid_center.r_pt);
+	VSUB2(dvec, grid_center.r_pt, center);
+	VUNITIZE(dvec);
+	delta = center_dist * VDOT(grid_center.r_dir, dvec);
+	backout_amt = rtip->rti_radius + delta;
+	VJOIN1(grid_center.r_pt, grid_center.r_pt, -backout_amt, grid_center.r_dir);
+    }
+
+    /* Set up the application/application bundle structs for
+     * raytracing the initial grid */
+    data.rays = rays;
+    data.count = 0;
+
+    RT_APPLICATION_INIT(&ap);
+    ap.a_ray = grid_center;
+    ap.a_hit = gen_stencil_ray_hit;
+    ap.a_miss = NULL;
+    ap.a_rt_i = rtip;
+    ap.a_uptr = (void *)&data;
+    b.b_ap = ap;
+    BU_LIST_INIT(&b.b_rays.l);
+    b.b_hit = NULL;
+    b.b_miss = NULL;
+
+    /* Generate the initial set of rays */
+    VMOVE(avec, up_vector);
+    VUNITIZE(avec);
+    VSCALE(avec, avec, 2 * rtip->rti_radius);
+    VCROSS(bvec, ap.a_ray.r_dir, avec);
+    VUNITIZE(bvec);
+    VSCALE(bvec, bvec, 2 * rtip->rti_radius);
+
+    rt_gen_rect(&b.b_rays, &grid_center, avec, bvec, grid_size, grid_size);
+
+    rt_shootrays(&b);
+
+    while (BU_LIST_WHILE(xr, xrays, &(b.b_rays.l))) {
+	BU_LIST_DEQUEUE(&(xr->l));
+	printf("original ray: %f,%f,%f\n", V3ARGS(xr->ray.r_pt));
+	bu_free(xr, "bundled ray");
+    }
+    return BU_LIST_PREV(xrays, &rays->l)->ray.index + 1;
+}
+
+
+
 /*
  * Local Variables:
  * mode: C
