1 /* Copyright (C) 1992-1998 The Geometry Center
2 * Copyright (C) 1998-2000 Stuart Levy, Tamara Munzner, Mark Phillips
3 *
4 * This file is part of Geomview.
5 *
6 * Geomview is free software; you can redistribute it and/or modify it
7 * under the terms of the GNU Lesser General Public License as published
8 * by the Free Software Foundation; either version 2, or (at your option)
9 * any later version.
10 *
11 * Geomview 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 Geomview; see the file COPYING. If not, write
18 * to the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139,
19 * USA, or visit http://www.gnu.org.
20 */
21
22 #if 0
23 static char copyright[] = "Copyright (C) 1992-1998 The Geometry Center\n\
24 Copyright (C) 1998-2000 Stuart Levy, Tamara Munzner, Mark Phillips";
25 #endif
26
27 #include "geom.h"
28 #include "list.h"
29 #include "plutil.h"
30 #include "Clip.h" /* Note: must include this *after* oogl includes */
31
32
33 Clip clip;
34
35 extern int span_vertices(Clip *clip, float *minp, float *maxp);
36
37 static char Usage[] = "\
38 Usage: clip [-v axisx,y,z,...] [-g value] [-l value] [-s nslices[,fraction]]\n\
39 [-sph centerx,y,z,...] [-cyl centerx,y,z,...] [-e] [file.oogl]\n\
40 Reads an OOGL object from file.oogl (or stdin if omitted).\n\
41 Slices it against a (series of) planes whose normal vectors are given\n\
42 by the -v argument (default 1,0,0,0,...).\n\
43 Given -g, selects portions where <point> dot <vector> > <pvalue>.\n\
44 Given -l, selects portions where <point> dot <vector> < <nvalue>.\n\
45 Given both, takes the portion lying between those two values.\n\
46 With -s, emits a series of slices spaced <spacing> units apart.\n\
47 With -e, emits just two numbers instead of geometry: range of function-values\n\
48 which the object spans.\n\
49 -sph : Slice against a sphere, not a plane. -g/-l gives radius.\n\
50 -cyl x,y,z... : Slice against a cylinder whose axis passes through x,y,z...\n\
51 -v gives axis direction, -g/-l radius.\n";
52
parsevector(char * str,float data[],int maxdata)53 static int parsevector(char *str, float data[], int maxdata)
54 {
55 int i;
56 char *p, *op;
57
58 memset(data, 0, maxdata*sizeof(float));
59 if(str == NULL)
60 return -1;
61 for(i = 0, p = str; i < maxdata && *p; i++) {
62 op = p;
63 data[i] = strtod(p, &p);
64 if(p == op)
65 break;
66 while(*p == ',' || *p == ' ') p++;
67 }
68 return i;
69 }
70
71 void
setclipat(Clip * clip,char * pstr,int dim,float * surf,void (* prepfunc)())72 setclipat(Clip *clip, char *pstr, int dim, float *surf, void (*prepfunc)())
73 {
74 float point[MAXDIM];
75 float level;
76 int n;
77
78 n = parsevector(pstr, point, dim);
79
80 if(n > 1) {
81 /* Given several numbers, consider them
82 * to be a point; choose a level-set
83 * going through that point.
84 */
85 setClipPlane(clip, surf, 0);
86 level = (*clip->clipfunc)(clip, point);
87 } else {
88 level = point[0]; /* If just one number, use exactly it */
89 }
90 setClipPlane(clip, surf, level);
91 if(prepfunc)
92 (*prepfunc)(clip);
93 }
94
95 /*
96 * Clipping functions.
97 * The default function, to clip against a plane, is in clip.c.
98 */
99 static float origin[MAXDIM];
100
101 /*
102 * Clip against a sphere.
103 * Uses static variable "origin".
104 */
sphere(Clip * clip,float * point)105 static float sphere(Clip *clip, float *point)
106 {
107 /* clip->surf[0..dim-1] are the coordinates of the center of the sphere. */
108 int i;
109 float r2 = 0, *s, *p;
110 for(i = clip->dim, p = point, s = origin; --i >= 0; p++, s++)
111 r2 += (*p - *s) * (*p - *s);
112 return r2;
113 }
114
115
116 /*
117 * Clip against a cylinder.
118 * Uses both clip->surf, the cylinder's axis vector,
119 * and static variable "origin", a point on the cylinder.
120 */
121
cylinder(Clip * clip,float * point)122 static float cylinder(Clip *clip, float *point)
123 {
124 /* clip->surf[0..dim-1] are the direction vector of the cylinder's axis;
125 * clip->surf[dim..2*dim-1] are coordinates of a point on the axis.
126 */
127 int dim = clip->dim;
128 int i;
129 float r2 = 0, dot = 0, dx, *s, *p, *c;
130
131 for(i=dim, p=point, c=origin, s=clip->surf; --i >= 0; p++, s++, c++)
132 dot += (*p - *c) * *s;
133
134 for(i=dim, p=point, c=origin, s=clip->surf; --i >= 0; p++, s++, c++) {
135 dx = (*p - *c) - dot * *s;
136 r2 += dx * dx;
137 }
138 return r2;
139 }
140
141
142 /* Special case for Mark Levi */
leviarctan(Clip * clip,float * point)143 float leviarctan(Clip *clip, float *point)
144 {
145 /* On input, point[0] .. point[2] are x,y,z coordinates of the point to test.
146 * Assigning w <- point[1] and T <- point[2]
147 * we compute x(w,T) on the surface
148 * x = constant - 2T +/- arctan(w)
149 * and return the difference between x(w,t) and our point's x.
150 *
151 * The sign of arctan is taken from clip->surf[0]; this is
152 * the number given to the "-v" command-line argument. So, using
153 * "-a -v 1,21.99" means x(w,T) = 7pi - 2T - arctan(w), and
154 * "-a -v -1,21.99" means x(w,T) = 7pi - 2T + arctan(w).
155 */
156
157 return clip->surf[1] - 2*point[2] - clip->surf[0] * atan(point[1])
158 - point[0];
159 }
160
161 /*
162 * For some level-set functions we'd rather precompute the square of the
163 * given value. They set prepfunc = squared.
164 */
squared(Clip * clip)165 static void squared(Clip *clip)
166 {
167 clip->level *= clip->level;
168 }
169
170 /* For the cylinder, we want to make the axis be a unit vector,
171 * so we can easily project it out. Want to square the level value, too.
172 */
squared_normalized(Clip * clip)173 static void squared_normalized(Clip *clip)
174 {
175 int dim = clip->dim;
176 int i;
177 float len = 0, *p;
178 for(i = dim, p = clip->surf; --i >= 0; p++)
179 len += *p * *p;
180 len = sqrt(len);
181 if(len > 0) {
182 for(i = dim, p = clip->surf; --i >= 0; p++)
183 *p /= len;
184 }
185
186 /* And square the radius */
187 clip->level *= clip->level;
188 }
189
190
main(int argc,char * argv[])191 int main(int argc, char *argv[])
192 {
193 Geom *g, *clipped;
194 float surf[MAXDIM];
195 float (*func)() = NULL;
196 void (*prepfunc)() = NULL;
197 char *lestr = NULL, *gestr = NULL;
198 char *type;
199 int slices = 0;
200 int nonlinear = 0;
201 int extent = 0;
202 float slicev[2];
203 IOBFILE *inf = iobfileopen(stdin);
204
205
206 if(argc <= 1) {
207 fputs(Usage, stderr);
208 exit(1);
209 }
210
211 clip_init(&clip);
212
213 memset(surf, 0, sizeof(surf));
214 surf[0] = 1;
215
216 while(argc > 1) {
217 if(strncmp(argv[1], "-g", 2) == 0 && argc>2) {
218 if(gestr != NULL) {
219 fprintf(stderr, "clip: can only specify one -g option.\n");
220 exit(1);
221 }
222 gestr = argv[2];
223 argc -= 2, argv += 2;
224 } else if(strncmp(argv[1], "-l", 2) == 0 && argc>2) {
225 if(lestr != NULL) {
226 fprintf(stderr, "clip: can only specify one -l option.\n");
227 exit(1);
228 }
229 lestr = argv[2];
230 argc -= 2, argv += 2;
231
232 } else if(strcmp(argv[1], "-v") == 0 && argc > 2) {
233 if(parsevector(argv[2], surf, MAXDIM) <= 0) {
234 fprintf(stderr, "-v: expected x,y,z,... surface parameters\n\
235 (plane normal or cylinder axis)\n");
236 exit(1);
237 }
238 argc -= 2, argv += 2;
239
240 } else if(strcmp(argv[1], "-e") == 0) {
241 extent = 1;
242 argc--, argv++;
243
244 } else if(strcmp(argv[1], "-a") == 0) {
245 func = leviarctan;
246 nonlinear = 1;
247 argc--, argv++;
248
249 } else if(strncmp(argv[1], "-sph", 4) == 0) {
250 func = sphere;
251 prepfunc = squared;
252 nonlinear = 1;
253 if(parsevector(argv[2], origin, COUNT(origin)) <= 0) {
254 fprintf(stderr, "clip -sph: expected x,y,z,... of sphere's center.\n");
255 exit(1);
256 }
257 argc -= 2, argv += 2;
258
259 } else if(strncmp(argv[1], "-cyl", 4) == 0) {
260 func = cylinder;
261 prepfunc = squared_normalized;
262 nonlinear = 1;
263 if(parsevector(argv[2], origin, COUNT(origin)) <= 0) {
264 fprintf(stderr, "clip -cyl: expected x,y,z,... of a point on the cylinder's axis\n");
265 exit(1);
266 }
267 argc -= 2, argv += 2;
268
269 } else if(strcmp(argv[1], "-s") == 0) {
270 slices = parsevector(argv[2], slicev, 2);
271 if(slices <= 0 || slicev[0] <= 1) {
272 fprintf(stderr, "clip -s: expected -s nribbons or nribbons,gap-fraction; must have >= 2 ribbons\n");
273 exit(1);
274 } else if(slices == 1)
275 slicev[1] = .5;
276 argc -= 2, argv += 2;
277
278 } else if(argv[1][0] == '-') {
279 fprintf(stderr, "Unknown option '%s'. Run 'clip' with no arguments for help.\n",
280 argv[1]);
281 exit(1);
282 } else {
283 break;
284 }
285 }
286
287 if(lestr == NULL && gestr == NULL && slices == 0 && extent == 0) {
288 fprintf(stderr, "clip: must specify at least one of -l or -g or -s or -e options.\n");
289 exit(2);
290 }
291
292 if(argc > 2) {
293 fprintf(stderr, "clip: can only handle one object\n");
294 exit(2);
295 }
296
297 if(argc > 1 && strcmp(argv[1], "-") != 0) {
298 g = GeomLoad(argv[1]);
299 } else {
300 argv[1] = "standard input";
301 g = GeomFLoad(inf, "stdin");
302 }
303
304 type = GeomName(g);
305 if(type!=NULL && strcmp(type, "polylist") && strcmp(type, "npolylist")) {
306 Geom *newg = AnyToPL(g, TM3_IDENTITY);
307 GeomDelete(g);
308 g = newg;
309 }
310
311 if(g == NULL) {
312 fprintf(stderr, "clip: Error loading OOGL object from %s\n", argv[1]);
313 exit(1);
314 }
315
316 if(func != NULL)
317 clip.clipfunc = func;
318 clip.nonlinear = nonlinear;
319
320 if(lestr != NULL) {
321 setGeom(&clip, g);
322 setclipat(&clip, lestr, clip.dim, surf, prepfunc);
323 setSide(&clip, CLIP_LE);
324 do_clip(&clip);
325 clipped = getGeom(&clip);
326 GeomDelete(g);
327 g = clipped;
328 }
329
330 if(gestr != NULL) {
331 setGeom(&clip, g);
332 setclipat(&clip, gestr, clip.dim, surf, prepfunc);
333 setSide(&clip, CLIP_GE);
334 do_clip(&clip);
335 clipped = getGeom(&clip);
336 GeomDelete(g);
337 g = clipped;
338 }
339
340
341 if(slices > 0) {
342 float min, max, v, step;
343 int i;
344 Geom *piece, *whole = NULL;
345 char lim[64];
346
347 setGeom(&clip, g);
348 setclipat(&clip, "0", clip.dim, surf, prepfunc);
349 if(span_vertices(&clip, &min, &max) != 0) {
350 /* Nonempty */
351 /* Interpolate linearly in real space -- invert prepfunc */
352 if(prepfunc == squared || prepfunc == squared_normalized) {
353 min = sqrt(min);
354 max = sqrt(max);
355 }
356 step = (max - min) / (slicev[0] + slicev[1] - 1);
357 sprintf(lim, "%g", min + slicev[1]*step);
358 setclipat(&clip, lim, clip.dim, surf, prepfunc);
359 setSide(&clip, CLIP_LE);
360 do_clip(&clip);
361 whole = GeomCreate("list", CR_GEOM, getGeom(&clip), CR_END);
362 for(i = 1, v = min + step; i < slicev[0]; i++, v += step) {
363 setGeom(&clip, g);
364 sprintf(lim, "%g", v + slicev[1]*step);
365 setclipat(&clip, lim, clip.dim, surf, prepfunc);
366 setSide(&clip, CLIP_LE);
367 do_clip(&clip);
368 piece = getGeom(&clip);
369
370 setGeom(&clip, piece);
371 sprintf(lim, "%g", v);
372 setclipat(&clip, lim, clip.dim, surf, prepfunc);
373 setSide(&clip, CLIP_GE);
374 do_clip(&clip);
375 piece = getGeom(&clip);
376 if(piece)
377 whole = ListAppend(whole, piece);
378 }
379 }
380 g = whole;
381 }
382
383 if(extent != 0) {
384 float min, max;
385 setGeom(&clip, g);
386 setclipat(&clip, "0", clip.dim, surf, prepfunc);
387 if(span_vertices(&clip, &min, &max) != 0) {
388 if(prepfunc == squared || prepfunc == squared_normalized) {
389 min = sqrt(min);
390 max = sqrt(max);
391 }
392 }
393 printf("%g %g\n", min, max);
394 } else {
395 GeomFSave(g, stdout, "stdout");
396 }
397 exit(0);
398 }
399