summaryrefslogtreecommitdiff
path: root/graphics/sketch/geometry.h
blob: c7863fd836cf748f591679dd5cfd64cea0c304b0 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
/* geometry.h
   Copyright (C) 2005,2006,2007 Eugene K. Ressler, Jr.

This file is part of Sketch, a small, simple system for making 
3d drawings with LaTeX and the PSTricks or TikZ package.

Sketch is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 3, or (at your option)
any later version.

Sketch is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with Sketch; see the file COPYING.txt.  If not, see
http://www.gnu.org/copyleft */

#ifndef __GEOMETRY_H
#define __GEOMETRY_H

// ---- memory -----------------------------------------------------------------

#include <float.h>		// floating point definitions
#include "dynarray.h"		// dynamic arrays

// bit N
#define bit(N)  (1u << (N))

// size of a static or auto declared array
#define ARRAY_SIZE(A) (sizeof (A) / sizeof (A)[0])

// checking memory allocators
void *safe_malloc (unsigned size);
void *safe_realloc (void *p, unsigned size);
char *safe_strdup (char *str);
void safe_free (void *p);
#ifdef _DEBUG
#define malloc(N)    __call_safe_malloc_instead()
#define realloc(P,N) __call_safe_alloc_instead()
#define strdup(S)    __call_safe_alloc_instead()
#define free(P)      __call_safe_free_instead()
#endif

// ---- numerics ---------------------------------------------------------------

// float declarations to ease compilation 
// with either single or double precision
typedef unsigned int SIZE, INDEX;
typedef double FLOAT;

#define FLOAT_SCAN_FMT "%lf"
#define FLOAT_EPS (8*DBL_EPSILON)
#define FLOAT_MIN FLT_MIN
#define FLOAT_MAX FLT_MAX

#ifdef _MSC_VER
// kill loss of precision warnings for case where FLOAT is float
#pragma warning(disable:4244 4305)
#endif

#define PI ((FLOAT)3.1415926535897932384626433832795028841971693993751)

// Max and min operators
FLOAT max_float (FLOAT x, FLOAT y);
FLOAT min_float (FLOAT x, FLOAT y);

// ---- points -----------------------------------------------------------------

// indices
#define X 0
#define Y 1
#define Z 2
#define W 3

// points
typedef FLOAT POINT_2D[2], POINT_3D[3];
void copy_pt_2d (POINT_2D r, POINT_2D s);
void copy_pt_3d (POINT_3D r, POINT_3D s);
void find_pt_3d_from_2d (POINT_3D r, POINT_2D pt);

// ---- polylines --------------------------------------------------------------

// polylines are just dynamic arrays of points

typedef struct polyline_2d_t
{
  DYNAMIC_2D_ARRAY_FIELDS (POINT_2D, v, n_vertices);
  struct polyline_2d_t *next;
}
POLYLINE_2D;

DECLARE_DYNAMIC_2D_ARRAY_PROTOS (POLYLINE_2D, POINT_2D, FLOAT, polyline_2d,
				 v, n_vertices)
     typedef struct polyline_3d_t
     {
       DYNAMIC_2D_ARRAY_FIELDS (POINT_3D, v, n_vertices);
       struct polyline_3d_t *next;
     }
POLYLINE_3D;

DECLARE_DYNAMIC_2D_ARRAY_PROTOS (POLYLINE_3D, POINT_3D, FLOAT, polyline_3d,
				 v, n_vertices)
// ---- polygons ---------------------------------------------------------------
// polygons are just a dynamic arrays of points; chains represent complex polygons
     typedef struct polygon_2d_t
     {
       DYNAMIC_2D_ARRAY_FIELDS (POINT_2D, v, n_sides);
       struct polygon_2d_t *next;
     }
POLYGON_2D;

DECLARE_DYNAMIC_2D_ARRAY_PROTOS (POLYGON_2D, POINT_2D, FLOAT, polygon_2d, v,
				 n_sides)
     typedef struct polygon_3d_t
     {
       DYNAMIC_2D_ARRAY_FIELDS (POINT_3D, v, n_sides);
       struct polygon_3d_t *next;
     }
POLYGON_3D;

DECLARE_DYNAMIC_2D_ARRAY_PROTOS (POLYGON_3D, POINT_3D, FLOAT, polygon_3d, v,
				 n_sides)
// ---- vectors ----------------------------------------------------------------
     typedef FLOAT *VECTOR;

// vectors of dynamic length
     void init_vec (VECTOR * v);
     void clear_vec (VECTOR * v);
     void setup_vec (VECTOR * v, SIZE n);
     void init_and_setup_vec (VECTOR * v, SIZE n);
     void zero_vec (VECTOR r, SIZE n);
     void copy_vec (VECTOR r, VECTOR v, SIZE n);

// vectors of useful static length.
     typedef FLOAT VECTOR_2D[2], VECTOR_3D[3], VECTOR_4D[4];

     FLOAT length_vec_2d (VECTOR_2D v);
     FLOAT length_vec_3d (VECTOR_3D v);
     FLOAT dist_2d (POINT_2D p1, POINT_2D p2);
     FLOAT dist_3d (POINT_3D p1, POINT_3D p2);
     FLOAT length_vec_2d_sqr (VECTOR_2D v);
     FLOAT length_vec_3d_sqr (VECTOR_3D v);
     FLOAT dist_2d_sqr (POINT_2D p1, POINT_2D p2);
     FLOAT dist_3d_sqr (POINT_3D p1, POINT_3D p2);
     void zero_vec_2d (VECTOR_2D v);
     void zero_vec_3d (VECTOR_3D v);
     void negate_vec_2d (VECTOR_2D r, VECTOR_2D v);
     void negate_vec_3d (VECTOR_3D r, VECTOR_3D v);
     void copy_vec_2d (VECTOR_2D r, VECTOR_2D s);
     void copy_vec_3d (VECTOR_3D r, VECTOR_3D s);
     void scale_vec_2d (VECTOR_2D r, VECTOR_2D v, FLOAT s);
     void scale_vec_3d (VECTOR_3D r, VECTOR_3D v, FLOAT s);
     int find_unit_vec_2d (VECTOR_2D r, VECTOR_2D v);
     int find_unit_vec_3d (VECTOR_3D r, VECTOR_3D v);
     void add_vecs_2d (VECTOR_2D r, VECTOR_2D a, VECTOR_2D b);
     void add_vecs_3d (VECTOR_3D r, VECTOR_3D a, VECTOR_3D b);
     void sub_vecs_2d (VECTOR_2D r, VECTOR_2D a, VECTOR_2D b);
     void sub_vecs_3d (VECTOR_3D r, VECTOR_3D a, VECTOR_3D b);
     void add_vec_to_pt_2d (POINT_2D r, POINT_2D pt, VECTOR_2D v);
     void add_vec_to_pt_3d (POINT_3D r, POINT_3D pt, VECTOR_3D v);
     void add_scaled_vec_to_pt_2d (POINT_2D r, POINT_2D pt, VECTOR_2D v,
				   FLOAT s);
     void add_scaled_vec_to_pt_3d (POINT_3D r, POINT_3D pt, VECTOR_3D v,
				   FLOAT s);
     void sub_pts_2d (VECTOR_2D r, POINT_2D a, POINT_2D b);
     void sub_pts_3d (VECTOR_3D r, POINT_3D a, POINT_3D b);
     void fold_min_pt_2d (POINT_2D min, POINT_2D new_pt);
     void fold_min_pt_3d (POINT_3D min, POINT_3D new_pt);
     void fold_max_pt_2d (POINT_2D max, POINT_2D new_pt);
     void fold_max_pt_3d (POINT_3D max, POINT_3D new_pt);

     FLOAT dot_2d (VECTOR_2D a, VECTOR_2D b);
     FLOAT dot_3d (VECTOR_3D a, VECTOR_3D b);
     void cross (VECTOR_3D r, VECTOR_3D a, VECTOR_3D b);

// linear interpolation operators
     void lerp_2d (POINT_2D r, FLOAT t, POINT_2D p1, POINT_2D p2);
     void lerp_3d (POINT_3D r, FLOAT t, POINT_3D p1, POINT_3D p2);

// find parameters of intersection point of two line segments
     int line_intersect_2d (POINT_2D a, POINT_2D b, POINT_2D c, POINT_2D d,
			    FLOAT eps, FLOAT * t_ab, FLOAT * t_cd);

// ---- planes -----------------------------------------------------------------
     typedef struct plane_t
     {
       VECTOR_3D n;
       POINT_3D p;
       FLOAT c;
     }
PLANE;

// return description of the plane of a polygon using Newell's method
     void find_polygon_plane (PLANE * plane, POLYGON_3D * polygon);

#define S_IN      (1)
#define S_ON      (2)
#define S_OUT     (4)
#define S_IN_ON   (S_ON | 8)
#define S_OUT_ON  (S_ON | 16)
#define S_SPLIT   (32)

// #define PLANE_HALF_THICKNESS  (10.0 * FLOAT_EPS)
#define PLANE_HALF_THICKNESS  (.001/2)

// given a plane of thickness 2 * half_thickness, return:
//   S_IN or S_OUT if the point is resp. inside or outside the thickness of the plane
//   S_IN_ON or S_OUT_ON if the point is within half_thickness of the plane on the resp. side
//   S_ON if the point is precisely on the plane; no IN or OUT determination can be made
     int pt_side_of_plane (PLANE * plane, POINT_3D p);

// given a polygon and a plane, return:
//   S_IN if all the verices are IN or ON the thickened plane
//   S_OUT if all the verices are OUTside or ON the thickened plane
//   S_ON if all vertice are ON the thickened plane
//   S_SPLIT otherwise
     int polygon_side_of_plane (POLYGON_3D * polygon, PLANE * plane);

// given a polyline and a plane, return:
//   S_IN if all segments of the line are fully INside the thickened plane
//   S_OUT if all segments of the line are fully OUTside the thickened plane
//   S_ON if all vertice are ON the thickened plane
//   S_SPLIT otherwise
     int polyline_side_of_plane (POLYLINE_3D * polyline, PLANE * plane);

// ---- boxes ------------------------------------------------------------------

     typedef struct box_2d_t
     {
       POINT_2D min, max;
     }
BOX_2D;

     typedef struct box_3d_t
     {
       POINT_3D min, max;
     }
BOX_3D;

     void init_box_2d (BOX_2D * b);
     void init_box_3d (BOX_3D * b);
     void fold_min_max_pt_2d (BOX_2D * b, POINT_2D p);
     void fold_min_max_pt_3d (BOX_3D * b, POINT_3D p);
     void fold_min_max_polygon_2d (BOX_2D * b, POLYGON_2D * polygon);
     void fold_min_max_polygon_3d (BOX_3D * b, POLYGON_3D * polygon);
     void fold_min_max_polyline_2d (BOX_2D * b, POLYLINE_2D * polyline);
     void fold_min_max_polyline_3d (BOX_3D * b, POLYLINE_3D * polyline);
     void copy_box_2d (BOX_2D * r, BOX_2D * s);
     void copy_box_3d (BOX_3D * r, BOX_3D * s);
     int boxes_2d_intersect_p (BOX_2D * a, BOX_2D * b);
     int boxes_3d_intersect_p (BOX_2D * a, BOX_2D * b);

// ---- transformations --------------------------------------------------------

// homogeneous transform stored in column major order
     typedef FLOAT TRANSFORM[16];

// for initializations of identity transforms
#define IDENT_TRANSFORM  \
{ 1.0, 0.0, 0.0, 0.0, \
  0.0, 1.0, 0.0, 0.0, \
  0.0, 0.0, 1.0, 0.0, \
  0.0, 0.0, 0.0, 1.0 }

// ---- global contstants ------------------------------------------------------

     extern TRANSFORM identity;
     extern POINT_2D origin_2d;
     extern POINT_3D origin_3d;
     extern VECTOR_2D I_2d;
     extern VECTOR_2D J_2d;
     extern VECTOR_3D I_3d;
     extern VECTOR_3D J_3d;
     extern VECTOR_3D K_3d;

// row-column tranform indexing matches OpenGL convention: column major
#define IT(I,J) (4 * ((J) - 1) + ((I) - 1))

// copy source to result transform
     void copy_transform (TRANSFORM r, TRANSFORM s);

// set the result transform to the identity
     void set_ident (TRANSFORM r);

// create a rotation transform thru angle theta about axis u (must be unit vec)
     void set_angle_axis_rot (TRANSFORM r, FLOAT theta, VECTOR_3D u);

// create a rotation transform thru angle theta
// u is optional axis which need not be a unit vector (default is [0,0,1])
// p is optional center of rotation (default is (0,0,0))
     void set_angle_axis_rot_about_point (TRANSFORM r, FLOAT theta,
					  POINT_3D p, VECTOR_3D u);

// create a scale transform
     void set_scale (TRANSFORM r, FLOAT sx, FLOAT sy, FLOAT sz);

// create a translation transform
     void set_translation (TRANSFORM r, FLOAT dx, FLOAT dy, FLOAT dz);

// create a true perspective projection (depth = p for all projected points)
     void set_perspective_projection (TRANSFORM r, FLOAT p);

// create a perspective transformation (depth is a pseudodepth)
     void set_perspective_transform (TRANSFORM r, FLOAT p);

// create a true parallel projection (depth = 0 for all projected points)
     void set_parallel_projection (TRANSFORM r);

// create an OpenGL-like view transformation matrix
     void set_view_transform (TRANSFORM r, POINT_3D eye, VECTOR_3D vd,
			      VECTOR_3D up);
     void set_view_transform_with_look_at (TRANSFORM r, POINT_3D eye,
					   POINT_3D look_at, VECTOR_3D up);

// invert a given transform m; return its determinant; we give up if the 
// determinant is too small
     void invert (TRANSFORM r, FLOAT * det_rtn, TRANSFORM m, FLOAT min_det);

// compose two transforms, but result cannot be the same as either operand
     void compose_unsafe (TRANSFORM r, TRANSFORM a, TRANSFORM b);

// same as above, but safe to use either operand to hold result.
     void compose (TRANSFORM r, TRANSFORM a, TRANSFORM b);

     void transform_pt_3d (POINT_3D r, TRANSFORM m, POINT_3D p);
     void transform_vec_3d (VECTOR_3D r, TRANSFORM m, VECTOR_3D p);

// ---- quaternions ------------------------------------------------------------

     typedef FLOAT QUATERNION[4];

// for initializations of identity quaternions
#define IDENT_QUAT { 0.0, 0.0, 0.0, 1.0 }

     void set_ident_quat (QUATERNION q);
     void set_angle_axis_quat (QUATERNION q, FLOAT theta, VECTOR_3D axis);
     void find_rot_from_quat (TRANSFORM r, QUATERNION q);
     void find_quat_from_rot (QUATERNION q, TRANSFORM r);
     void mult_quat (QUATERNION r, QUATERNION a, QUATERNION b);

// clear any storage for vertices in a polygon; after this,
// its state is the same as after init_polygon_2d()
     void clear_polygon_2d (POLYGON_2D * poly);

// compute minkowski difference B - A with distinguished point p
     void make_cso_polygon_2d (POLYGON_2D * r, POLYGON_2D * a, POINT_2D p,
			       POLYGON_2D * b);

// checks to see if p is left of or on all the edges of polygon a.
     int point_inside_convex_polygon_2d_p (POINT_2D p, POLYGON_2D * a);

// checks to see if p is no more than eps right of all the edges of polygon a.
     int point_near_convex_polygon_2d_p (POINT_2D p, POLYGON_2D * a,
					 FLOAT eps);

#endif