BRL-CAD
superell_brep.cpp
Go to the documentation of this file.
1 /* S U P E R E L L _ B R E P . C P P
2  * BRL-CAD
3  *
4  * Copyright (c) 2012-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 superell_brep.cpp
21  *
22  * Convert a Superquadratic Ellipsoid 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 extern "C" {
33  void rt_ell_brep(ON_Brep **b, struct rt_db_internal *ip, const struct bn_tol *tol);
34 }
35 
36 extern "C" void
37 rt_superell_brep(ON_Brep **b, const struct rt_db_internal *ip, const struct bn_tol *tol)
38 {
39  struct rt_superell_internal *sip;
40 
42  sip = (struct rt_superell_internal *)ip->idb_ptr;
43  RT_SUPERELL_CK_MAGIC(sip);
44 
45  // First, create brep of an ellipsoid
46  struct rt_ell_internal *eip;
47  BU_ALLOC(eip, struct rt_ell_internal);
48  eip->magic = RT_ELL_INTERNAL_MAGIC;
49  VMOVE(eip->v, sip->v);
50  VMOVE(eip->a, sip->a);
51  VMOVE(eip->b, sip->b);
52  VMOVE(eip->c, sip->c);
53 
54  struct rt_db_internal tmp_internal;
55  RT_DB_INTERNAL_INIT(&tmp_internal);
56  tmp_internal.idb_major_type = DB5_MAJORTYPE_BRLCAD;
57  tmp_internal.idb_ptr = (void *)eip;
58  tmp_internal.idb_minor_type = ID_ELL;
59  tmp_internal.idb_meth = &OBJ[ID_ELL];
60 
61  ON_Brep *tmp_brep = ON_Brep::New();
62  rt_ell_brep(&tmp_brep, &tmp_internal, tol);
63 
64  ON_NurbsSurface *surf = ON_NurbsSurface::New();
65  tmp_brep->m_S[0]->GetNurbForm(*surf);
66 
67  // Calculate the weight value
68  // See: Joe F. Thompson, B. K. Soni, N. P. Weatherill. Handbook of Grid Generation.
69  // Ch.30.3.5 Superellipse to NURBS Curve.
70  const fastf_t cos45 = cos(ON_PI/4.0);
71  fastf_t weight_e = (pow(cos45, sip->e) - 0.5) / (1 - pow(cos45, sip->e));
72  fastf_t weight_n = (pow(cos45, sip->n) - 0.5) / (1 - pow(cos45, sip->n));
73  // When the center (origin) is used as a control point
74  if (weight_e < 0) {
75  weight_e = (pow(cos45, sip->e) - 0.5) / pow(cos45, sip->e);
76  }
77  if (weight_n < 0) {
78  weight_n = (pow(cos45, sip->n) - 0.5) / pow(cos45, sip->n);
79  }
80 
81  // Change the weight of the control points. rt_ell_brep() creates a NURBS surface with 9*5
82  // control points.
83  // When both e and n <= 2, this can generate an ideal b-rep for superell, otherwise cannot.
84  if (sip->e <= 2 && sip->n <= 2) {
85  for (int i = 0; i < surf->CVCount(0); i++) {
86  fastf_t tmp_weight = 0.0;
87  switch (i) {
88  case 0:
89  case 4:
90  case 8: // the XZ plane
91  case 2:
92  case 6: // the YZ plane
93  tmp_weight = 1.0;
94  break;
95  case 1:
96  case 3:
97  case 5:
98  case 7:
99  tmp_weight = weight_e;
100  break;
101  default:
102  bu_log("Should not reach here!\n");
103  }
104  for (int j = 0; j < surf->CVCount(1); j++) {
105  fastf_t new_weight = 0.0;
106  switch (j) {
107  case 0:
108  case 2:
109  case 4:
110  new_weight = tmp_weight;
111  break;
112  case 1:
113  case 3:
114  new_weight = tmp_weight * weight_n;
115  break;
116  default:
117  bu_log("Should not reach here!\n");
118  }
119  ON_3dPoint ctrlpt;
120  surf->GetCV(i, j, ctrlpt);
121  ON_4dPoint newctrlpt(ctrlpt[0]*new_weight, ctrlpt[1]*new_weight, ctrlpt[2]*new_weight, new_weight);
122  surf->SetCV(i, j, newctrlpt);
123  }
124  }
125 
126  surf->SetDomain(0, 0.0, 1.0);
127  surf->SetDomain(1, 0.0, 1.0);
128 
129  // Make final BREP structure
130  (*b)->m_S.Append(surf);
131  int surfindex = (*b)->m_S.Count();
132  (*b)->NewFace(surfindex - 1);
133  int faceindex = (*b)->m_F.Count();
134  (*b)->NewOuterLoop(faceindex-1);
135  } else {
136  // When one or more than one of e and n are greater than 2, generate 8 surfaces
137  // independently to represent the superell
138 
139  // First, we create a b-rep superell centered at the origin, with radius 1.
140  ON_NurbsSurface *surface[8];
141 
142  surface[0] = ON_NurbsSurface::New(3, true, 3, 3, 3, 3);
143  surface[0]->SetKnot(0, 0, 0);
144  surface[0]->SetKnot(0, 1, 0);
145  surface[0]->SetKnot(0, 2, 1);
146  surface[0]->SetKnot(0, 3, 1);
147  surface[0]->SetKnot(1, 0, 0);
148  surface[0]->SetKnot(1, 1, 0);
149  surface[0]->SetKnot(1, 2, 1);
150  surface[0]->SetKnot(1, 3, 1);
151 
152  // Surface 1
153  surface[0]->SetCV(0, 0, ON_4dPoint(0, 0, -1, 1));
154  surface[0]->SetCV(0, 2, ON_4dPoint(1, 0, 0, 1));
155  if (weight_e > 0) {
156  surface[0]->SetCV(1, 0, ON_4dPoint(0, 0, -weight_e, weight_e));
157  surface[0]->SetCV(1, 2, ON_4dPoint(weight_e, weight_e, 0, weight_e));
158  } else {
159  surface[0]->SetCV(1, 0, ON_4dPoint(0, 0, weight_e, -weight_e));
160  surface[0]->SetCV(1, 2, ON_4dPoint(0, 0, 0, -weight_e));
161  }
162  surface[0]->SetCV(2, 0, ON_4dPoint(0, 0, -1, 1));
163  surface[0]->SetCV(2, 2, ON_4dPoint(0, 1, 0, 1));
164 
165  if (weight_n >= 0) {
166  surface[0]->SetCV(0, 1, ON_4dPoint(weight_n, 0, -weight_n, weight_n));
167  surface[0]->SetCV(1, 1, ON_4dPoint(0, 0, -fabs(weight_n*weight_e), fabs(weight_n*weight_e)));
168  surface[0]->SetCV(2, 1, ON_4dPoint(0, weight_n, -weight_n, weight_n));
169  } else {
170  // use the center as a control point
171  surface[0]->SetCV(0, 1, ON_4dPoint(0, 0, 0, -weight_n));
172  if (weight_e > 0) {
173  surface[0]->SetCV(1, 1, ON_4dPoint(0, 0, 0, -weight_n*weight_e));
174  } else {
175  surface[0]->SetCV(1, 1, ON_4dPoint(0, 0, 0, -weight_n));
176  }
177  surface[0]->SetCV(2, 1, ON_4dPoint(0, 0, 0, -weight_n));
178  }
179 
180  surface[0]->SetDomain(0, 0.0, 1.0);
181  surface[0]->SetDomain(1, 0.0, 1.0);
182 
183  ON_3dPoint axis_a(0, 0, -1);
184  ON_3dPoint axis_b(0, 0, 0);
185  ON_3dVector axis;
186  axis = axis_a - axis_b;
187 
188  // Surface 2
189  surface[1] = surface[0]->Duplicate();
190  surface[1]->Rotate(ON_PI/2.0, axis, axis_b);
191 
192  // Surface 3
193  surface[2] = surface[0]->Duplicate();
194  surface[2]->Rotate(ON_PI, axis, axis_b);
195 
196  // Surface 4
197  surface[3] = surface[0]->Duplicate();
198  surface[3]->Rotate(ON_PI*3.0/2.0, axis, axis_b);
199 
200  axis_a = ON_3dPoint(0, 1, 0);
201  ON_3dVector axis2;
202  axis2 = axis_a - axis_b;
203 
204  // Surface 5
205  surface[4] = surface[0]->Duplicate();
206  surface[4]->Rotate(ON_PI, axis2, axis_b);
207 
208  // Surface 6
209  surface[5] = surface[4]->Duplicate();
210  surface[5]->Rotate(ON_PI/2.0, axis, axis_b);
211 
212  // Surface 7
213  surface[6] = surface[4]->Duplicate();
214  surface[6]->Rotate(ON_PI, axis, axis_b);
215 
216  // Surface 8
217  surface[7] = surface[4]->Duplicate();
218  surface[7]->Rotate(ON_PI*3.0/2.0, axis, axis_b);
219 
220  // Do the transformations and add elements to the b-rep object
221  for (int i = 0; i < 8; i++) {
222  mat_t R, invR, invS, invRinvS;
223  MAT_IDN(R);
224  VMOVE(&R[0], sip->a);
225  VUNITIZE(&R[0]);
226  VMOVE(&R[4], sip->b);
227  VUNITIZE(&R[4]);
228  VMOVE(&R[8], sip->c);
229  VUNITIZE(&R[8]);
230  bn_mat_trn(invR, R);
231 
232  MAT_IDN(invS);
233  invS[0] = MAGNITUDE(sip->a);
234  invS[5] = MAGNITUDE(sip->b);
235  invS[10] = MAGNITUDE(sip->c);
236  bn_mat_mul(invRinvS, invR, invS);
237 
238  for (int j = 0; j < surface[j]->CVCount(0); j++) {
239  for (int k = 0; k < surface[k]->CVCount(1); k++) {
240  point_t cvpt;
241  ON_4dPoint ctrlpt;
242  surface[i]->GetCV(j, k, ctrlpt);
243  MAT3X3VEC(cvpt, invRinvS, ctrlpt);
244  point_t scale_v;
245  VSCALE(scale_v, eip->v, ctrlpt.w);
246  VADD2(cvpt, scale_v, cvpt);
247  ON_4dPoint newpt = ON_4dPoint(cvpt[0], cvpt[1], cvpt[2], ctrlpt.w);
248  surface[i]->SetCV(j, k, newpt);
249  }
250  }
251 
252  (*b)->m_S.Append(surface[i]);
253  int surfindex = (*b)->m_S.Count();
254  (*b)->NewFace(surfindex - 1);
255  int faceindex = (*b)->m_F.Count();
256  (*b)->NewOuterLoop(faceindex-1);
257  }
258  }
259 }
260 
261 
262 // Local Variables:
263 // tab-width: 8
264 // mode: C++
265 // c-basic-offset: 4
266 // indent-tabs-mode: t
267 // c-file-style: "stroustrup"
268 // End:
269 // ex: shiftwidth=4 tabstop=8
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
void rt_superell_brep(ON_Brep **b, const struct rt_db_internal *ip, const struct bn_tol *tol)
Header file for the BRL-CAD common definitions.
int idb_major_type
Definition: raytrace.h:192
#define RT_CK_DB_INTERNAL(_p)
Definition: raytrace.h:207
#define RT_ELL_INTERNAL_MAGIC
Definition: magic.h:91
#define BU_ALLOC(_ptr, _type)
Definition: malloc.h:223
#define RT_DB_INTERNAL_INIT(_p)
Definition: raytrace.h:199
const struct rt_functab * idb_meth
for ft_ifree(), etc.
Definition: raytrace.h:194
void bn_mat_trn(mat_t om, const mat_t im)
void bn_mat_mul(mat_t o, const mat_t a, const mat_t b)
Support for uniform tolerances.
Definition: tol.h:71
void rt_ell_brep(ON_Brep **b, struct rt_db_internal *ip, const struct bn_tol *tol)
void * idb_ptr
Definition: raytrace.h:195
const struct rt_functab OBJ[]
Definition: table.c:159
axis
Definition: color.c:48
#define R
Definition: msr.c:55
int idb_minor_type
ID_xxx.
Definition: raytrace.h:193
double fastf_t
Definition: defines.h:300
#define ID_ELL
Ellipsoid.
Definition: raytrace.h:461