/* ---------------------------------------------------------------------

Copyright 1994    Christopher A. Huson

Octree representation for n-cube for doing inertial cut algorithm.
We have two types of nodes in the octree :

    non-terminal - an array of 8 pointers to the next level of
		   non-terminal or terminal.

    terminal     - . array of 8 counters for population counts
		   . array of 8 indicating what subpart of the
		     cube the point is in.
		   . three fields indicating the common bits of
		     the points this terminal represents.

The index at level I is the I'th bits of the (R,G,B) values combined.
In this case, the combination is

	      4 * R[i] + 2 * G[i] + B[i]

To give an example, suppose the depth of the tree is four (that is,
three non-terminals.)  We are going to increment the count for the
color 

    R(0100) G(0011) B(1110)

Level   R   G   B    Structure

                     +----+----+----+----+----+----+----+----+
  0     0   0   1    |    | 1  |    |    |    |    |    |    |
                     +----+----+----+----+----+----+----+----+
                            |
                            +------------+
                                         V
                     +----+----+----+----+----+----+----+----+
  1     1   0   1    |    |    |    |    |    | 5  |    |    |
                     +----+----+----+----+----+----+----+----+
                                                |
                                         +------+
                                         V
                     +----+----+----+----+----+----+----+----+
  2     0   1   1    |    |    |    | 3  |    |    |    |    |
                     +----+----+----+----+----+----+----+----+
                                      |
                                      +--+
                                         V
                     +----+----+----+----+----+----+----+----+
  3     0   1   0    |    |    | XX |    |    |    |    |    |
                     +----+----+----+----+----+----+----+----+

The array element marked "XX" in the last (or terminal) row, is
the element to be incremented.

This representation gives three advantages :

   *  It is not as sparse as a simple 3-D array implementation.  The
      minimum density at the terminal level is 1/8, and if there are
      N distinct colors with this minimum density, the number of
      non-terminal nodes to represent the data is N - 1, (if N is a
      power of 2), so the total structure size at this minimum density
      is

	  sizeof(nonterminal) * (N - 1) + N * sizeof(terminal)


      One hopes this doesn't grow too quickly.  With coherence, the
      number will be smaller.

   *  It is possible to gracefully degrade the accuracy of the
      representation.  If there are N bits per color, we can start
      constructing an N x N x N cube.  If space is not sufficient, we
      can summarize the data at that level and continue building an
      (N-1) x (N-1) x (N-1) cube, and so on.

   *  The cube can be traversed more efficiently for the purposes of
      doing the inertial cut.  We can walk just the terminal nodes,
      and all the data we need is stored there.

We will allocate space in chunk size which is a common multiple of the
sizes of the non-terminal and terminal nodes, so any space allocated
can be used for either.
--------------------------------------------------------------------- */
#define COLOR_S unsigned short     /* maximum of 65535 colors */
#define COLOR_REP_S unsigned short  /* object which holds the largest color representation */
/* In this case, we will be able to handle up to 48-bit color representations (if short is 16 bits) */
#define COLOR_REP_BITS 16
#define APPROXIMATE_BLOCK_SIZE 2048
#define BYTE char
#define UBYTE unsigned char

#ifndef TRUE
#define TRUE 1
#endif

#ifndef FALSE
#define FALSE 0
#endif

#ifdef MAIN
#define DEFINE 
#define INIT(x) = x
#else
#define DEFINE extern
#define INIT(x)
#endif

#ifdef DEBUG
DEFINE long debug_flag INIT(0);
#define DHEADER printf("<none>-")
#define DEBUGMSG(n,dlist) if(debug_flag>=n){DHEADER;printf dlist;};
#else
#define DEBUGMSG(n,dlist)
#endif

struct octree_nonterminal {
    struct octree_nonterminal *next_level[ 8 ];    /* Will be cast at terminal level */
    };

struct octree_terminal {
    unsigned long OT_term_cnt[ 8 ];
    COLOR_S OT_term_assign[ 8 ];
    COLOR_REP_S OT_red_val, OT_green_val, OT_blue_val;    /* because we cannot derive this at the terminal level */
    };

DEFINE struct octree_nonterminal *current_position[ COLOR_REP_BITS ];

DEFINE long blue_values[8];
DEFINE long green_values[8];
DEFINE long red_values[8];

struct octree_space {
    struct octree_space *next_bit;
    unsigned long OS_length;
    unsigned long OS_size;     /* size of the structures being parcelled out of this space */
    unsigned long OS_last_used;    /* index must be multiplied by OS_size to get byte offset. */
    UBYTE OS_space[ 1 ];   /* real size will be much larger. */
    };

DEFINE struct octree_space *free_space;
DEFINE unsigned long OS_block_size INIT( 0 );       /* will be set in initialize_octree */

DEFINE struct octree_space *octree_levels[ COLOR_REP_BITS ];
DEFINE unsigned long free_nodes INIT( 0 );          /* Number of nodes in the freespace */

DEFINE unsigned long gl_current_nbits INIT( 0 );    /* will be decremented as we contract the cube */
DEFINE unsigned long gl_bits_per_color INIT( 0 );   /* How many bits in each color are significant */
						    /* (from right).                               */
						    /* this should always be true :                */
						    /* gl_current_nbits <= gl_bits_per_color       */

DEFINE struct octree_nonterminal first_level;

DEFINE struct {
    struct octree_space *next_bit;
    unsigned long OS_length;
    unsigned long OS_size;     /* size of the structures being parcelled out of this space */
    unsigned long OS_last_used;    /* index must be multiplied by OS_size to get byte offset. */
    UBYTE OS_space[ sizeof(struct octree_nonterminal) * 8 ];
    } second_level;

struct inert_data {
    double xavg, yavg, zavg;
    double a, b, c, f, g, h;
    double xdot, ydot, zdot;
    double sigmasquare;
    double maxmoment;
    long consider_region;
    long plane_data_index;
    };

DEFINE struct inert_data *all_data INIT(NULL);

struct ary_data {
    long xd, yd, zd;
    long pixelcnt;
    struct ary_data *next_data;
    };

DEFINE struct ary_data *line_data;
DEFINE struct ary_data **color_levels;

/* ----------------------------------------------
Cut plane data handed back to ppmquant for deciding
which final color to map a particular pixel to.

Consists of the following fields :

      avgr, avgg, avgb -- in final color, the color
			  mapped to.  Always the center
			  of gravity of that stage of
			  split.
      cosr, cosg, cosb -- The direction cosines we use
			  to split the space.
      ifneg, ifnonneg  -- index of next planes if this
			  space was split.

If we are mapping color (r,g,b), then

    if (r - avgr) * cosr + 
       (g - avgg) * cosg +
       (b - avgb) * cosb < 0,
    then
	next plane to split by is indexed by ifneg
    else
	next plane to split by is indexed by ifnonneg
	endif
---------------------------------------------- */

struct plane_data {
    double avgr, avgg, avgb, cosr, cosg, cosb;
    long ifneg, ifnonneg;
    };

DEFINE struct plane_data *all_plane_data;
DEFINE long last_plane_used INIT(0);

DEFINE long last_inert_data INIT(0);
DEFINE long free_space_cnt INIT(0);
DEFINE long nregs;
DEFINE double red_weight INIT( 1.0 );
DEFINE double green_weight INIT( 1.0 );
DEFINE double blue_weight INIT( 1.0 );
DEFINE long do_max_moment INIT( 0 );
DEFINE long do_best_possible INIT( 0 );
DEFINE long picturesize INIT( 0 );

/* prototypes */
#ifdef __STDC__
extern void initialize_octree( long nbits_per_color, long nbits_per_side, long num_of_registers );
extern void increment_cell( long red_val, long green_val, long blue_val, long inc_value );
extern void clean_up_octree( );
extern void reduce_octree( );
extern colorhist_vector assign_registers( );
extern void sum_group(long group_no, struct inert_data *data_ptr );
extern void inertia(double a, double b, double c, double f, double g, double h,
    double *xaxis, double *yaxis, double *zaxis, double *maxmom );
extern void split_group( long indx );
extern void cubic_roots( double a1, double a2, double a3, double *r0, double *r1r, double *r1i, double *r2r, double *r2i );
extern void cube_root( double rr, double ri, double *r0r, double *r0i );
extern void merge_group( long toregion, long fromregion );
extern long count_distinct_colors( );
extern void quadratic_roots(double a, double b, double c, double *x1r, double *x1i, double *x2r, double *x2i, double error);
#else /* __STDC__ */
extern VOID initialize_octree( );
extern VOID increment_cell( );
extern VOID clean_up_octree( );
extern VOID reduce_octree( );
extern colorhist_vector assign_registers( );
extern VOID sum_group( );
extern VOID inertia( );
extern VOID split_group( );
extern VOID cubic_roots( );
extern VOID cube_root( );
extern VOID merge_group( );
extern long count_distinct_colors( );
extern VOID quadratic_roots( );
#endif  /* __STDC__ */