path: root/graphics/sketch/geometry.h

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