/* 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 // 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