BRL-CAD
hyp_brep.cpp
Go to the documentation of this file.
1 /* H Y P _ 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 hyp_brep.cpp
21  *
22  * Convert a Hyperboloid of One Sheet 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_hyp_brep(ON_Brep **b, const struct rt_db_internal *ip, const struct bn_tol *)
35 {
36  struct rt_hyp_internal *eip;
37 
39  eip = (struct rt_hyp_internal *)ip->idb_ptr;
40  RT_HYP_CK_MAGIC(eip);
41 
42  point_t p1_origin, p2_origin;
43  ON_3dPoint plane1_origin, plane2_origin;
44  ON_3dVector plane_x_dir, plane_y_dir;
45 
46  // First, find planes corresponding to the top and bottom faces - initially
47 
48  vect_t x_dir, y_dir;
49  VMOVE(x_dir, eip->hyp_A);
50  VCROSS(y_dir, eip->hyp_A, eip->hyp_Hi);
51  VREVERSE(y_dir, y_dir);
52 
53  VMOVE(p1_origin, eip->hyp_Vi);
54  plane1_origin = ON_3dPoint(p1_origin);
55  plane_x_dir = ON_3dVector(x_dir);
56  plane_y_dir = ON_3dVector(y_dir);
57  const ON_Plane hyp_bottom_plane(plane1_origin, plane_x_dir, plane_y_dir);
58 
59  VADD2(p2_origin, eip->hyp_Vi, eip->hyp_Hi);
60  plane2_origin = ON_3dPoint(p2_origin);
61  const ON_Plane hyp_top_plane(plane2_origin, plane_x_dir, plane_y_dir);
62 
63  // Next, create ellipses in the planes corresponding to the edges of the hyp
64 
65  ON_Ellipse b_ell(hyp_bottom_plane, MAGNITUDE(eip->hyp_A), eip->hyp_b);
66  ON_NurbsCurve* bcurve = ON_NurbsCurve::New();
67  b_ell.GetNurbForm((*bcurve));
68  bcurve->SetDomain(0.0, 1.0);
69 
70  ON_Ellipse t_ell(hyp_top_plane, MAGNITUDE(eip->hyp_A), eip->hyp_b);
71  ON_NurbsCurve* tcurve = ON_NurbsCurve::New();
72  t_ell.GetNurbForm((*tcurve));
73  tcurve->SetDomain(0.0, 1.0);
74 
75  // Generate the bottom cap
76  ON_SimpleArray<ON_Curve*> boundary;
77  boundary.Append(ON_Curve::Cast(bcurve));
78  ON_PlaneSurface* bp = new ON_PlaneSurface();
79  bp->m_plane = hyp_bottom_plane;
80  bp->SetDomain(0, -100.0, 100.0);
81  bp->SetDomain(1, -100.0, 100.0);
82  bp->SetExtents(0, bp->Domain(0));
83  bp->SetExtents(1, bp->Domain(1));
84  (*b)->m_S.Append(bp);
85  const int bsi = (*b)->m_S.Count() - 1;
86  ON_BrepFace& bface = (*b)->NewFace(bsi);
87  (*b)->NewPlanarFaceLoop(bface.m_face_index, ON_BrepLoop::outer, boundary, true);
88  const ON_BrepLoop* bloop = (*b)->m_L.Last();
89  bp->SetDomain(0, bloop->m_pbox.m_min.x, bloop->m_pbox.m_max.x);
90  bp->SetDomain(1, bloop->m_pbox.m_min.y, bloop->m_pbox.m_max.y);
91  bp->SetExtents(0, bp->Domain(0));
92  bp->SetExtents(1, bp->Domain(1));
93  (*b)->FlipFace(bface);
94  (*b)->SetTrimIsoFlags(bface);
95  boundary.Empty();
96  delete bcurve;
97 
98  // Generate the top cap
99  boundary.Append(ON_Curve::Cast(tcurve));
100  ON_PlaneSurface* tp = new ON_PlaneSurface();
101  tp->m_plane = hyp_top_plane;
102  tp->SetDomain(0, -100.0, 100.0);
103  tp->SetDomain(1, -100.0, 100.0);
104  tp->SetExtents(0, bp->Domain(0));
105  tp->SetExtents(1, bp->Domain(1));
106  (*b)->m_S.Append(tp);
107  int tsi = (*b)->m_S.Count() - 1;
108  ON_BrepFace& tface = (*b)->NewFace(tsi);
109  (*b)->NewPlanarFaceLoop(tface.m_face_index, ON_BrepLoop::outer, boundary, true);
110  ON_BrepLoop* tloop = (*b)->m_L.Last();
111  tp->SetDomain(0, tloop->m_pbox.m_min.x, tloop->m_pbox.m_max.x);
112  tp->SetDomain(1, tloop->m_pbox.m_min.y, tloop->m_pbox.m_max.y);
113  tp->SetExtents(0, bp->Domain(0));
114  tp->SetExtents(1, bp->Domain(1));
115  (*b)->SetTrimIsoFlags(tface);
116  delete tcurve;
117 
118  // Now, the hard part. Need an elliptical hyperbolic NURBS surface.
119  // First step is to create a nurbs curve.
120 
121  double MX = eip->hyp_b * eip->hyp_bnr;
122  point_t ep1, ep2, ep3;
123  VSET(ep1, -eip->hyp_b, 0, 0.5*MAGNITUDE(eip->hyp_Hi));
124  VSET(ep2, -MX*eip->hyp_bnr, 0, 0);
125  VSET(ep3, -eip->hyp_b, 0, -0.5*MAGNITUDE(eip->hyp_Hi));
126 
127  ON_3dPoint onp1 = ON_3dPoint(ep1);
128  ON_3dPoint onp2 = ON_3dPoint(ep2);
129  ON_3dPoint onp3 = ON_3dPoint(ep3);
130 
131  ON_3dPointArray cpts(3);
132  cpts.Append(onp1);
133  cpts.Append(onp2);
134  cpts.Append(onp3);
135  ON_BezierCurve *bezcurve = new ON_BezierCurve(cpts);
136  bezcurve->MakeRational();
137  bezcurve->SetWeight(1, bezcurve->Weight(0)/eip->hyp_bnr);
138 
139  ON_NurbsCurve* tnurbscurve = ON_NurbsCurve::New();
140  bezcurve->GetNurbForm(*tnurbscurve);
141  delete bezcurve;
142 
143  ON_3dPoint revpnt1 = ON_3dPoint(0, 0, -0.5*MAGNITUDE(eip->hyp_Hi));
144  ON_3dPoint revpnt2 = ON_3dPoint(0, 0, 0.5*MAGNITUDE(eip->hyp_Hi));
145 
146  ON_Line revaxis = ON_Line(revpnt1, revpnt2);
147  ON_RevSurface* hyp_surf = ON_RevSurface::New();
148  hyp_surf->m_curve = tnurbscurve;
149  hyp_surf->m_axis = revaxis;
150  hyp_surf->m_angle = ON_Interval(0, 2*ON_PI);
151 
152  // Get the NURBS form of the surface
153  ON_NurbsSurface *hypcurvedsurf = ON_NurbsSurface::New();
154  hyp_surf->GetNurbForm(*hypcurvedsurf, 0.0);
155  delete hyp_surf;
156 
157  for (int i = 0; i < hypcurvedsurf->CVCount(0); i++) {
158  for (int j = 0; j < hypcurvedsurf->CVCount(1); j++) {
159  point_t cvpt;
160  ON_4dPoint ctrlpt;
161  hypcurvedsurf->GetCV(i, j, ctrlpt);
162 
163  // Scale and shear
164  vect_t proj_ah;
165  vect_t proj_ax;
166  fastf_t factor;
167 
168  VPROJECT(eip->hyp_A, eip->hyp_Hi, proj_ah, proj_ax);
169  VSET(cvpt, ctrlpt.x * MAGNITUDE(proj_ax)/eip->hyp_b, ctrlpt.y, ctrlpt.z);
170  factor = VDOT(eip->hyp_A, eip->hyp_Hi)>0 ? 1.0 : -1.0;
171  cvpt[2] += factor*cvpt[0]/MAGNITUDE(proj_ax)*MAGNITUDE(proj_ah) + 0.5*MAGNITUDE(eip->hyp_Hi)*ctrlpt.w;
172 
173  // Rotate
174  vect_t Au, Bu, Hu;
175  mat_t R;
176  point_t new_cvpt;
177 
178  VSCALE(Bu, y_dir, 1/MAGNITUDE(y_dir));
179  VSCALE(Hu, eip->hyp_Hi, 1/MAGNITUDE(eip->hyp_Hi));
180  VCROSS(Au, Bu, Hu);
181  VUNITIZE(Au);
182  MAT_IDN(R);
183  VMOVE(&R[0], Au);
184  VMOVE(&R[4], Bu);
185  VMOVE(&R[8], Hu);
186  VEC3X3MAT(new_cvpt, cvpt, R);
187  VMOVE(cvpt, new_cvpt);
188 
189  // Translate
190  vect_t scale_v;
191  VSCALE(scale_v, eip->hyp_Vi, ctrlpt.w);
192  VADD2(cvpt, cvpt, scale_v);
193  ON_4dPoint newpt = ON_4dPoint(cvpt[0], cvpt[1], cvpt[2], ctrlpt.w);
194  hypcurvedsurf->SetCV(i, j, newpt);
195  }
196  }
197 
198  (*b)->m_S.Append(hypcurvedsurf);
199  int surfindex = (*b)->m_S.Count();
200  ON_BrepFace& face = (*b)->NewFace(surfindex - 1);
201  (*b)->FlipFace(face);
202  int faceindex = (*b)->m_F.Count();
203  (*b)->NewOuterLoop(faceindex-1);
204 
205 }
206 
207 
208 // Local Variables:
209 // tab-width: 8
210 // mode: C++
211 // c-basic-offset: 4
212 // indent-tabs-mode: t
213 // c-file-style: "stroustrup"
214 // End:
215 // ex: shiftwidth=4 tabstop=8
void rt_hyp_brep(ON_Brep **b, const struct rt_db_internal *ip, const struct bn_tol *)
Definition: hyp_brep.cpp:34
#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
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
double fastf_t
Definition: defines.h:300