BRL-CAD
ehy_brep.cpp
Go to the documentation of this file.
1 /* E H Y _ B R E P . C P P
2  * BRL-CAD
3  *
4  * Copyright (c) 2008-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 /** @file ehy_brep.cpp
21  *
22  * Convert an Elliptical Hyperboloid to b-rep form
23  *
24  */
25 
26 #include "common.h"
27 
28 #include "raytrace.h"
29 #include "rtgeom.h"
30 #include "brep.h"
31 
32 
33 extern "C" void
34 rt_ehy_brep(ON_Brep **b, const struct rt_db_internal *ip, const struct bn_tol *)
35 {
36  struct rt_ehy_internal *eip;
37 
39  eip = (struct rt_ehy_internal *)ip->idb_ptr;
40  RT_EHY_CK_MAGIC(eip);
41 
42  // Check the parameters
43  if (!NEAR_ZERO(VDOT(eip->ehy_Au, eip->ehy_H), RT_DOT_TOL)) {
44  bu_log("rt_ehy_brep: Au and H are not perpendicular!\n");
45  return;
46  }
47 
48  if (!NEAR_EQUAL(MAGNITUDE(eip->ehy_Au), 1.0, RT_LEN_TOL)) {
49  bu_log("rt_ehy_brep: Au not a unit vector!\n");
50  return;
51  }
52 
53  if (MAGNITUDE(eip->ehy_H) < RT_LEN_TOL
54  || eip->ehy_c < RT_LEN_TOL
55  || eip->ehy_r1 < RT_LEN_TOL
56  || eip->ehy_r2 < RT_LEN_TOL) {
57  bu_log("rt_ehy_brep: not all dimensions positive!\n");
58  return;
59  }
60 
61  if (eip->ehy_r2 > eip->ehy_r1) {
62  bu_log("rt_ehy_brep: semi-minor axis cannot be longer than semi-major axis!\n");
63  return;
64  }
65 
66  point_t p1_origin;
67  ON_3dPoint plane1_origin, plane2_origin;
68  ON_3dVector plane_x_dir, plane_y_dir;
69 
70  // First, find plane in 3 space corresponding to the bottom face of the EPA.
71 
72  vect_t x_dir, y_dir;
73 
74  VMOVE(x_dir, eip->ehy_Au);
75  VCROSS(y_dir, eip->ehy_Au, eip->ehy_H);
76  VUNITIZE(y_dir);
77 
78  VMOVE(p1_origin, eip->ehy_V);
79  plane1_origin = ON_3dPoint(p1_origin);
80  plane_x_dir = ON_3dVector(x_dir);
81  plane_y_dir = ON_3dVector(y_dir);
82  const ON_Plane ehy_bottom_plane(plane1_origin, plane_x_dir, plane_y_dir);
83 
84  // Next, create an ellipse in the plane corresponding to the edge of the ehy.
85 
86  ON_Ellipse ellipse1(ehy_bottom_plane, eip->ehy_r1, eip->ehy_r2);
87  ON_NurbsCurve* ellcurve1 = ON_NurbsCurve::New();
88  ellipse1.GetNurbForm((*ellcurve1));
89  ellcurve1->SetDomain(0.0, 1.0);
90 
91  // Generate the bottom cap
92  ON_SimpleArray<ON_Curve*> boundary;
93  boundary.Append(ON_Curve::Cast(ellcurve1));
94  ON_PlaneSurface* bp = new ON_PlaneSurface();
95  bp->m_plane = ehy_bottom_plane;
96  bp->SetDomain(0, -100.0, 100.0);
97  bp->SetDomain(1, -100.0, 100.0);
98  bp->SetExtents(0, bp->Domain(0));
99  bp->SetExtents(1, bp->Domain(1));
100  (*b)->m_S.Append(bp);
101  const int bsi = (*b)->m_S.Count() - 1;
102  ON_BrepFace& bface = (*b)->NewFace(bsi);
103  (*b)->NewPlanarFaceLoop(bface.m_face_index, ON_BrepLoop::outer, boundary, true);
104  const ON_BrepLoop* bloop = (*b)->m_L.Last();
105  bp->SetDomain(0, bloop->m_pbox.m_min.x, bloop->m_pbox.m_max.x);
106  bp->SetDomain(1, bloop->m_pbox.m_min.y, bloop->m_pbox.m_max.y);
107  bp->SetExtents(0, bp->Domain(0));
108  bp->SetExtents(1, bp->Domain(1));
109  (*b)->SetTrimIsoFlags(bface);
110  delete ellcurve1;
111 
112  // Now, the hard part. Need an elliptical hyperbolic NURBS surface
113  // First step is to create a nurbs curve.
114 
115  double intercept_calc = (eip->ehy_c)*(eip->ehy_c)/(MAGNITUDE(eip->ehy_H) + eip->ehy_c);
116  double intercept_dist = MAGNITUDE(eip->ehy_H) + eip->ehy_c - intercept_calc;
117  double intercept_length = intercept_dist - MAGNITUDE(eip->ehy_H);
118  double MX = MAGNITUDE(eip->ehy_H);
119  double MP = MX + intercept_length;
120  double w = (MX/MP)/(1-MX/MP);
121 
122  point_t ep1, ep2, ep3;
123  VSET(ep1, -eip->ehy_r1, 0, 0);
124  VSET(ep2, 0, 0, w*intercept_dist);
125  VSET(ep3, eip->ehy_r1, 0, 0);
126  ON_3dPoint onp1 = ON_3dPoint(ep1);
127  ON_3dPoint onp2 = ON_3dPoint(ep2);
128  ON_3dPoint onp3 = ON_3dPoint(ep3);
129 
130  ON_3dPointArray cpts(3);
131  cpts.Append(onp1);
132  cpts.Append(onp2);
133  cpts.Append(onp3);
134  ON_BezierCurve *bcurve = new ON_BezierCurve(cpts);
135  bcurve->MakeRational();
136  bcurve->SetWeight(1, w);
137 
138  ON_NurbsCurve* tnurbscurve = ON_NurbsCurve::New();
139  bcurve->GetNurbForm(*tnurbscurve);
140  ON_NurbsCurve* hypbnurbscurve = ON_NurbsCurve::New();
141  const ON_Interval subinterval = ON_Interval(0, 0.5);
142  tnurbscurve->GetNurbForm(*hypbnurbscurve, 0.0, &subinterval);
143 
144  // Next, rotate that curve around the height vector.
145 
146  point_t revpoint1, revpoint2;
147  VSET(revpoint1, 0, 0, 0);
148  VSET(revpoint2, 0, 0, MX);
149  ON_3dPoint rpnt1 = ON_3dPoint(revpoint1);
150  ON_3dPoint rpnt2 = ON_3dPoint(revpoint2);
151 
152  ON_Line revaxis = ON_Line(rpnt1, rpnt2);
153  ON_RevSurface* hyp_surf = ON_RevSurface::New();
154  hyp_surf->m_curve = hypbnurbscurve;
155  hyp_surf->m_axis = revaxis;
156  hyp_surf->m_angle = ON_Interval(0, 2*ON_PI);
157 
158  // Get the NURBS form of the surface
159  ON_NurbsSurface *ehycurvedsurf = ON_NurbsSurface::New();
160  hyp_surf->GetNurbForm(*ehycurvedsurf, 0.0);
161 
162  delete hyp_surf;
163  delete tnurbscurve;
164  delete bcurve;
165 
166  // Transformations
167 
168  for (int i = 0; i < ehycurvedsurf->CVCount(0); i++) {
169  for (int j = 0; j < ehycurvedsurf->CVCount(1); j++) {
170  point_t cvpt;
171  ON_4dPoint ctrlpt;
172  ehycurvedsurf->GetCV(i, j, ctrlpt);
173 
174  // Scale the control points of the
175  // resulting surface to map to the shorter axis.
176  VSET(cvpt, ctrlpt.x, ctrlpt.y * eip->ehy_r2/eip->ehy_r1, ctrlpt.z);
177 
178  // Rotate according to the directions of Au and H
179  vect_t Hu;
180  mat_t R;
181  point_t new_cvpt;
182 
183  VSCALE(Hu, eip->ehy_H, 1/MAGNITUDE(eip->ehy_H));
184  MAT_IDN(R);
185  VMOVE(&R[0], eip->ehy_Au);
186  VMOVE(&R[4], y_dir);
187  VMOVE(&R[8], Hu);
188  VEC3X3MAT(new_cvpt, cvpt, R);
189  VMOVE(cvpt, new_cvpt);
190 
191  // Translate according to V
192  vect_t scale_v;
193  VSCALE(scale_v, eip->ehy_V, ctrlpt.w);
194  VADD2(cvpt, cvpt, scale_v);
195 
196  ON_4dPoint newpt = ON_4dPoint(cvpt[0], cvpt[1], cvpt[2], ctrlpt.w);
197  ehycurvedsurf->SetCV(i, j, newpt);
198  }
199  }
200 
201  (*b)->m_S.Append(ehycurvedsurf);
202  int surfindex = (*b)->m_S.Count();
203  ON_BrepFace& face = (*b)->NewFace(surfindex - 1);
204  (*b)->FlipFace(face);
205  int faceindex = (*b)->m_F.Count();
206  (*b)->NewOuterLoop(faceindex-1);
207 }
208 
209 
210 // Local Variables:
211 // tab-width: 8
212 // mode: C++
213 // c-basic-offset: 4
214 // indent-tabs-mode: t
215 // c-file-style: "stroustrup"
216 // End:
217 // ex: shiftwidth=4 tabstop=8
#define RT_LEN_TOL
Definition: raytrace.h:169
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
#define RT_DOT_TOL
Definition: raytrace.h:170
#define VSET(a, b, c, d)
Definition: color.c:53
Header file for the BRL-CAD common definitions.
#define RT_CK_DB_INTERNAL(_p)
Definition: raytrace.h:207
#define NEAR_ZERO(val, epsilon)
Definition: color.c:55
unsigned char * bp
Definition: rot.c:56
Support for uniform tolerances.
Definition: tol.h:71
void * idb_ptr
Definition: raytrace.h:195
#define R
Definition: msr.c:55
void rt_ehy_brep(ON_Brep **b, const struct rt_db_internal *ip, const struct bn_tol *)
Definition: ehy_brep.cpp:34