v47i121: inertial_cut - Inertial Cut Algorithm for Color Quantization, Part01/02

[ch...@kai.com]

Submitted-by: ch...@kai.com
Posting-number: Volume 47, Issue 121
Archive-name: inertial_cut/part01
Environment: PBM

This program implements a new method of color quantization. Imagine
the data in the color cube as a rigid set of particles. The algorithm
finds an axis which has a minimal moment of inertia (that is, if the
data points were a solid which could be spun, it finds the axis for
which that solid has the smallest resistance to spinning.) The solid
is then split by a plane through the center of gravity and
perpendicular to this axis.

The algorithm is nearly always better than median cut (of the hundred
or so pictures I have processed, I have found one which median cut does
a better job on for three colors.) It is probably provably-optimal
for splitting any space in two. It is better-behaved as the number of
colors is increased.

The program is implemented as a part of the ppmquant program in the PBM
Toolkit, and requires that package.
--------
#! /bin/sh
# This is a shell archive. Remove anything before this line, then feed it
# into a shell via "sh file" or similar. To overwrite existing files,
# type "sh file -c".
# Contents: README Mfile.patch man.patch median.uu octree.h
# patchlevel.h quant.patch
# Wrapped by kent@ftp on Wed Apr 12 22:21:17 1995
PATH=/bin:/usr/bin:/usr/ucb:/usr/local/bin:/usr/lbin:$PATH ; export PATH
echo If this archive is complete, you will see the following message:
echo ' "shar: End of archive 1 (of 2)."'
if test -f 'README' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'README'\"
else
echo shar: Extracting \"'README'\" \(4203 characters\)
sed "s/^X//" >'README' <<'END_OF_FILE'
XINERTIAL CUT ALGORITHM FOR COLOR QUANTIZATION
Xby Christopher A. Huson <ch...@kai.com>
X
XIn the last decade or more median cut has been the best method
Xavailable for color quantization. Heckbert hinted other methods were
Xpossible. This is one such method.
X
XThe motivations for developing the algorithm were two:
X
X . I had the idea while studying statics as an undergraduate.
X
X . I have an Amiga computer which could only display 16 or
X 32 distinct colors. This was depressing.
X
X- - - - - - - - - - - -
X
XThis is an addition to the ppmquant command in PBM. The original
Xfunctionality of ppmquant is retained.
X
XThe version of PBM this is a patch to is
X PBMPLUS_VERSION "10dec91"
X
XThe register allocation scheme is called "inertial cut", from its
Xderivation (from a problem in statics.) At any stage this procedure
Xshould be provably-optimal, given we only allow splitting a region in
Xhalf, each space is considered independently, and given the traditional
Xmeasure of error
X
X 2 2 2
X error = sigma( dx + dy + dz )
X
XI have a paper describing the algorithm, available upon request. (I
Xdid get tired trying to prove the optimality; maybe someone with more
Xrecent math experience can succeed where I failed.)
X
XThis shar file contains
X
X README - this file
X
X quant.patch - The patch to the file ppmquant.c in the ppm
X subdirectory to incorporate the new function.
X
X man.patch - Patch to ppmquant.1 to add the new switches and
X descriptions.
X
X Mfile.patch - Patch to ppm/Makefile to build the new version
X of ppmquant.
X
X octree.h - Declarations for the structures.
X
X inert_cut.c - The algorithm.
X
X median.uu - Example file where (for three colors) median cut
X does better than inertial cut. (This is quite
X possible because each space is treated as an
X independent problem, but the final pixel mapping
X depends on the nearness of _all_ colors.)
X
X patchlevel.h - Version number for inertial cut algorithm.
X
X
XIn addition to the inertial cut algorithm, this source contains an
Xalternate representation of the "color cube"; it is a "gracefully-degrading
Xoctree". At interior nodes of the octree each corresponding bit of
Xthe (rgb) values is(are?) used to construct an index to an array which
Xpoints to the next nonterminal node, or to the population counts at the
Xterminal level. As long as space can be malloced the octree will
Xretain the full precision of the original data. When this is not
Xpossible, the terminal level of the octree is "harvested" and the least
Xsignificant bit of the octree is discarded. On our system I have never
Xhad this happen, and I have processed files with as many as a
Xquarter-million distinct colors. Of course our systems have 500-meg
Xswap partitions.
X
XThe program has been successfully installed on SGI Indigos. It has
Xbeen compiled on other systems (IBM, HP, Sun and Sequent) but has not
Xbeen installed (I did not have PBM libraries on the other systems.)
X
XThe only non-ANSI-C code I am aware of in my new stuff is that I use
Xpointers to store integers for the "harvesting" operation. This may
Xcause problems on machines with 65K segment memory addressing.
X
XThe possibilities for future development are:
X
X . Someone could make a HAM-specific colormap generator by
X weighting each pixel by the apparent nearness of the leftward
X pixel.
X
X . If we store each of the cutplanes as we generate them, each
X pixel can be assigned a final value based on which subspace it
X appears in, instead of checking for nearness of the pixel to
X each of the final colors. This will be slightly less "good"
X than checking each color, but for large numbers of colors would
X be phenomenally faster. (It would make assigning 32,000 colors
X possible, for example.)
X
XThose interested in the algorithm should take a look at the following
Xtwo references:
X
X Becker, Robert A., _Introduction to Theoretical Mechanics_,
X McGraw-Hill Book Co., Inc, 1954, pp 272-288.
X
X Synge, John L. & Griffith, Byron A., _Principles of Mechanics_,
X Second Edition, McGraw-Hill Book Co., 1949, pp 313-320.
X
END_OF_FILE
if test 4203 -ne `wc -c <'README'`; then
echo shar: \"'README'\" unpacked with wrong size!
fi
# end of 'README'
fi
if test -f 'Mfile.patch' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'Mfile.patch'\"
else
echo shar: Extracting \"'Mfile.patch'\" \(2710 characters\)
sed "s/^X//" >'Mfile.patch' <<'END_OF_FILE'
X*** Makefile.orig Tue Mar 7 00:21:14 1995
X--- Makefile Tue Mar 7 01:47:08 1995
X***************
X*** 45,54 ****
X ALLCFLAGS = $(CFLAGS) $(RGBDEF) $(INCLUDE)
X LIBPPM = libppm.a
X
X PORTBINARIES = giftoppm gouldtoppm ilbmtoppm imgtoppm mtvtoppm \
X pcxtoppm pgmtoppm pi1toppm picttoppm \
X pjtoppm ppmdither ppmhist ppmmake \
X! ppmquant ppmrelief ppmtoacad ppmtogif ppmtoicr \
X ppmtoilbm ppmtopcx ppmtopgm ppmtopi1 ppmtopict \
X ppmtopj ppmtopuzz ppmtorgb3 ppmtosixel \
X ppmtotga ppmtouil ppmtoxpm ppmtoyuv qrttoppm \
X--- 45,55 ----
X ALLCFLAGS = $(CFLAGS) $(RGBDEF) $(INCLUDE)
X LIBPPM = libppm.a
X
X+ QUANTBINARY = ppmquant
X PORTBINARIES = giftoppm gouldtoppm ilbmtoppm imgtoppm mtvtoppm \
X pcxtoppm pgmtoppm pi1toppm picttoppm \
X pjtoppm ppmdither ppmhist ppmmake \
X! ppmrelief ppmtoacad ppmtogif ppmtoicr \
X ppmtoilbm ppmtopcx ppmtopgm ppmtopi1 ppmtopict \
X ppmtopj ppmtopuzz ppmtorgb3 ppmtosixel \
X ppmtotga ppmtouil ppmtoxpm ppmtoyuv qrttoppm \
X***************
X*** 55,61 ****
X rawtoppm rgb3toppm sldtoppm spctoppm sputoppm \
X tgatoppm ximtoppm xpmtoppm yuvtoppm
X MATHBINARIES = ppmforge ppmpat
X! BINARIES = $(PORTBINARIES) $(MATHBINARIES)
X SCRIPTS = ppmquantall
X
X OBJECTS = giftoppm.o gouldtoppm.o ilbmtoppm.o imgtoppm.o mtvtoppm.o \
X--- 56,62 ----
X rawtoppm rgb3toppm sldtoppm spctoppm sputoppm \
X tgatoppm ximtoppm xpmtoppm yuvtoppm
X MATHBINARIES = ppmforge ppmpat
X! BINARIES = $(PORTBINARIES) $(MATHBINARIES) $(QUANTBINARY)
X SCRIPTS = ppmquantall
X
X OBJECTS = giftoppm.o gouldtoppm.o ilbmtoppm.o imgtoppm.o mtvtoppm.o \
X***************
X*** 67,73 ****
X ppmtotga.o ppmtouil.o ppmtoxpm.o ppmtoyuv.o qrttoppm.o \
X rawtoppm.o rgb3toppm.o sldtoppm.o spctoppm.o sputoppm.o \
X tgatoppm.o ximtoppm.o xpmtoppm.o yuvtoppm.o \
X! ppmforge.o ppmpat.o
X
X MANUALS1 = $(BINARIES) $(SCRIPTS)
X MANUALS3 = libppm
X--- 68,74 ----
X ppmtotga.o ppmtouil.o ppmtoxpm.o ppmtoyuv.o qrttoppm.o \
X rawtoppm.o rgb3toppm.o sldtoppm.o spctoppm.o sputoppm.o \
X tgatoppm.o ximtoppm.o xpmtoppm.o yuvtoppm.o \
X! ppmforge.o ppmpat.o inert_cut.o
X
X MANUALS1 = $(BINARIES) $(SCRIPTS)
X MANUALS3 = libppm
X***************
X*** 127,132 ****
X--- 128,137 ----
X # Rule for math-dependent programs.
X $(MATHBINARIES): ppm.h $(DEFPGM) $(DEFPBM) $(LIBPPM) $(LIBPGM) $(LIBPBM)
X $(CC) $(ALLCFLAGS) $(LDFLAGS) -o $@ $@.c -lm $(LIBPPM) $(LIBPGM) $(LIBPBM)
X+
X+ #Special rule for ppmquant (needs both inert_cut and math libraries)
X+ $(QUANTBINARY): ppm.h inert_cut.o $(DEFPGM) $(DEFPBM) $(LIBPPM) $(LIBPGM) $(LIBPBM)
X+ $(CC) $(ALLCFLAGS) $(LDFLAGS) -o $@ $@.c inert_cut.o -lm $(LIBPPM) $(LIBPGM) $(LIBPBM)
X
X # Rule for objects.
X $(OBJECTS): ppm.h $(DEFPGM) $(DEFPBM)
END_OF_FILE
if test 2710 -ne `wc -c <'Mfile.patch'`; then
echo shar: \"'Mfile.patch'\" unpacked with wrong size!
fi
# end of 'Mfile.patch'
fi
if test -f 'man.patch' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'man.patch'\"
else
echo shar: Extracting \"'man.patch'\" \(3061 characters\)
sed "s/^X//" >'man.patch' <<'END_OF_FILE'
X*** orig.ppmquant.1 Mon Mar 6 22:41:41 1995
X--- ppmquant.1 Tue Mar 7 01:38:18 1995
X***************
X*** 4,15 ****
X--- 4,21 ----
X ppmquant - quantize the colors in a portable pixmap down to a specified number
X .SH SYNOPSIS
X .B ppmquant
X+ .RB [ -inertial ]
X .RB [ -floyd | -fs ]
X+ .RB [ -[no]fsweight ]
X+ .RB [ -weight ]
X+ .RB [ -maxmoment ]
X+ .RB [ -best ]
X .I ncolors
X .RI [ ppmfile ]
X .br
X .B ppmquant
X .RB [ -floyd | -fs ]
X+ .RB [ -[no]fsweight ]
X .B -map
X .I mapfile
X .RI [ ppmfile ]
X***************
X*** 22,30 ****
X to the new ones, and writes a portable pixmap as output.
X .IX "colormap reduction"
X .PP
X! The quantization method is Heckbert's "median cut".
X .IX "median cut"
X .PP
X Alternately, you can skip the color-choosing step by
X specifying your own set of colors with the
X .B -map
X--- 28,40 ----
X to the new ones, and writes a portable pixmap as output.
X .IX "colormap reduction"
X .PP
X! The default quantization method is Heckbert's "median cut".
X .IX "median cut"
X .PP
X+ If the
X+ .B -inertial
X+ switch is specified, the "inertial cut" algorithm will be used to pick colors.
X+ .PP
X Alternately, you can skip the color-choosing step by
X specifying your own set of colors with the
X .B -map
X***************
X*** 62,68 ****
X--- 72,101 ----
X quantization has banding or other artifacts, especially when going to a
X small number of colors such as the above IBM set.
X However, it does take substantially more CPU time, so the default is off.
X+ The
X+ .B -fsweight
X+ switch will use the apparent brightness of each primary in weighting the
X+ error function.
X .PP
X+ When the "inertial cut" algorithm is specified, the following switches may
X+ be used to modify its operation:
X+ .TP 11
X+ .B -weight
X+ will use apparent brightness of each primary in sizing the dataspace to
X+ determine the cut plane. This usually results in a distinct green tinge to the
X+ picture, but is useful when large amounts of fleshtones or blues appear in the
X+ original.
X+ .TP 11
X+ .B -maxmoment
X+ picks the space to divide which has the largest maximum moment of inertia.
X+ By default the space with the largest sigma squared error is chosen.
X+ .TP 11
X+ .B -best
X+ will, when picking the last few colors, choose between the space with the
X+ largest sigma-squared error and the largest maximal moment of inertia
X+ (if they differ) based on which cut will result in greater reduction of
X+ the sigma-squared error of the mapping.
X+ .PP
X All flags can be abbreviated to their shortest unique prefix.
X .SH REFERENCES
X "Color Image Quantization for Frame Buffer Display" by Paul Heckbert,
X***************
X*** 71,76 ****
X--- 104,111 ----
X ppmquantall(1), pnmdepth(1), ppmdither(1), ppm(5)
X .SH AUTHOR
X Copyright (C) 1989, 1991 by Jef Poskanzer.
X+ .PP
X+ Inertial Cut Algorithm Copyright (C) 1994 by Christopher A. Huson
X .\" Permission to use, copy, modify, and distribute this software and its
X .\" documentation for any purpose and without fee is hereby granted, provided
X .\" that the above copyright notice appear in all copies and that both that
END_OF_FILE
if test 3061 -ne `wc -c <'man.patch'`; then
echo shar: \"'man.patch'\" unpacked with wrong size!
fi
# end of 'man.patch'
fi
if test -f 'median.uu' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'median.uu'\"
else
echo shar: Extracting \"'median.uu'\" \(650 characters\)
sed "s/^X//" >'median.uu' <<'END_OF_FILE'
Xbegin 400 median_better.ppm
XM4#8*,3(@,3(*,C4U"@$``0$``08%!3PT*WIF39I]69!S4&Y8/"PB&`$``0$`
XM`0$``0$``1<3$7]O6K*3;;^:;<":;+".8[Z7:+6/86Q6/0(!`@$``0D'!XU]
XM9Y2`9*.(9J:'8K219:6&7J^-8;^8:=2H<V!..`$``4M&/9N)<Y)^9(=Q5II]
XM6II]5YI\6)Q_6J^-8\*<;+^;;@@'"'!I8(-V97AL77YM5J""6YZ`6IY_6:6%
XM7*N)7[208\FA<2<E(W%L96)@6E=74F-=49=\6HYW685P59E]6)]_6*:#6)Q_
XM6C$P+E=965565$I-3$]/3%!.2%=325!-1D9%06A;2))U4(YS4"TK*#(V.6->
XM56A?4&M?36Y@3G-C35E21U-,05Y51VI=2W9E3!D7%@P.$D=)2&E=3'QI3X1M
XM3X-K3'ID1WAC1WAF2WAH4&-71P(!`@$``14;($I'06YA3G5D3'9C2GEF2WQI
XM3W)C3GMR8AT;&`$``0$``0$``0\3%S<X-EI2171H57]Q7H9Y9WQR8BHG(P$`
XH`0$``0$``0$``0$``0("!!`3%B(D)R<G)Q44%`(!`@$``0$``0$``2HG
X`
Xend
END_OF_FILE
if test 650 -ne `wc -c <'median.uu'`; then
echo shar: \"'median.uu'\" unpacked with wrong size!
fi
# end of 'median.uu'
fi
if test -f 'octree.h' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'octree.h'\"
else
echo shar: Extracting \"'octree.h'\" \(9232 characters\)
sed "s/^X//" >'octree.h' <<'END_OF_FILE'
X/* ---------------------------------------------------------------------
X
XCopyright 1994 Christopher A. Huson
X
XOctree representation for n-cube for doing inertial cut algorithm.
XWe have two types of nodes in the octree :
X
X non-terminal - an array of 8 pointers to the next level of
X non-terminal or terminal.
X
X terminal - . array of 8 counters for population counts
X . array of 8 indicating what subpart of the
X cube the point is in.
X . three fields indicating the common bits of
X the points this terminal represents.
X
XThe index at level I is the I'th bits of the (R,G,B) values combined.
XIn this case, the combination is
X
X 4 * R[i] + 2 * G[i] + B[i]
X
XTo give an example, suppose the depth of the tree is four (that is,
Xthree non-terminals.) We are going to increment the count for the
Xcolor
X
X R(0100) G(0011) B(1110)
X
XLevel R G B Structure
X
X +----+----+----+----+----+----+----+----+
X 0 0 0 1 | | 1 | | | | | | |
X +----+----+----+----+----+----+----+----+
X |
X +------------+
X V
X +----+----+----+----+----+----+----+----+
X 1 1 0 1 | | | | | | 5 | | |
X +----+----+----+----+----+----+----+----+
X |
X +------+
X V
X +----+----+----+----+----+----+----+----+
X 2 0 1 1 | | | | 3 | | | | |
X +----+----+----+----+----+----+----+----+
X |
X +--+
X V
X +----+----+----+----+----+----+----+----+
X 3 0 1 0 | | | XX | | | | | |
X +----+----+----+----+----+----+----+----+
X
XThe array element marked "XX" in the last (or terminal) row, is
Xthe element to be incremented.
X
XThis representation gives three advantages :
X
X * It is not as sparse as a simple 3-D array implementation. The
X minimum density at the terminal level is 1/8, and if there are
X N distinct colors with this minimum density, the number of
X non-terminal nodes to represent the data is N - 1, (if N is a
X power of 2), so the total structure size at this minimum density
X is
X
X sizeof(nonterminal) * (N - 1) + N * sizeof(terminal)
X
X
X One hopes this doesn't grow too quickly. With coherence, the
X number will be smaller.
X
X * It is possible to gracefully degrade the accuracy of the
X representation. If there are N bits per color, we can start
X constructing an N x N x N cube. If space is not sufficient, we
X can summarize the data at that level and continue building an
X (N-1) x (N-1) x (N-1) cube, and so on.
X
X * The cube can be traversed more efficiently for the purposes of
X doing the inertial cut. We can walk just the terminal nodes,
X and all the data we need is stored there.
X
XWe will allocate space in chunk size which is a common multiple of the
Xsizes of the non-terminal and terminal nodes, so any space allocated
Xcan be used for either.
X--------------------------------------------------------------------- */
X#define COLOR_S unsigned short /* maximum of 65535 colors */
X#define COLOR_REP_S unsigned short /* object which holds the largest color representation */
X/* In this case, we will be able to handle up to 48-bit color representations (if short is 16 bits) */
X#define COLOR_REP_BITS 16
X#define APPROXIMATE_BLOCK_SIZE 2048
X#define BYTE char
X#define UBYTE unsigned char
X
X#ifndef TRUE
X#define TRUE 1
X#endif
X
X#ifndef FALSE
X#define FALSE 0
X#endif
X
X#ifdef MAIN
X#define DEFINE
X#define INIT(x) = x
X#else
X#define DEFINE extern
X#define INIT(x)
X#endif
X
X#ifdef DEBUG
XDEFINE long debug_flag INIT(0);
X#define DHEADER printf("<none>-")
X#define DEBUGMSG(n,dlist) if(debug_flag>=n){DHEADER;printf dlist;};
X#else
X#define DEBUGMSG(n,dlist)
X#endif
X
Xstruct octree_nonterminal {
X struct octree_nonterminal *next_level[ 8 ]; /* Will be cast at terminal level */
X };
X
Xstruct octree_terminal {
X unsigned long OT_term_cnt[ 8 ];
X COLOR_S OT_term_assign[ 8 ];
X COLOR_REP_S OT_red_val, OT_green_val, OT_blue_val; /* because we cannot derive this at the terminal level */
X };
X
XDEFINE struct octree_nonterminal *current_position[ COLOR_REP_BITS ];
X
XDEFINE long blue_values[8];
XDEFINE long green_values[8];
XDEFINE long red_values[8];
X
Xstruct octree_space {
X struct octree_space *next_bit;
X unsigned long OS_length;
X unsigned long OS_size; /* size of the structures being parcelled out of this space */
X unsigned long OS_last_used; /* index must be multiplied by OS_size to get byte offset. */
X UBYTE OS_space[ 1 ]; /* real size will be much larger. */
X };
X
XDEFINE struct octree_space *free_space;
XDEFINE unsigned long OS_block_size INIT( 0 ); /* will be set in initialize_octree */
X
XDEFINE struct octree_space *octree_levels[ COLOR_REP_BITS ];
XDEFINE unsigned long free_nodes INIT( 0 ); /* Number of nodes in the freespace */
X
XDEFINE unsigned long gl_current_nbits INIT( 0 ); /* will be decremented as we contract the cube */
XDEFINE unsigned long gl_bits_per_color INIT( 0 ); /* How many bits in each color are significant */
X /* (from right). */
X /* this should always be true : */
X /* gl_current_nbits <= gl_bits_per_color */
X
XDEFINE struct octree_nonterminal first_level;
X
XDEFINE struct {
X struct octree_space *next_bit;
X unsigned long OS_length;
X unsigned long OS_size; /* size of the structures being parcelled out of this space */
X unsigned long OS_last_used; /* index must be multiplied by OS_size to get byte offset. */
X UBYTE OS_space[ sizeof(struct octree_nonterminal) * 8 ];
X } second_level;
X
Xstruct inert_data {
X double xavg, yavg, zavg;
X double a, b, c, f, g, h;
X double xdot, ydot, zdot;
X double sigmasquare;
X double maxmoment;
X long consider_region;
X long plane_data_index;
X };
X
XDEFINE struct inert_data *all_data INIT(NULL);
X
Xstruct ary_data {
X long xd, yd, zd;
X long pixelcnt;
X struct ary_data *next_data;
X };
X
XDEFINE struct ary_data *line_data;
XDEFINE struct ary_data **color_levels;
X
X/* ----------------------------------------------
XCut plane data handed back to ppmquant for deciding
Xwhich final color to map a particular pixel to.
X
XConsists of the following fields :
X
X avgr, avgg, avgb -- in final color, the color
X mapped to. Always the center
X of gravity of that stage of
X split.
X cosr, cosg, cosb -- The direction cosines we use
X to split the space.
X ifneg, ifnonneg -- index of next planes if this
X space was split.
X
XIf we are mapping color (r,g,b), then
X
X if (r - avgr) * cosr +
X (g - avgg) * cosg +
X (b - avgb) * cosb < 0,
X then
X next plane to split by is indexed by ifneg
X else
X next plane to split by is indexed by ifnonneg
X endif
X---------------------------------------------- */
X
Xstruct plane_data {
X double avgr, avgg, avgb, cosr, cosg, cosb;
X long ifneg, ifnonneg;
X };
X
XDEFINE struct plane_data *all_plane_data;
XDEFINE long last_plane_used INIT(0);
X
XDEFINE long last_inert_data INIT(0);
XDEFINE long free_space_cnt INIT(0);
XDEFINE long nregs;
XDEFINE double red_weight INIT( 1.0 );
XDEFINE double green_weight INIT( 1.0 );
XDEFINE double blue_weight INIT( 1.0 );
XDEFINE long do_max_moment INIT( 0 );
XDEFINE long do_best_possible INIT( 0 );
XDEFINE long picturesize INIT( 0 );
X
X/* prototypes */
X#ifdef __STDC__
Xextern void initialize_octree( long nbits_per_color, long nbits_per_side, long num_of_registers );
Xextern void increment_cell( long red_val, long green_val, long blue_val, long inc_value );
Xextern void clean_up_octree( );
Xextern void reduce_octree( );
Xextern colorhist_vector assign_registers( );
Xextern void sum_group(long group_no, struct inert_data *data_ptr );
Xextern void inertia(double a, double b, double c, double f, double g, double h,
X double *xaxis, double *yaxis, double *zaxis, double *maxmom );
Xextern void split_group( long indx );
Xextern void cubic_roots( double a1, double a2, double a3, double *r0, double *r1r, double *r1i, double *r2r, double *r2i );
Xextern void cube_root( double rr, double ri, double *r0r, double *r0i );
Xextern void merge_group( long toregion, long fromregion );
Xextern long count_distinct_colors( );
Xextern void quadratic_roots(double a, double b, double c, double *x1r, double *x1i, double *x2r, double *x2i, double error);
X#else /* __STDC__ */
Xextern VOID initialize_octree( );
Xextern VOID increment_cell( );
Xextern VOID clean_up_octree( );
Xextern VOID reduce_octree( );
Xextern colorhist_vector assign_registers( );
Xextern VOID sum_group( );
Xextern VOID inertia( );
Xextern VOID split_group( );
Xextern VOID cubic_roots( );
Xextern VOID cube_root( );
Xextern VOID merge_group( );
Xextern long count_distinct_colors( );
Xextern VOID quadratic_roots( );
X#endif /* __STDC__ */
END_OF_FILE
if test 9232 -ne `wc -c <'octree.h'`; then
echo shar: \"'octree.h'\" unpacked with wrong size!
fi
# end of 'octree.h'
fi
if test -f 'patchlevel.h' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'patchlevel.h'\"
else
echo shar: Extracting \"'patchlevel.h'\" \(76 characters\)
sed "s/^X//" >'patchlevel.h' <<'END_OF_FILE'
X/* Revision of Inertial Cut software */
X
X#define INERTIAL_VERSION 1995mar01
END_OF_FILE
if test 76 -ne `wc -c <'patchlevel.h'`; then
echo shar: \"'patchlevel.h'\" unpacked with wrong size!
fi
# end of 'patchlevel.h'
fi
if test -f 'quant.patch' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'quant.patch'\"
else
echo shar: Extracting \"'quant.patch'\" \(9167 characters\)
sed "s/^X//" >'quant.patch' <<'END_OF_FILE'
X*** orig.ppmquant.c Mon Mar 6 22:41:23 1995
X--- ppmquant.c Mon Mar 6 22:36:36 1995
X***************
X*** 10,19 ****
X--- 10,24 ----
X ** implied warranty.
X */
X
X+ #define MAIN
X #include "ppm.h"
X #include "ppmcmap.h"
X+ #include "octree.h"
X
X #define MAXCOLORS 32767
X+ #define FS_RED_WEIGHT 0.299
X+ #define FS_GREEN_WEIGHT 0.587
X+ #define FS_BLUE_WEIGHT 0.114
X
X /* #define LARGE_NORM */
X #define LARGE_LUM
X***************
X*** 52,58 ****
X register int ind;
X colorhist_vector chv, colormap;
X colorhash_table cht;
X! int floyd;
X int usehash;
X long* thisrerr;
X long* nextrerr;
X--- 57,64 ----
X register int ind;
X colorhist_vector chv, colormap;
X colorhash_table cht;
X! int floyd, do_weight, do_inertial_cut, do_fs_weight;
X! double fs_green_weight, fs_red_weight, fs_blue_weight;
X int usehash;
X long* thisrerr;
X long* nextrerr;
X***************
X*** 61,75 ****
X long* thisberr;
X long* nextberr;
X long* temperr;
X register long sr, sg, sb, err;
X #define FS_SCALE 1024
X int fs_direction;
X! char* usage = "[-floyd|-fs] <ncolors> [ppmfile]\n [-floyd|-fs] -map mapfile [ppmfile]";
X
X ppm_init( &argc, argv );
X
X argn = 1;
X floyd = 0;
X mappixels = (pixel**) 0;
X
X while ( argn < argc && argv[argn][0] == '-' && argv[argn][1] != '\0' )
X--- 67,93 ----
X long* thisberr;
X long* nextberr;
X long* temperr;
X+ long ii,jj,rval,gval,bval;
X register long sr, sg, sb, err;
X #define FS_SCALE 1024
X int fs_direction;
X! char* usage = "[-floyd|-fs] [-[no]fsweight] <ncolors> [ppmfile]\n\
X! -inertial [-floyd|-fs] [-weight] [-[no]fsweight] [-maxmoment] [-best] <ncolors> [ppmfile]\n\
X! [-floyd|-fs] [-[no]fsweight] -map mapfile [ppmfile]\n\
X! -fsweight specifies weighted distances for Floyd-Steinberg dithering.\n\
X! For inertial cut algorithm :\n\
X! -maxmoment chooses by maximal maximum moment of inertia\n\
X! -best chooses the best by either moment or sigma near the end.\n\
X! -weight modifies values by apparent brightness.";
X
X ppm_init( &argc, argv );
X
X argn = 1;
X floyd = 0;
X+ do_inertial_cut = 0;
X+ do_weight = do_fs_weight = 0;
X+ do_max_moment = 0;
X+ do_best_possible = 0;
X mappixels = (pixel**) 0;
X
X while ( argn < argc && argv[argn][0] == '-' && argv[argn][1] != '\0' )
X***************
X*** 80,85 ****
X--- 98,115 ----
X else if ( pm_keymatch( argv[argn], "-nofs", 2 ) ||
X pm_keymatch( argv[argn], "-nofloyd", 2 ) )
X floyd = 0;
X+ else if ( pm_keymatch( argv[argn], "-inertial", 2 ) )
X+ do_inertial_cut = 1;
X+ else if ( pm_keymatch( argv[argn], "-weight", 2 ) )
X+ do_weight = do_fs_weight = 1;
X+ else if ( pm_keymatch( argv[argn], "-maxmoment", 3 ) )
X+ do_max_moment = 1;
X+ else if ( pm_keymatch( argv[argn], "-best", 2 ) )
X+ do_best_possible = 1;
X+ else if ( pm_keymatch( argv[argn], "-fsweight", 3 ) )
X+ do_fs_weight = 1;
X+ else if ( pm_keymatch( argv[argn], "-nofsweight", 5 ) )
X+ do_fs_weight = 0;
X else if ( pm_keymatch( argv[argn], "-map", 2 ) )
X {
X ++argn;
X***************
X*** 127,163 ****
X
X if ( mappixels == (pixel**) 0 )
X {
X! /*
X! ** Step 2: attempt to make a histogram of the colors, unclustered.
X! ** If at first we don't succeed, lower maxval to increase color
X! ** coherence and try again. This will eventually terminate, with
X! ** maxval at worst 15, since 32^3 is approximately MAXCOLORS.
X! */
X! for ( ; ; )
X {
X! pm_message( "making histogram..." );
X! chv = ppm_computecolorhist(
X! pixels, cols, rows, MAXCOLORS, &colors );
X! if ( chv != (colorhist_vector) 0 )
X! break;
X! pm_message( "too many colors!" );
X! newmaxval = maxval / 2;
X! pm_message(
X! "scaling colors from maxval=%d to maxval=%d to improve clustering...",
X! maxval, newmaxval );
X! for ( row = 0; row < rows; ++row )
X! for ( col = 0, pP = pixels[row]; col < cols; ++col, ++pP )
X! PPM_DEPTH( *pP, *pP, maxval, newmaxval );
X! maxval = newmaxval;
X! }
X! pm_message( "%d colors found", colors );
X
X! /*
X! ** Step 3: apply median-cut to histogram, making the new colormap.
X! */
X! pm_message( "choosing %d colors...", newcolors );
X! colormap = mediancut( chv, colors, rows * cols, maxval, newcolors );
X! ppm_freecolorhist( chv );
X }
X else
X {
X--- 157,230 ----
X
X if ( mappixels == (pixel**) 0 )
X {
X! if(!do_inertial_cut)
X {
X! /*
X! ** Step 2: attempt to make a histogram of the colors, unclustered.
X! ** If at first we don't succeed, lower maxval to increase color
X! ** coherence and try again. This will eventually terminate, with
X! ** maxval at worst 15, since 32^3 is approximately MAXCOLORS.
X! */
X! for ( ; ; )
X! {
X! pm_message( "making histogram..." );
X! chv = ppm_computecolorhist(
X! pixels, cols, rows, MAXCOLORS, &colors );
X! if ( chv != (colorhist_vector) 0 )
X! break;
X! pm_message( "too many colors!" );
X! newmaxval = maxval / 2;
X! pm_message(
X! "scaling colors from maxval=%d to maxval=%d to improve clustering...",
X! maxval, newmaxval );
X! for ( row = 0; row < rows; ++row )
X! for ( col = 0, pP = pixels[row]; col < cols; ++col, ++pP )
X! PPM_DEPTH( *pP, *pP, maxval, newmaxval );
X! maxval = newmaxval;
X! }
X! pm_message( "%d colors found", colors );
X!
X! /*
X! ** Step 3: apply median-cut to histogram, making the new colormap.
X! */
X! pm_message( "choosing %d colors...", newcolors );
X! colormap = mediancut( chv, colors, rows * cols, maxval, newcolors );
X! ppm_freecolorhist( chv );
X! }
X! else
X! {
X! /* Now, create octree */
X
X! if(do_weight) {
X! set_color_weights( FS_RED_WEIGHT * 3.0, FS_GREEN_WEIGHT * 3.0, FS_BLUE_WEIGHT * 3.0 );
X! }
X! else {
X! set_color_weights( 1.0, 1.0, 1.0 );
X! }
X! for(ii=maxval,jj=0;ii>0;ii>>=1,jj++);
X! initialize_octree( jj, jj, newcolors );
X!
X! for(ii=0;ii<rows;ii++) {
X! for(jj=0;jj<cols;jj++) {
X! rval = PPM_GETR(pixels[ii][jj]);
X! gval = PPM_GETG(pixels[ii][jj]);
X! bval = PPM_GETB(pixels[ii][jj]);
X! increment_cell(rval,gval,bval,1);
X! }
X! }
X!
X! if(picturesize == 0) {
X! pm_error( "Source file is empty." );
X! }
X! assign_colors( );
X! colormap = (colorhist_vector)assign_registers( );
X! if(colormap){
X! }
X! else {
X! pm_error("out of memory" );
X! }
X! clean_up_octree( );
X! }
X }
X else
X {
X***************
X*** 185,190 ****
X--- 252,258 ----
X ** Step 4: map the colors in the image to their closest match in the
X ** new colormap, and write 'em out.
X */
X+
X pm_message( "mapping image to new colors..." );
X cht = ppm_alloccolorhash( );
X usehash = 1;
X***************
X*** 248,269 ****
X { /* No; search colormap for closest match. */
X register int i, r1, g1, b1, r2, g2, b2;
X register long dist, newdist;
X r1 = PPM_GETR( *pP );
X g1 = PPM_GETG( *pP );
X b1 = PPM_GETB( *pP );
X! dist = 2000000000;
X! for ( i = 0; i < newcolors; ++i )
X {
X! r2 = PPM_GETR( colormap[i].color );
X! g2 = PPM_GETG( colormap[i].color );
X! b2 = PPM_GETB( colormap[i].color );
X! newdist = ( r1 - r2 ) * ( r1 - r2 ) +
X! ( g1 - g2 ) * ( g1 - g2 ) +
X! ( b1 - b2 ) * ( b1 - b2 );
X! if ( newdist < dist )
X! {
X! ind = i;
X! dist = newdist;
X }
X }
X if ( usehash )
X--- 316,359 ----
X { /* No; search colormap for closest match. */
X register int i, r1, g1, b1, r2, g2, b2;
X register long dist, newdist;
X+ register double ddist, newddist;
X r1 = PPM_GETR( *pP );
X g1 = PPM_GETG( *pP );
X b1 = PPM_GETB( *pP );
X! if(!do_fs_weight )
X! {
X! dist = 2000000000;
X! for ( i = 0; i < newcolors; ++i )
X! {
X! r2 = PPM_GETR( colormap[i].color );
X! g2 = PPM_GETG( colormap[i].color );
X! b2 = PPM_GETB( colormap[i].color );
X! newdist = ( r1 - r2 ) * ( r1 - r2 ) +
X! ( g1 - g2 ) * ( g1 - g2 ) +
X! ( b1 - b2 ) * ( b1 - b2 );
X! if ( newdist < dist )
X! {
X! ind = i;
X! dist = newdist;
X! }
X! }
X! }
X! else
X {
X! ddist = 2000000000.0;
X! for ( i = 0; i < newcolors; ++i )
X! {
X! r2 = PPM_GETR( colormap[i].color );
X! g2 = PPM_GETG( colormap[i].color );
X! b2 = PPM_GETB( colormap[i].color );
X! newddist = (double)(( r1 - r2 ) * ( r1 - r2 ))*FS_RED_WEIGHT*FS_RED_WEIGHT +
X! (double)(( g1 - g2 ) * ( g1 - g2 ))*FS_GREEN_WEIGHT*FS_GREEN_WEIGHT +
X! (double)(( b1 - b2 ) * ( b1 - b2 ))*FS_BLUE_WEIGHT*FS_BLUE_WEIGHT;
X! if(newddist < ddist)
X! {
X! ind = i;
X! ddist = newddist;
X! }
X }
X }
X if ( usehash )
END_OF_FILE
if test 9167 -ne `wc -c <'quant.patch'`; then
echo shar: \"'quant.patch'\" unpacked with wrong size!
fi
# end of 'quant.patch'
fi
echo shar: End of archive 1 \(of 2\).
cp /dev/null ark1isdone
MISSING=""
for I in 1 2 ; do
if test ! -f ark${I}isdone ; then
MISSING="${MISSING} ${I}"
fi
done
if test "${MISSING}" = "" ; then
echo You have unpacked both archives.
rm -f ark[1-9]isdone
else
echo You still must unpack the following archives:
echo " " ${MISSING}
fi
exit 0
exit 0 # Just in case...

[ch...@kai.com]

Submitted-by: ch...@kai.com
Posting-number: Volume 47, Issue 122
Archive-name: inertial_cut/part02
Environment: PBM

#! /bin/sh
# This is a shell archive. Remove anything before this line, then feed it
# into a shell via "sh file" or similar. To overwrite existing files,
# type "sh file -c".

# Contents: inert_cut.c

# Wrapped by kent@ftp on Wed Apr 12 22:21:17 1995
PATH=/bin:/usr/bin:/usr/ucb:/usr/local/bin:/usr/lbin:$PATH ; export PATH
echo If this archive is complete, you will see the following message:

echo ' "shar: End of archive 2 (of 2)."'
if test -f 'inert_cut.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'inert_cut.c'\"
else
echo shar: Extracting \"'inert_cut.c'\" \(64461 characters\)
sed "s/^X//" >'inert_cut.c' <<'END_OF_FILE'
X/* --------------------------------------------------------------------\
X| |
X| Implementation of inertial-cut algorithm. The structure used to |
X| store the population counts is an octree addressed at level i by |
X| the i-th bits of the red, blue and green values. When space can |
X| no longer be allocated, the octree will gracefully degrade (if we |
X| start with 8 bits, we can reduce to a 7-bit octree and so on). |
X| |
X| After the structure is filled, a dense list-oriented representation |
X| is built. If we do not have space to create this structure, we can |
X| use the octree itself, though the process will be slower. |
X| |
X| Then, starting with the space as a whole, at each stage the largest |
X| space (either by maximum moment of inertia or by maximum sigma-- |
X| squared error) is chosen, and split in half by a plane through the |
X| "center of gravity" and perpendicular to the axis of minimum moment. |
X| If the resulting spaces will be mapped into the same output colors, |
X| the spaces are merged back to one space, and marked so we no longer |
X| consider splitting it. |
X| |
X| This code is meant to be called by ppmquant under the proper |
X| circumstances. |
X| |
X| Copyright 1994 Christopher A. Huson |
X| |
X| Permission to use, copy, modify, and distribute this software |
X| and its documentation for any purpose and without fee is hereby |
X| granted, provided that the above copyright notice appear in all |
X| copies and that both that copyright notice and this permission |
X| notice appear in supporting documentation. This software is |
X| provided "as is" without express or implied warranty. |
X| |
X| Thanks, Jef. |
X| |
X\-------------------------------------------------------------------- */
X
X/* ---------------------------------------------------------------------
X
XCALLING SEQUENCE :
X
X See ppmquant for an example of use of this package.
X
X . Set weights for red, green and blue by set_color_weights.
X
X set_color_weights( 1.0, 1.0, 1.0 );
X
X will result in an unweighted analysis;
X
X set_color_weights( 0.299, 0.587, 0.114 );
X
X will weight by luminosity;
X
X set_color_weights(0.897, 1.761, 0.342);
X
X (three times the luminosities) gives the same results, but makes
X the sums of the error weights comparable to the unweighted case.
X
X . Call initialize_octree( nbits, bits_per_side, num_of_registers),
X where nbits is the number of bits in each color
X
X ceil(log2(maximum-color-value))
X
X num_of_registers is the number of registers we are targeting,
X and bits_per_side is the number of bits initially used in making
X the octree. bits_per_side <= nbits, and usually will be equal.
X If bits_per_side is smaller, the space used by the octree is
X smaller. If bits_per_side is zero, the octree will be nbits per
X side. As we run out of space, the bits per side of the octree
X will be reduced dynamically.
X
X . For each pixel, call increment_cell( rval, gval, bval, incval ).
X If the use requires a non-one increment, it can be passed in
X incval, otherwise pass in one. (The non-one case I was thinking
X of was a treatment for HAM ILBMs (Amiga format.))
X
X . After this process completes, call assign_colors( ). This will
X set the colors for the terminal level of the octree.
X
X . Do
X
X colormap = (colorhist_vector)assign_registers( );
X
X The color map generated will be allocated on the heap.
X
X--------------------------------------------------------------------- */
X
X#include <stdio.h>
X#include <stdlib.h>
X
X#ifdef __STDC__
X#define VOID void
X#define VOIDPTR void *
X#else /* __STDC__ */
X#define VOID
X#define VOIDPTR char *
X#endif /* __STDC__ */
X
X#ifdef COPROCESSOR
X#include <m68881.h>
X#else /* COPROCESSOR */
X#include <math.h>
X#endif /* COPROCESSOR */
X
X#include "ppm.h"
X#include "ppmcmap.h"
X#include "octree.h"
X
X/* ---------------------------------------------------------------------
XThis procedure assigns the registers for the octree. It will
X
X 1) Set up the first space data structure.
X 2) While we have registers left to assign
X
X a) pick the space with the largest maximum moment of
X inertia.
X b) Walk through the terminal level of the tree structure,
X assigning a new group number to the color values
X which fall on the "wrong" side of the dividing
X plane,
X c) Re-calculate the data structures for the new spaces.
X
XWhen this process is done, for j in [0:nregs-1], the colors for the
Xnew register j will be all_data[j].{xavg,yavg,zavg}.
X--------------------------------------------------------------------- */
X
X#ifdef DEBUG
X#undef DHEADER
X#define DHEADER printf("\nassign_registers-");
X#endif
Xcolorhist_vector
Xassign_registers( )
X{
X colorhist_vector newmap;
X long new_x, new_y, new_z;
X long new_xval, new_yval, new_zval;
X long i, j, maxi, maxmomi, color_cnt;
X long index_to_split, report_mod;
X double maxval, maxmom, errorsum;
X double max_sigma, max_error;
X char tostring[100];
X /* First, adjust red_values, blue_values and green_values by the
X actual number of bits represented in the octree */
X
X i = gl_bits_per_color - gl_current_nbits;
X if ( i > 0) {
X for(j=0;j<8;j++) {
X red_values[j] <<= i;
X green_values[j] <<= i;
X blue_values[j] <<= i;
X }
X }
X
X /* now, set up the first of the data structures */
X
X all_data = (struct inert_data *)malloc(sizeof(struct inert_data) * nregs );
X if(!all_data)
X {
X pm_error( "Out of space" );
X }
X last_inert_data = 1;
X
X /* Free all the octree data except the terminal level. */
X
X clean_up_octree( TRUE );
X
X color_cnt = count_distinct_colors( );
X pm_message( "Total of %d distinct colors found", color_cnt );
X report_mod = nregs / 8;
X if(report_mod == 0) report_mod = 1;
X
X if( nregs > 8 ) {
X color_levels = (struct ary_data **)malloc( nregs * sizeof( struct ary_data **) );
X line_data = (struct ary_data *)malloc( color_cnt * sizeof( struct ary_data ) );
X if( line_data == NULL ) {
X pm_message( "Must use the octree structure" );
X }
X else {
X pm_message( "Restructuring data" );
X restructure_info( );
X clean_up_octree( FALSE );
X }
X }
X else {
X line_data = NULL;
X }
X
X all_plane_data = malloc( 2 * nregs * sizeof(struct plane_data));
X if (all_plane_data == NULL) {
X pm_error( "Out of space" );
X }
X last_plane_used = 0;
X
X sum_group( 0, &(all_data[0]));
X all_data[0].consider_region = TRUE;
X while( last_inert_data < nregs) {
X#ifdef DEBUG
X for(j=0;j<last_inert_data;j++) {
X DHEADER
X printf( "data[ %d ] : ( %g, %g, %g ), max=%g, consider=%s", j,
X all_data[j].xavg, all_data[j].yavg, all_data[j].zavg,
X all_data[j].maxmoment, (all_data[j].consider_region) ? "TRUE" : "FALSE" );
X }
X printf( "\n" );
X#endif
X maxi = -1; maxval = errorsum = 0;
X maxmom = 0; maxmomi = -1;
X for(j=0;j<last_inert_data;j++) {
X errorsum += all_data[j].sigmasquare;
X if(all_data[j].consider_region) {
X if(maxval < all_data[j].sigmasquare) {
X maxi = j;
X maxval = all_data[j].sigmasquare;
X }
X if(maxmom < all_data[j].maxmoment) {
X maxmomi = j;
X maxmom = all_data[j].maxmoment;
X }
X }
X }
X
X DEBUGMSG(10,( "maxi=%d, maxval=%g, maxmomi=%d, maxmom=%g\n", maxi, maxval, maxmomi, maxmom ));
X errorsum = errorsum / (double)picturesize;
X if(last_inert_data % report_mod == 0 ) {
X sprintf( tostring, "%d %s %s error per pixel %g",
X last_inert_data,
X (last_inert_data == 1) ? "color" : "colors",
X (last_inert_data == 1) ? "has" : "have",
X errorsum );
X pm_message( tostring );
X }
X if(maxval == 0.0) {
X pm_message( "There is no reason to split to more than %d colors", last_inert_data );
X for(j=last_inert_data;j<nregs;j++) {
X all_data[j].xavg = 0;
X all_data[j].yavg = 0;
X all_data[j].zavg = 0;
X }
X break;
X }
X else {
X
X /* We will only do the computationally-difficult part near the end of the
X conversion process. We suppose if one region has better chance to be
X split at any stage before the end, it will become the best choice by
X either method very soon. */
X
X if(do_best_possible && (maxmomi != maxi) && (last_inert_data > (nregs-7)) ) {
X max_error = all_data[maxmomi].sigmasquare;
X split_group( maxmomi );
X sum_group( maxmomi, &(all_data[maxmomi]));
X sum_group( last_inert_data, &(all_data[last_inert_data]));
X max_error = max_error - all_data[maxmomi].sigmasquare - all_data[last_inert_data].sigmasquare;
X merge_group( maxmomi, last_inert_data );
X sum_group( maxmomi, &(all_data[maxmomi]));
X max_sigma = all_data[maxi].sigmasquare;
X split_group( maxi );
X sum_group( maxi, &(all_data[maxi]));
X sum_group( last_inert_data, &(all_data[last_inert_data]));
X max_sigma = max_sigma - all_data[maxi].sigmasquare - all_data[last_inert_data].sigmasquare;
X merge_group( maxi, last_inert_data );
X sum_group( maxi, &(all_data[maxi]));
X if ( max_sigma < max_error ) {
X pm_message( "Split %d by maximum moment.", last_inert_data );
X index_to_split = maxmomi;
X }
X else {
X pm_message( "Split %d by maximum sigma-squared error.", last_inert_data );
X index_to_split = maxi;
X }
X }
X else {
X if(do_max_moment) {
X index_to_split = maxmomi;
X }
X else {
X index_to_split = maxi;
X }
X }
X split_group( index_to_split );
X sum_group( index_to_split, &(all_data[index_to_split]));
X sum_group( last_inert_data, &(all_data[last_inert_data]));
X new_xval = (all_data[last_inert_data].xavg/red_weight) + 0.5;
X new_yval = (all_data[last_inert_data].yavg/green_weight) + 0.5;
X new_zval = (all_data[last_inert_data].zavg/blue_weight) + 0.5;
X new_x = (all_data[index_to_split].xavg/red_weight) + 0.5;
X new_y = (all_data[index_to_split].yavg/green_weight) + 0.5;
X new_z = (all_data[index_to_split].zavg/blue_weight) + 0.5;
X if((new_xval==new_x)&&(new_yval==new_y)&&(new_zval==new_z)) {
X DEBUGMSG(10,( "Whoops! Shouldn't have split that group up.\n" ));
X merge_group( index_to_split, last_inert_data );
X sum_group( index_to_split, &(all_data[index_to_split]));
X all_data[index_to_split].consider_region = FALSE;
X pm_message( "Re-merging group; resolution of colormap exceeded." );
X }
X else {
X all_data[last_inert_data].consider_region = TRUE;
X last_inert_data += 1;
X }
X }
X }
X#ifdef DEBUG
X DEBUGMSG(10,( "After final assignment,\n"));
X for(j=0;j<last_inert_data;j++) {
X DHEADER
X printf( "data[ %d ] : ( %g, %g, %g ), max=%g", j,
X all_data[j].xavg, all_data[j].yavg, all_data[j].zavg,
X all_data[j].sigmasquare );
X }
X printf( "\n" );
X#endif
X errorsum = 0;
X for(j=0;j<last_inert_data;j++) {
X errorsum += all_data[j].sigmasquare;
X }
X errorsum = errorsum / (double) picturesize;
X pm_message( "Final error per pixel %g.", errorsum );
X
X clean_up_octree( FALSE );
X
X newmap = (colorhist_vector)malloc(sizeof(struct colorhist_item) * nregs);
X for(j=0;j<nregs;j++) {
X new_x = (all_data[j].xavg/red_weight) + 0.5;
X new_y = (all_data[j].yavg/green_weight) + 0.5;
X new_z = (all_data[j].zavg/blue_weight) + 0.5;
X PPM_ASSIGN(newmap[j].color, new_x, new_y, new_z );
X }
X
X free( all_data );
X
X return(newmap);
X } /* assign_registers */
X
X/* ---------------------------------------------------------------------
Xsplit group with identifier indx into two groups. The new identifier
X(the value in OT_term_assign[]) will be last_indx_used before it is
Xincremented.)
X--------------------------------------------------------------------- */
X#ifdef DEBUG
X#undef DHEADER
X#define DHEADER printf("\nsplit_group-");
X#endif
X#ifdef __STDC__
XVOID
Xsplit_group( long indx )
X#else /* __STDC__ */
XVOID
Xsplit_group( indx )
Xlong indx;
X#endif /* __STDC__ */
X{
X struct octree_space *tptr;
X struct octree_terminal *termptr;
X double rval, gval, bval, sum;
X long i, j;
X struct ary_data *lowptr, *highptr, *lowlast, *highlast, *ppptr, *nextppptr;
X
X if( line_data == NULL ) {
X for(tptr=octree_levels[gl_current_nbits-1];tptr;tptr = tptr->next_bit) {
X termptr = (struct octree_terminal *)&(tptr->OS_space[0]);
X for(i=0;i<=tptr->OS_last_used;i++) {
X for(j=0;j<8;j++) {
X if(indx == termptr[i].OT_term_assign[j]) {
X rval = (termptr[i].OT_red_val + red_values[j])*red_weight;
X gval = (termptr[i].OT_green_val + green_values[j])*green_weight;
X bval = (termptr[i].OT_blue_val + blue_values[j])*blue_weight;
X rval -= all_data[indx].xavg;
X gval -= all_data[indx].yavg;
X bval -= all_data[indx].zavg;
X sum = rval * all_data[indx].xdot +
X gval * all_data[indx].ydot +
X bval * all_data[indx].zdot;
X if(sum < 0.0) {
X termptr[i].OT_term_assign[j] = last_inert_data;
X }
X }
X }
X }
X }
X }
X else {
X lowptr = NULL;
X highptr = NULL;
X lowlast = NULL;
X highlast = NULL;
X for( ppptr = color_levels[indx];ppptr;ppptr=nextppptr) {
X nextppptr = ppptr->next_data;
X rval = ((ppptr->xd) * red_weight) - all_data[indx].xavg;
X gval = ((ppptr->yd) * green_weight) - all_data[indx].yavg;
X bval = ((ppptr->zd) * blue_weight) - all_data[indx].zavg;
X sum = rval * all_data[indx].xdot +
X gval * all_data[indx].ydot +
X bval * all_data[indx].zdot;
X if(sum < 0.0) {
X if(lowptr) {
X lowlast->next_data = ppptr;
X lowlast = ppptr;
X }
X else {
X lowlast = lowptr = ppptr;
X }
X lowlast->next_data = NULL;
X }
X else {
X if(highptr) {
X highlast->next_data = ppptr;
X highlast = ppptr;
X }
X else {
X highlast = highptr = ppptr;
X }
X highlast->next_data = NULL;
X }
X }
X color_levels[indx] = lowptr;
X color_levels[last_inert_data] = highptr;
X }
X } /* split_group */
X
X/* ---------------------------------------------------------------------
XIf we are working in the octree representation, rename the records
Xwith the value fromgroup to togroup. If we are working in the linked
Xlist structure, put the list fromgroup at the end of the list togroup.
X--------------------------------------------------------------------- */
X#ifdef DEBUG
X#undef DHEADER
X#define DHEADER printf("\nmerge_group-");
X#endif
X#ifdef __STDC__
XVOID
Xmerge_group( long togroup, long fromgroup )
X#else /* __STDC__ */
XVOID
Xmerge_group( togroup, fromgroup )
Xlong togroup;
Xlong fromgroup;
X#endif /* __STDC__ */
X{
X struct octree_space *tptr;
X struct octree_terminal *termptr;
X long i, j;
X
X struct ary_data *lowptr, *highptr, *resptr, *lastptr;
X
X if ( line_data == NULL ) {
X for(tptr=octree_levels[gl_current_nbits-1];tptr;tptr = tptr->next_bit) {
X termptr = (struct octree_terminal *)&(tptr->OS_space[0]);
X for(i=0;i<=tptr->OS_last_used;i++) {
X for(j=0;j<8;j++) {
X if(fromgroup == termptr[i].OT_term_assign[j]) {
X termptr[i].OT_term_assign[j] = togroup;
X }
X }
X }
X }
X }
X else {
X lowptr = color_levels[ togroup ];
X if(lowptr == NULL) {
X color_levels[ togroup ] = color_levels[ fromgroup ];
X return;
X }
X highptr = color_levels[ fromgroup ];
X if(highptr == NULL) {
X return;
X }
X while(lowptr->next_data != NULL) {
X lowptr = lowptr->next_data;
X }
X lowptr->next_data = highptr;
X color_levels[ fromgroup ] = NULL;
X }
X } /* merge_group */
X
X/* ---------------------------------------------------------------------
XFigure the average of the pixels in a particular group.
X--------------------------------------------------------------------- */
X
X#ifdef DEBUG
X#undef DHEADER
X#define DHEADER printf("\nsum_group-");
X#endif
X#ifdef __STDC__
XVOID
Xsum_group( long group_no, struct inert_data *data_ptr )
X#else /* __STDC__ */
XVOID
Xsum_group( group_no, data_ptr )
Xlong group_no;
Xstruct inert_data *data_ptr;
X#endif /* __STDC__ */
X{
X struct octree_terminal *termptr;
X struct octree_space *tptr;
X long i, popcnt, j, mass;
X double rval, gval, bval;
X double xmid, ymid, zmid;
X
X struct ary_data *ppptr;
X
X data_ptr->xavg = 0.0;
X data_ptr->yavg = 0.0;
X data_ptr->zavg = 0.0;
X popcnt = 0;
X
X if(line_data == NULL) {
X for(tptr=octree_levels[gl_current_nbits - 1];tptr;tptr=tptr->next_bit) {
X termptr = (struct octree_terminal *)&(tptr->OS_space[0]);
X for(i=0;i<=tptr->OS_last_used;i++) {
X for(j=0;j<8;j++){
X if(group_no == termptr[i].OT_term_assign[j]) {
X if(termptr[i].OT_term_cnt[j]>0) {
X rval = (termptr[i].OT_red_val + red_values[j])*red_weight;
X gval = (termptr[i].OT_green_val + green_values[j])*green_weight;
X bval = (termptr[i].OT_blue_val + blue_values[j])*blue_weight;
X mass = termptr[i].OT_term_cnt[j];
X data_ptr->xavg += (double)(rval*mass);
X data_ptr->yavg += (double)(gval*mass);
X data_ptr->zavg += (double)(bval*mass);
X popcnt += mass;
X }
X }
X }
X }
X }
X }
X else {
X for(ppptr=color_levels[group_no];ppptr;ppptr=ppptr->next_data) {
X rval = (ppptr->xd) * red_weight;
X gval = (ppptr->yd) * green_weight;
X bval = (ppptr->zd) * blue_weight;
X mass = (ppptr->pixelcnt);
X data_ptr->xavg += (double)(rval*mass);
X data_ptr->yavg += (double)(gval*mass);
X data_ptr->zavg += (double)(bval*mass);
X popcnt += mass;
X }
X }
X if(popcnt > 0) {
X data_ptr->xavg /= (double)(popcnt);
X data_ptr->yavg /= (double)(popcnt);
X data_ptr->zavg /= (double)(popcnt);
X }
X#ifdef DEBUG
X printf("For group %d, we have average ( %g, %g, %g )\n", group_no,
X data_ptr->xavg, data_ptr->yavg, data_ptr->zavg );
X#endif
X xmid = data_ptr->xavg;
X ymid = data_ptr->yavg;
X zmid = data_ptr->zavg;
X data_ptr->a = 0.0;
X data_ptr->b = 0.0;
X data_ptr->c = 0.0;
X data_ptr->f = 0.0;
X data_ptr->g = 0.0;
X data_ptr->h = 0.0;
X if( line_data == NULL) {
X for(tptr=octree_levels[gl_current_nbits - 1];tptr;tptr=tptr->next_bit) {
X termptr = (struct octree_terminal *)&(tptr->OS_space[0]);
X for(i=0;i<=tptr->OS_last_used;i++) {
X for(j=0;j<8;j++){
X if(group_no == termptr[i].OT_term_assign[j]) {
X if(termptr[i].OT_term_cnt[j]>0) {
X rval = (termptr[i].OT_red_val + red_values[j])*red_weight - xmid;
X gval = (termptr[i].OT_green_val + green_values[j])*green_weight - ymid;
X bval = (termptr[i].OT_blue_val + blue_values[j])*blue_weight - zmid;
X mass = termptr[i].OT_term_cnt[j];
X data_ptr->a += mass * (gval*gval+bval*bval);
X data_ptr->b += mass * (bval*bval+rval*rval);
X data_ptr->c += mass * (rval*rval+gval*gval);
X data_ptr->f += mass * gval * bval;
X data_ptr->g += mass * bval * rval;
X data_ptr->h += mass * rval * gval;
X }
X }
X }
X }
X }
X }
X else {
X for(ppptr=color_levels[group_no];ppptr;ppptr=ppptr->next_data) {
X rval = (ppptr->xd) * red_weight - xmid;
X gval = (ppptr->yd) * green_weight - ymid;
X bval = (ppptr->zd) * blue_weight - zmid;
X mass = (ppptr->pixelcnt);
X data_ptr->a += mass * (gval*gval+bval*bval);
X data_ptr->b += mass * (bval*bval+rval*rval);
X data_ptr->c += mass * (rval*rval+gval*gval);
X data_ptr->f += mass * gval * bval;
X data_ptr->g += mass * bval * rval;
X data_ptr->h += mass * rval * gval;
X }
X }
X
X data_ptr->sigmasquare = ( data_ptr->a + data_ptr->b + data_ptr->c ) / 2.0;
X inertia(data_ptr->a, data_ptr->b, data_ptr->c,
X data_ptr->f, data_ptr->g, data_ptr->h,
X &(data_ptr->xdot), &(data_ptr->ydot), &(data_ptr->zdot),
X &(data_ptr->maxmoment));
X } /* sum_group */
X
X/* ---------------------------------------------------------------------
XCount the number of non-zero entries in the octree
X--------------------------------------------------------------------- */
X#ifdef DEBUG
X#undef DHEADER
X#define DHEADER printf("\ncount_distinct_colors-");
X#endif
Xlong
Xcount_distinct_colors( )
X{
X struct octree_space *tptr;
X struct octree_terminal *termptr;
X long i, j, cnt;
X
X cnt = 0;
X for(tptr=octree_levels[gl_current_nbits-1];tptr;tptr = tptr->next_bit) {
X termptr = (struct octree_terminal *)&(tptr->OS_space[0]);
X for(i=0;i<=tptr->OS_last_used;i++) {
X for(j=0;j<8;j++) {
X if ( termptr[i].OT_term_cnt[ j ] > 0 ) cnt++;
X }
X }
X }
X return( cnt );
X } /* count_distinct_colors */
X
X/* ----------------------------------------------------------------------
XThis records the data at the terminal level of the octree in the array
Xline_data. Builds the linked list of cells.
X---------------------------------------------------------------------- */
X#ifdef DEBUG
X#undef DHEADER
X#define DHEADER printf("\nrestructure_info-");
X#endif
Xrestructure_info( )
X{
X struct octree_space *tptr;
X struct octree_terminal *termptr;
X double rval, gval, bval, sum;
X long i, j;
X
X long num_data;
X
X num_data = 0;
X for(tptr=octree_levels[gl_current_nbits-1];tptr;tptr = tptr->next_bit) {
X termptr = (struct octree_terminal *)&(tptr->OS_space[0]);
X for(i=0;i<=tptr->OS_last_used;i++) {
X for(j=0;j<8;j++) {
X if ( termptr[i].OT_term_cnt[j] > 0) {
X line_data[num_data].xd = termptr[i].OT_red_val + red_values[j];
X line_data[num_data].yd = termptr[i].OT_green_val + green_values[j];
X line_data[num_data].zd = termptr[i].OT_blue_val + blue_values[j];
X line_data[num_data].pixelcnt = termptr[i].OT_term_cnt[j];
X line_data[num_data].next_data = NULL;
X if(num_data>0) {
X line_data[num_data-1].next_data = &(line_data[num_data]);
X }
X num_data++;
X }
X }
X }
X }
X color_levels[ 0 ] = &(line_data[0]);
X }
X
X
X/* ---------------------------------------------------------------------
XSupport procedures for octree representation. We have
X
X initialize_octree - sets up initial pointers, and calculates the
X needed values.
X
X increment_cell(R,G,B,val) - Procedure to increment the count of a cell
X of the octree
X
X clean_up_octree - frees data spaces and sweeps up
X
XINTERNAL PROCEDURES (used only by support procedures)
X
X allocate_data_space - mallocs and initializes a block of space.
X gets something off the free list if it is
X there.
X
X reduce_octree - gets rid of one level of the tree (sums up
X the data at the following level and makes
X new level the terminal level)
X
XDEBUG PROCEDURES
X
X put_location - print pointers in readable format.
X
X--------------------------------------------------------------------- */
X
X/* Simple greatest common denominator algorithm */
X#ifdef DEBUG
X#undef DHEADER
X#define DHEADER printf("\ngcd-");
X#endif
X#ifdef __STDC__
Xstatic long
Xgcd( long i, long j )
X#else /* __STDC__ */
Xstatic long
Xgcd( i, j )
Xlong i;
Xlong j;
X#endif /* __STDC__ */
X{
X long ibig, ilittle, k;
X
X if ( i < j ) {
X ilittle = i;
X ibig = j;
X }
X else {
X ilittle = j;
X ibig = i;
X }
X
X if( ilittle == 0 ) {
X return( 0 );
X }
X k = ibig % ilittle;
X if ( k == 0 ) {
X return( ilittle );
X }
X else {
X return( gcd( ilittle, k ) );
X }
X }
X
X#ifdef DEBUG
X#undef DHEADER
X#define DHEADER printf("\nallocate_data_space-")
X#endif
X/* ---------------------------------------------------------------------
XAllocate a new bit of space using malloc. set up the data structure,
Xdo not initialize the space.
X
XIf a reduction in octree representation was required, return NULL
X--------------------------------------------------------------------- */
X#ifdef __STDC__
Xstatic struct octree_space *
Xallocate_data_space( unsigned long element_size )
X#else /* __STDC__ */
Xstatic struct octree_space *
Xallocate_data_space( element_size )
Xunsigned long element_size;
X#endif /* __STDC__ */
X{
X struct octree_space *temp_ptr = NULL;
X long header_size = sizeof(struct octree_space ) - sizeof( UBYTE );
X
X /* We have to be sure we keep enough stuff in the free_space list to reduce the */
X /* outermost non-terminal level. */
X DEBUGMSG(10,("called, element_size=%d, header_size=%d", element_size, header_size) );
X if ( free_space ) {
X DEBUGMSG(10,("Getting space off the free list"));
X temp_ptr = free_space;
X free_space = temp_ptr->next_bit;
X }
X else {
X DEBUGMSG(10,("Attempting to malloc some space"));
X if ( !(temp_ptr = (struct octree_space *)malloc( header_size + OS_block_size ))) {
X DEBUGMSG(10,("Could not malloc space; must reduce"));
X if (gl_current_nbits > 3) {
X DEBUGMSG(10,("Returning NULL; let someone else worry" ));
X return( NULL );
X }
X else {
X printf( "Could not malloc space.\n" );
X exit(1);
X }
X }
X }
X
X /* now we have space pointed to by temp_ptr */
X
X DEBUGMSG(10,("We have pointer %d\n", temp_ptr ) );
X temp_ptr->OS_size = element_size;
X temp_ptr->OS_length = OS_block_size / element_size;
X temp_ptr->OS_last_used = -1; /* none used yet. */
X DEBUGMSG(100,("Setting OS_size=%d, OS_length=%d, OS_last_used=%d", temp_ptr->OS_size, temp_ptr->OS_length, temp_ptr->OS_last_used ));
X temp_ptr->next_bit = NULL;
X
X return( temp_ptr );
X } /* allocate_data_space */
X
X#ifdef DEBUG
X#undef DHEADER
X#define DHEADER printf("\nnew_record_at_level-")
X#endif
X#ifdef __STDC__
Xstatic VOIDPTR
Xnew_record_at_level( long at_level )
X#else /* __STDC__ */
Xstatic VOIDPTR
Xnew_record_at_level( at_level )
Xlong at_level;
X#endif /* __STDC__ */
X{
X long i;
X struct octree_space *temp_ptr, *new_ptr;
X VOIDPTR new_rec;
X struct octree_nonterminal *ntermptr;
X struct octree_terminal *termptr;
X
X DEBUGMSG(10,("at_level=%d", at_level));
X temp_ptr = octree_levels[ at_level ];
X if ( temp_ptr->OS_length - 1 == temp_ptr->OS_last_used ) {
X DEBUGMSG(10,("No more room here; try to get more" ));
X if(!(new_ptr = allocate_data_space( temp_ptr->OS_size ))) {
X DEBUGMSG(10,("Could not allocate; return NULL" ));
X return(NULL);
X }
X new_ptr->next_bit = temp_ptr;
X octree_levels[ at_level ] = new_ptr;
X temp_ptr = new_ptr;
X }
X
X temp_ptr->OS_last_used += 1;
X DEBUGMSG(10,("Add record %d", temp_ptr->OS_last_used));
X new_rec = (VOIDPTR)(&(temp_ptr->OS_space[ temp_ptr->OS_size * temp_ptr->OS_last_used ]));
X if (at_level < gl_current_nbits - 1) {
X DEBUGMSG(10,("Non-terminal record" ));
X ntermptr = (struct octree_nonterminal *)new_rec;
X for(i=0;i<8;i++) {
X ntermptr->next_level[ i ] = NULL;
X }
X }
X else {
X DEBUGMSG(10,("Terminal record"));
X termptr = (struct octree_terminal *)new_rec;
X for(i=0;i<8;i++) {
X termptr->OT_term_cnt[ i ] = 0;
X termptr->OT_term_assign[ i ] = 0;
X }
X }
X return ( new_rec );
X } /* new_record_at_level */
X
X#ifdef DEBUG
X#undef DHEADER
X#define DHEADER printf("\nreduce_octree-");
X#endif
XVOID
Xreduce_octree( )
X{
X struct octree_space *nonterminal_data, *terminal_data;
X struct octree_space *tptr;
X struct octree_nonterminal *nptr, *nnptr;
X struct octree_terminal *optr, *new_terminal_ptr;
X long i, j, k, gbitmask, summa;
X char tostring[100];
X
X pm_message( "Reducing octree to %d bits", gl_current_nbits - 1 );
X nonterminal_data = octree_levels[ gl_current_nbits - 2 ];
X terminal_data = octree_levels[ gl_current_nbits - 1 ];
X gl_current_nbits -= 1; /* reduce representation by one bit */
X gbitmask = !((1<<(gl_bits_per_color - gl_current_nbits + 1)) - 1);
X DEBUGMSG(100,("gl_current_nbits=%d, gbitmask=%x", gl_current_nbits, gbitmask ));
X octree_levels[ gl_current_nbits - 1 ] = allocate_data_space( (unsigned long) sizeof(struct octree_terminal) );
X
X /* First, we will sum the terminal level of data, and coerce that value to
X type (struct octree_nonterminal *) and store in the pointer next_level.
X Then, we can put the terminal level on the free list. At this point we
X know we have enough space on the free list to represent the new terminal
X level as terminal records, not non-terminal records.
X
X This has the problem we lose the (R,G,B) values in the terminal nodes.
X We can recreate that data from the tree structure, and will do so later.
X On the other hand, this has the advantage we need not insure we
X retain enough free space to represent the new terminal level. */
X
X for(tptr = nonterminal_data;tptr;tptr=tptr->next_bit) {
X nptr = (struct octree_nonterminal *)&(tptr->OS_space[0]);
X for(i=0;i<=tptr->OS_last_used;i++) {
X for(j=0;j<8;j++) {
X summa = 0;
X optr = (struct octree_terminal *)(nptr[i].next_level[j]);
X if(optr) {
X for(k=0;k<8;k++) {
X summa += optr->OT_term_cnt[k];
X }
X nptr[i].next_level[j] = (struct octree_nonterminal *)(summa);
X }
X else {
X nptr[i].next_level[j] = (struct octree_nonterminal *)0;
X }
X }
X }
X }
X
X /* Now, we can place the terminal level on the free list. */
X
X if(free_space) {
X for (tptr=terminal_data;tptr->next_bit;tptr=tptr->next_bit);
X tptr->next_bit = free_space;
X }
X free_space = terminal_data;
X
X /* now, we can replace the bogus "terminal" records with real terminal records. */
X for(tptr = octree_levels[ gl_current_nbits - 2 ];tptr;tptr = tptr->next_bit) {
X /* we will walk each data structure, treating each non-terminal; adding up the */
X /* sums at the terminal level, making new record for next level */
X nptr = (struct octree_nonterminal *)&(tptr->OS_space[0]);
X for(i=0;i<=tptr->OS_last_used;i++) {
X /* nptr[i] is the record we are looking at. */
X for(j=0;j<8;j++){
X /* nptr[i].next_level[ j ] is what we are reducing now */
X nnptr = nptr[i].next_level[j];
X if (nnptr) {
X new_terminal_ptr = (struct octree_terminal *)new_record_at_level( gl_current_nbits - 1 );
X for(k=0;k<8;k++){
X new_terminal_ptr->OT_term_cnt[k] =
X (long)(nnptr->next_level[k]);
X }
X nptr[i].next_level[j] = (struct octree_nonterminal *)new_terminal_ptr;
X }
X }
X }
X }
X for(tptr=nonterminal_data;tptr->next_bit;tptr = tptr->next_bit);
X tptr->next_bit = free_space;
X free_space = nonterminal_data;
X } /* reduce_octree */
X
X#ifdef DEBUG
X#undef DHEADER
X#define DHEADER printf("\ninitialize_octree-");
X#endif
X#ifdef __STDC__
XVOID
Xinitialize_octree( long nbits_per_color, long nbits_per_side, long num_of_registers )
X#else /* __STDC__ */
XVOID
Xinitialize_octree( nbits_per_color, nbits_per_side, num_of_registers )
Xlong nbits_per_color, nbits_per_side, num_of_registers;
X#endif /* __STDC__ */
X{
X long i, j, k;
X
X DEBUGMSG(10,("nbits_per_color=%d, nbits_per_side=%d", nbits_per_color, nbits_per_side));
X if(nbits_per_color<3)
X pm_error( "Too few bits for octree" );
X
X if(nbits_per_color > COLOR_REP_BITS)
X pm_error( "Too many bits for octree" );
X
X if(nbits_per_side>0) {
X if(nbits_per_side < 3) {
X pm_error( "Cannot do fewer than 3 bits" );
X }
X }
X if ( OS_block_size == 0 ) {
X /* Make OS_block_size divisible by either the nonterminal or terminal size */
X /* and approximately equal to APPROXIMATE_BLOCK_SIZE */
X i = sizeof(struct octree_nonterminal );
X j = sizeof(struct octree_terminal );
X k = gcd( i,j );
X i = i * j / k; /* LCM */
X j = APPROXIMATE_BLOCK_SIZE % i;
X if ( j < (i / 2)) {
X OS_block_size = APPROXIMATE_BLOCK_SIZE - j;
X }
X else {
X OS_block_size = APPROXIMATE_BLOCK_SIZE + i - j;
X }
X }
X
X /* Initialize the data structures for creating the cube */
X
X DEBUGMSG(10,("We now have OS_block_size=%d", OS_block_size));
X nregs = num_of_registers;
X gl_bits_per_color = nbits_per_color;
X if ( nbits_per_side > 0 ) {
X gl_current_nbits = nbits_per_side;
X }
X else {
X gl_current_nbits = nbits_per_color;
X }
X
X /* Allocate a chunk of space for each needed level, so we do not have to
X make that test when traversing the tree */
X pm_message( "Initializing octree data structure" );
X for(i=0;i<gl_current_nbits;i++) {
X if ( i<2 ) {
X /* Don't allocate for the first two levels; we use static data for them */
X octree_levels[ i ] = NULL;
X }
X else if (i < (gl_current_nbits - 1) ) {
X octree_levels[ i ] = allocate_data_space( (unsigned long)sizeof(struct octree_nonterminal ) );
X }
X else {
X octree_levels[ i ] = allocate_data_space( (unsigned long)sizeof(struct octree_terminal ) );
X }
X }
X
X /* Set up the first two levels of the tree, so we don't have to */
X /* do that part. */
X for(i=0;i<8;i++) {
X first_level.next_level[ i ] = NULL;
X }
X second_level.OS_length = 8;
X second_level.OS_size = sizeof(struct octree_nonterminal);
X second_level.OS_last_used = -1;
X octree_levels[ 1 ] = (struct octree_space *)&second_level;
X
X /* bit 0 ==> first_level
X bit 1 ==> second_level
X bit i (i>1) ==> octree_levels[ i ] */
X
X for(i=0;i<8;i++) {
X blue_values[i] = i & 1;
X green_values[i] = (i & 2) >>1;
X red_values[i] = (i & 4) >>2;
X }
X
X picturesize = 0;
X } /* initialize_octree */
X
X#ifdef DEBUG
X#undef DHEADER
X#define DHEADER printf("\nincrement_cell-");
X#endif
X/* ---------------------------------------------------------------------
XThis is the workhorse procedure in population counts. If there is a
Xproblem in space, the procedure is called recursively after the size
Xof the cube has been reduced.
X--------------------------------------------------------------------- */
X#ifdef __STDC__
XVOID
Xincrement_cell( long red_val, long green_val, long blue_val, long inc_value )
X#else /* __STDC__ */
XVOID
Xincrement_cell( red_val, green_val, blue_val, inc_value )
Xlong red_val, green_val, blue_val, inc_value;
X#endif /* __STDC__ */
X{
X long i, ii, j, rbit, gbit, bbit, indx, gbitmask;
X struct octree_nonterminal *current_pt = &first_level;
X struct octree_terminal *final_pt, *new_termptr;
X struct octree_nonterminal *new_nonptr;
X
X DEBUGMSG(100,("red_val=%d, green_val=%d, blue_val=%d, inc_value=%d", red_val, green_val, blue_val, inc_value));
X DEBUGMSG(100,("gl_bits_per_color=%d, gl_current_nbits=%d", gl_bits_per_color, gl_current_nbits));
X for(ii=0,i=gl_bits_per_color-1;i>= gl_bits_per_color - gl_current_nbits;ii++, i--) {
X rbit = (red_val>>i)& 1;
X gbit = (green_val>>i) & 1;
X bbit = (blue_val>>i) & 1;
X DEBUGMSG(100,("level %d(%d) : (%d,%d,%d)", ii, i, rbit, gbit, bbit));
X DEBUGMSG(999,("gl_bits_per_color - gl_current_nbits=%d", gl_bits_per_color - gl_current_nbits));
X indx = (rbit<<2) + (gbit<<1) + bbit;
X DEBUGMSG(100,("indx=%d", indx));
X
X if ( i > gl_bits_per_color - gl_current_nbits) {
X DEBUGMSG(100,("Non-terminal level"));
X if ( current_pt->next_level[ indx ] ) {
X DEBUGMSG(100,("Existing record"));
X current_pt = current_pt->next_level[ indx ];
X }
X else { /* missing level */
X DEBUGMSG(100,("Missing record"));
X if ( i > gl_bits_per_color - gl_current_nbits + 1) {
X DEBUGMSG(100,("Allocating new non-terminal at level %d", ii + 1));
X if (!(new_nonptr = (struct octree_nonterminal *)new_record_at_level( ii+1 ))) {
X DEBUGMSG(100,("Have to reduce tree and try again"));
X reduce_octree( );
X increment_cell( red_val, green_val, blue_val, inc_value );
X return;
X }
X current_pt->next_level[ indx ] = new_nonptr;
X current_pt = new_nonptr;
X }
X else {
X DEBUGMSG(100,("Allocating new terminal at level %d", ii + 1 ));
X if (!(new_termptr = ( struct octree_terminal *)new_record_at_level( ii+1 ))) {
X DEBUGMSG(100,("Reduce and try again"));
X reduce_octree( );
X increment_cell( red_val, green_val, blue_val, inc_value );
X return;
X }
X current_pt->next_level[ indx ] = (struct octree_nonterminal *)new_termptr;
X current_pt = (struct octree_nonterminal *)new_termptr;
X DEBUGMSG(100,("Have to initialize record"));
X for(j=0;j<8;j++) {
X new_termptr->OT_term_cnt[ j ] = 0;
X }
X gbitmask = !((1<<(gl_bits_per_color - gl_current_nbits + 1)) - 1);
X new_termptr->OT_red_val = red_val & gbitmask;
X new_termptr->OT_green_val = green_val & gbitmask;
X new_termptr->OT_blue_val = blue_val & gbitmask;
X }
X }
X }
X else {
X DEBUGMSG(100,("At final level; put in bucket"));
X final_pt = (struct octree_terminal *)current_pt;
X final_pt->OT_term_cnt[ indx ] += inc_value;
X picturesize++;
X return;
X }
X }
X } /* increment_cell */
X
X#ifdef DEBUG
X#undef DHEADER
X#define DHEADER printf("\nassign_at_level-");
X#endif
X#ifdef __STDC__
Xstatic VOID
Xassign_at_level( long red_in, long green_in, long blue_in, long level_in, VOIDPTR ptr )
X#else /* __STDC__ */
Xstatic VOID
Xassign_at_level( red_in, green_in, blue_in, level_in, ptr )
Xlong red_in, green_in, blue_in, level_in;
XVOIDPTR ptr;
X#endif
X{
X struct octree_nonterminal *nptr;
X struct octree_terminal *tptr;
X long i;
X
X red_in = red_in << 1;
X green_in = green_in << 1;
X blue_in = blue_in << 1;
X
X if (level_in < gl_current_nbits - 1) {
X nptr = (struct octree_nonterminal *)ptr;
X for(i=0;i<8;i++) {
X if (nptr->next_level[i]) {
X assign_at_level( red_in + red_values[i], green_in + green_values[i], blue_in + blue_values[i],
X level_in + 1, nptr->next_level[i] );
X }
X }
X }
X else {
X tptr = (struct octree_terminal *)ptr;
X
X i = gl_bits_per_color - gl_current_nbits;
X tptr->OT_red_val = red_in << i;
X tptr->OT_green_val = green_in << i;
X tptr->OT_blue_val = blue_in << i;
X }
X } /* assign_at_level */
X
X#ifdef DEBUG
X#undef DHEADER
X#define DHEADER printf("\nassign_colors-");
X#endif
XVOID
Xassign_colors( )
X{
X DEBUGMSG(10,("gl_current_nbits=%d", gl_current_nbits ));
X assign_at_level( 0, 0, 0, 0, (VOIDPTR)&(first_level));
X } /* assign_colors */
X
X/* ---------------------------------------------------------------------
XFree up the space occupied by the octree. If all_but_last is TRUE, will
Xleave that level. Can be called repeatedly; it is nice about how it
Xleaves the pointers.
X--------------------------------------------------------------------- */
X#ifdef DEBUG
X#undef DHEADER
X#define DHEADER printf("\nclean_up_octree-");
X#endif
XVOID
X#ifdef __STDC__
Xclean_up_octree( long all_but_last )
X#else /* __STDC__ */
Xclean_up_octree( all_but_last )
Xlong all_but_last;
X#endif /* __STDC__ */
X{
X struct octree_space *tptr, *ntptr;
X long i, last_lev;
X
X last_lev = ( all_but_last ) ? gl_current_nbits - 1 : gl_current_nbits;
X
X for(tptr=free_space;tptr;tptr=ntptr) {
X ntptr = tptr->next_bit;
X free(tptr);
X }
X free_space = NULL;
X /* We have to start at 2 because level 1 is a global structure */
X for(i=2;i<last_lev;i++) {
X for(tptr=octree_levels[i];tptr;tptr = ntptr) {
X ntptr = tptr->next_bit;
X free(tptr);
X }
X octree_levels[i] = NULL;
X }
X }
X
X#ifdef DEBUG
X/* ---------------------------------------------------------------------
XThis procedure will print a more readable version of pointers into the
Xoctree. The form will be
X
X Lxx.yyy.zzz where
X
X xx - index into octree_levels the pointer is at
X yyy - ordinal in the list of the structure the pointer is in
X zzz - index from the start of the pointer within the space.
X
XSo, if the pointer is in octree_level[1], in the second dataspace in
Xthat list, and the third record in that dataspace, we should get
X
X L01.001.002
X
XThe freespace pointers will be printed as
X
X FSyyy.zzz
X
Xwhere "yyy" and "zzz" are as above.
X
XIf we cannot find the pointer in the known space, we will print the
Xaddress itself.
X--------------------------------------------------------------------- */
X#ifdef DEBUG
X#undef DHEADER
X#define DHEADER printf("\nput_location-");
X#endif
XVOID
Xput_location( VOIDPTR ptr )
X{
X long i,j,k, lowaddr, highaddr, blksize;
X long longptr;
X struct octree_space *tptr;
X
X DEBUGMSG(9999,("ptr=%d", (long)ptr));
X longptr = (long)ptr;
X if(ptr == NULL) {
X printf( " <NULL> " );
X return;
X }
X for(i=0;i<gl_current_nbits;i++) {
X tptr = octree_levels[ i ];
X j = 0;
X for(j=0,tptr=octree_levels[ i ]; tptr; j++, tptr=tptr->next_bit) {
X lowaddr = (long)&(tptr->OS_space[0]);
X blksize = tptr->OS_size;
X highaddr = lowaddr + (tptr->OS_length)*blksize - 1;
X DEBUGMSG(9999,("lowaddr=%d, highaddr=%d, blksize=%d", lowaddr, highaddr, blksize ));
X if(((long)ptr>=lowaddr)&&((long)ptr<highaddr)) {
X k = (long)ptr - lowaddr;
X if(k%blksize) {
X printf("<er>" );
X };
X k /= blksize;
X printf("L%02d\.%03d\.%03d", i, j, k );
X return;
X }
X }
X }
X for(i=0, tptr = free_space;tptr;i++,tptr=tptr->next_bit) {
X lowaddr = (long)&(tptr->OS_space[0]);
X highaddr = lowaddr + OS_block_size;
X if(((long)ptr>=lowaddr)&&((long)ptr<highaddr)) {
X j = (long)ptr - lowaddr;
X if(!(j%(tptr->OS_size))) {
X printf("<er>");
X };
X j /= tptr->OS_size;
X printf("FS%03d\.%03d", i, j );
X return;
X }
X }
X /* We could not find the pointer in any known space. */
X printf("<<%d>>", (long)ptr);
X }
X
X#ifdef DEBUG
X#undef DHEADER
X#define DHEADER printf("\ndump_nonterminal-");
X#endif
XVOID
Xdump_nonterminal(struct octree_nonterminal *nptr, long array_len )
X{
X long k, l;
X for(k=0;k<=array_len;k++) {
X printf( " -- RECORD %d ", k );
X for(l=0;l<8;l++) {
X put_location( nptr[k].next_level[l] );
X printf( " " );
X }
X printf( "\n" );
X }
X }
X
X#ifdef DEBUG
X#undef DHEADER
X#define DHEADER printf("\ndump_terminal-");
X#endif
XVOID
Xdump_terminal(struct octree_terminal *tptr, long array_len )
X{
X long k, l;
X for(k=0;k<=array_len;k++) {
X printf( " -- RECORD %d (%2x:%2x:%2x) -:-", k, tptr[k].OT_red_val, tptr[k].OT_green_val, tptr[k].OT_blue_val );
X for(l=0;l<8;l++) {
X printf( " %d:%7d", l, tptr[k].OT_term_cnt[l] );
X }
X printf( "\n" );
X }
X }
X
X#ifdef DEBUG
X#undef DHEADER
X#define DHEADER printf("\ndump_tree-");
X#endif
XVOID
Xdump_tree()
X{
X long i,j,k, l;
X long array_len, low_addr, high_addr, last_addr;
X struct octree_nonterminal *ntermptr;
X struct octree_terminal *termptr;
X struct octree_space *tptr;
X
X printf("\lDump of Octree data\n");
X printf(" gl_current_nbits=%d, gl_bits_per_color=%d\n", gl_current_nbits, gl_bits_per_color );
X printf(" OS_block_size=%d sizeof(octree_nonterminal)=%d, sizeof(octree_terminal)=%d\n",
X OS_block_size,
X sizeof(struct octree_nonterminal),
X sizeof(struct octree_terminal) );
X for(i=0;i<gl_current_nbits;i++) {
X printf( "\n\noctree_level[%02d] ", i );
X if(i==gl_current_nbits - 1) {
X printf( "(TERMINAL)" );
X }
X printf( ":\n" );
X if(!octree_levels[i]) {
X if(i==0) {
X dump_nonterminal( &first_level, 0);
X }
X else {
X printf( " EMPTY\n" );
X }
X }
X else {
X for(j=0,tptr=octree_levels[i];tptr;j++,tptr=tptr->next_bit) {
X array_len = tptr->OS_last_used;
X low_addr = (long)&(tptr->OS_space[0]);
X high_addr = low_addr + (tptr->OS_length)*(tptr->OS_size) - 1;
X last_addr = low_addr + (tptr->OS_size)*(tptr->OS_last_used+1) - 1;
X printf( " -- STRUCTURE %d (OS_length=%d) (OS_last_used=%d) (OS_size=%d) low:%d last:%d high:%d\n", j,
X tptr->OS_length,
X array_len,
X tptr->OS_size, low_addr, last_addr, high_addr );
X if(i==gl_current_nbits-1) {
X /* dump terminal records */
X termptr = (struct octree_terminal *)&(tptr->OS_space[0]);
X dump_terminal( termptr, array_len );
X }
X else {
X /* dump nonterminal records */
X ntermptr = (struct octree_nonterminal *)&(tptr->OS_space[0]);
X dump_nonterminal( ntermptr, array_len );
X }
X }
X }
X }
X printf( "\n\n\n" );
X }
X#endif
X
X#define EPSILON 1e-6
X/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Xgiven the values
X 2 2 2 2
X a = sigma( m * (y + z )) b = sigma( m * (z + x ))
X
X 2 2
X c = sigma( m * (x + y )) f = sigma( myz )
X
X g = sigma( mzx ) h = sigma( mxy )
X
X(where m is the "mass" of the point, in our case the population
Xcount or some analog thereof.)
X
Xfinds the axis values (x, y, z) of the axis along which the
Xset of points has the minimal moment of inertia. When the
Xset of points is split along the plane this axis is perpendicular
Xto, the resulting two sets of points have the maximal distance
Xbetween them, and the maximal reduction in the error in register
Xassignment.
X
X This routine implements the solution of the matrix
X
X | A - K -H -G |
X | |
X | -H B - K -F | = 0
X | |
X | -G -F C - K |
X
Xsolving for the three roots K, substituting the values back to
Xfind the three vectors which are the principal axes of inertia
Xof the set of points being examined. The minimal axis is the
Xone returned.

X
XCopyright 1994 Christopher A. Huson

X- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
XVOID
Xinertia(a,b,c,f,g,h, xaxis, yaxis, zaxis, maxmom )
Xdouble a, b, c, f, g, h;
Xdouble *xaxis, *yaxis, *zaxis, *maxmom;
X{
X#define NEAR_ZERO 0.0001
X double p,q,r, r2i, r3i, ap, bp, cp;
X double x12, y12, z12, x13, y13, x23, vmax;
X double xsq, ysq, zsq;
X double newx, newy, newz;
X int i, j, icase, solve_for_max;
X double roots[ 3 ];
X#ifdef DEBUG
X fprintf( stderr, "a=%g, b=%g, c=%g, f=%g, g=%g, h=%g\n", a, b, c, f, g, h );
X#endif
X if((a+b+c) <= NEAR_ZERO) {
X#ifdef DEBUG
X fprintf( stderr, "all values are close to zero; choose a vector arbitrarily\n" );
X#endif
X *xaxis = 1.0;
X *yaxis = 0.0;
X *zaxis = 0.0;
X *maxmom = 0.0;
X return;
X };
X
X p = - ( a + b + c );
X q = (b*c + a*c + a*b - f*f - g*g - h*h);
X r = 2*f*g*h + g*g*b + h*h*c + f*f*a - a*b*c;
X
X cubic_roots( p, q, r, &roots[0], &roots[1], &r2i, &roots[2], &r3i );
X
X#ifdef DEBUG
X if( r2i != 0e0) {
X fprintf( stderr, "Non-zero imaginary parts of roots %g\n", r2i );
X };
X#endif
X /* sort the roots by increasing size */
X for(i=0;i<3;i++) {
X for(j=i+1;j<3;j++) {
X if(roots[i] > roots[j]) {
X r = roots[i];
X roots[i] = roots[j];
X roots[j] = r;
X };
X }
X }
X#ifdef DEBUG
X fprintf( stderr, "Moments about principal axes are %g, %g and %g\n",
X roots[0], roots[1], roots[2] );
X#endif
X
X *maxmom = roots[2];
X
X /*
X Remember, we have purposely set the determinant to 0, which is
X not soluble (that is, there is not one single solution). If we assume
X one value is given in terms of the other two, then by a combination of
X any two of the three equations, and introducing a change of variables,
X we can solve the other two values (as long as one of the minors of the
X matrix has a non-zero determinant).
X
X There are (I believe) three possibilities :
X
X 1) The minimum moment of inertia is 0 (or very close). In this case,
X the sample is collinear; we can discover the vector along the
X line by using the identities
X
X A = sigma(m(y**2 + z**2))
X = sigma(m(ry**2 * t**2 + rz**2 * t**2)) (all are collinear)
X = sigma(m*t**2) (ry**2 + rz**2)
X
X we can decompose B and C the same way, which results in
X
X rx**2 = (C - A + B) / 2 * sigma(m*t**2)
X ry**2 = (A - B + C) / 2 * sigma(m*t**2)
X rz**2 = (B - C + A) / 2 * sigma(m*t**2)
X
X and we have
X
X F = sigma(m*y*z)
X = sigma(m*ry*t*rz*t)
X = sigma(m*t**2)*ry*rz
X
X Similarly for G and H. So, we can get the absolute values of
X a ratio for rx, ry and rz using the A, B and C, and can find a
X set of signs for the vector using F, G and H.
X
X 2) Of the three moments of inertia, the smaller two are equal. This
X indicates symmetry at the minimal moment of inertia. To get a minimal
X vector, solve for the maximal, and choose a value for one of
X X, Y or Z, and then solve the equation
X
X dot_product( _MAX_, _vec_ ) = 0
X
X where _MAX_ is the maximal vector.
X
X 3) only one minimal moment of inertia. Just solve for the vector
X below.
X
X Now, we have to give ourselves the best chance of coming up with the
X most accurate axes. There are only six distinct values resulting from
X the above minors. If we number them {x,y,z}ab, where "a" is the first
X row in the solution, and "b" is the second, then we have (given the signs
X of the solutions below):
X
X x12 = -z23 y12 = z13 x13 = -y23
X
X We get that index with the largest absolute value, and figure the final
X axes given that case (stored in icase).
X */
X
X if( fabs( roots[0] ) <= NEAR_ZERO ) {
X#ifdef DEBUG
X fprintf( stderr, "smallest root is near zero\n" );
X#endif
X xsq = (c - a + b);
X ysq = (a - b + c);
X zsq = (b - c + a);
X#ifdef DEBUG
X fprintf( stderr, "xsq=%g, ysq=%g, zsq=%g\n", xsq, ysq, zsq );
X#endif
X xsq = sqrt(fabs(xsq));
X ysq = sqrt(fabs(ysq));
X zsq = sqrt(fabs(zsq));
X
X if(ysq+zsq>0) {
X if(xsq/(ysq+zsq) < EPSILON) {
X xsq = 0.0;
X }
X }
X
X if(xsq+zsq>0) {
X if(ysq/(xsq+zsq) < EPSILON) {
X ysq = 0.0;
X }
X }
X
X if(xsq+ysq>0) {
X if(zsq/(xsq+ysq) < EPSILON) {
X zsq = 0.0;
X }
X }
X /* we will choose x to be positive */
X if(xsq > 0.0) {
X if(g < 0.0) {
X zsq = -zsq;
X };
X if(h < 0.0) {
X ysq = -ysq;
X };
X }
X else if(ysq > 0.0) {
X if(f < 0.0) {
X zsq = -zsq;
X };
X };
X
X#ifdef DEBUG
X fprintf( stderr, "xsq=%g, ysq=%g, zsq=%g\n", xsq, ysq, zsq );
X#endif
X *xaxis = xsq;
X *yaxis = ysq;
X *zaxis = zsq;
X return;
X };
X
X solve_for_max = ( fabs(roots[0] - roots[1]) <= NEAR_ZERO );
X
X if(solve_for_max) {
X#ifdef DEBUG
X fprintf( stderr, "Solving for maximum rather than minimum\n" );
X#endif
X ap = a - roots[2]; bp = b - roots[2]; cp = c - roots[2];
X }
X else {
X#ifdef DEBUG
X fprintf( stderr, "Solving for minimum\n" );
X#endif
X ap = a - roots[0]; bp = b - roots[0]; cp = c - roots[0];
X };
X
X x12 = h*f + g*bp;
X y12 = ap*f + g*h;
X z12 = ap*bp - h*h;
X x13 = (-h)*cp - f*g;
X y13 = ap*cp - g*g;
X x23 = bp*cp - f*f;
X#ifdef DEBUG
X fprintf( stderr,"x12=%g, y12=%g, z12=%g, x13=%g, y13=%g, x23=%g\n",
X x12, y12, z12, x13, y13, x23 );
X#endif
X vmax = 0.0;
X icase = 0;
X if( fabs(x12) > vmax ) {
X icase = 1;
X vmax = fabs(x12);
X };
X if( fabs(y12) > vmax ) {
X icase = 2;
X vmax = fabs(y12);
X };
X if( fabs(z12) > vmax ) {
X icase = 3;
X vmax = fabs(z12);
X };
X if( fabs(x13) > vmax ) {
X icase = 4;
X vmax = fabs(x13);
X };
X if( fabs(y13) > vmax ) {
X icase = 5;
X vmax = fabs(y13);
X };
X if( fabs(x23) > vmax ) {
X icase = 6;
X vmax = fabs(x23);
X };
X#ifdef DEBUG
X fprintf( stderr, "icase = %d, vmax = %g\n", icase, vmax );
X#endif
X switch(icase) {
X case 0: *xaxis = 0.0;
X *yaxis = 0.0;
X *zaxis = 0.0;
X break;
X case 1: *xaxis = 1.0; /* x12 */
X *yaxis = y12 / x12;
X *zaxis = z12 / x12;
X break;
X case 2: *xaxis = x12 / y12; /* y12 */
X *yaxis = 1.0;
X *zaxis = z12 / y12;
X break;
X case 3: *xaxis = x12 / z12; /* z12 */
X *yaxis = y12 / z12;
X *zaxis = 1.0;
X break;
X case 4: *xaxis = 1.0; /* x13 */
X *yaxis = (-y13) / x13;
X *zaxis = (-y12) / x13; /* z13 = y12 */
X break;
X case 5: *xaxis = (-x13) / y13; /* y13 */
X *yaxis = 1.0;
X *zaxis = y12 / y13; /* z13 = y12 */
X break;
X case 6: *xaxis = 1.0; /* x23 */
X *yaxis = (-x13) / x23; /* y23 = -x13 */
X *zaxis = x12/x23; /* z23 = -x12 */
X break;
X default:
X break;
X };
X#ifdef DEBUG
X fprintf( stderr, "xaxis=%g, yaxis=%g, zaxis=%g\n", *xaxis, *yaxis, *zaxis );
X#endif
X
X if( solve_for_max ) {
X#ifdef DEBUG
X fprintf( stderr, "Have to grab a vector which has a dot prod with the _MAX_ of zero.\n" );
X#endif
X if(*zaxis != 0.0 ){
X newx = 1.0;
X newy = 0.0;
X newz = (-(*xaxis)) / *zaxis;
X }
X else if (*yaxis != 0.0) {
X newx = 1.0;
X newz = 0.0;
X newy = (-(*xaxis)) / *yaxis;
X }
X else if(*xaxis != 0.0) {
X newx = 0.0;
X newy = 1.0;
X newz = 1.0;
X }
X else {
X newx = 0.0;
X newy = 0.0;
X newz = 0.0;
X };
X *xaxis = newx;
X *yaxis = newy;
X *zaxis = newz;
X };
X#ifdef DEBUG
X fprintf( stderr, "INERTIA- returning ( %g, %g, %g )\n", *xaxis, *yaxis, *zaxis );
X#endif
X }
X
X/* analytic method of finding the roots of a cubic equation (Cardano's
X Formula). The general method is multi-step :
X
X 3 2
X given equation x + a x + a x + a = 0
X 1 2 3
X
X we change the equation to the form (by substituting x = (y - a1/3)
X
X 3
X y + p Y + q = 0
X
X where
X
X p = a2 - (a1**2)/3, q = (2*a1**3)/27 - a1*a2/3 + a3
X
X and then, subtituting ( u + v ) for y, and reducing, we get
X
X u1 = (-q/2 + sqrt((q/2)**2 + (p/3)**3))**(1/3)
X v1 = (-q/2 - sqrt((q/2)**2 + (p/3)**3))**(1/3)
X
X and if w1 and w2 are the complex cube roots of unity,
X
X y1 = u1 + v1
X y2 = u1 * w1 + v1 * w2
X y3 = u1 * w2 + v1 * w1
X
X and so
X
X x1 = y1 - a1/3
X x2 = y2 - a1/3
X x3 = y3 - a1/3
X
X The special cases are where one or more of the roots are zero. These
X are handled first (a3==0 => one zero root, a3==0 && a2==0 => two zero roots)
X
X given the above equations, we only need to calculate one cube root of a complex
X number such that -PI/3 <= atan2(ri,rr) <= PI/3.
X
X Copyright 1994 Christopher A. Huson
X */
X
X#define COMPLEXMUL(resr, resi, ar, ai, br, bi) resr = (ar*br)-(ai*bi); resi = (ar*bi)+(ai*br);
X#define COMPLEXADD(resr, resi, ar, ai, br, bi) resr = ar + br; resi = ai + bi;
X
XVOID
Xcubic_roots(a1, a2, a3, r0, r1r, r1i, r2r, r2i )
Xdouble a1, a2, a3, *r0, *r1r, *r1i, *r2r, *r2i;
X{ double p, q, t0r, t0i, t1r, t1i, discr, discri, discrr;
X double u13r, u13i, v13r, v13i, u1r, u1i;
X double v1r, v1i, y1r, y1i, w1r, w1i, w2r, w2i, y2r, y2i, y3r, y3i;
X double x0r, x0i, x1r, x1i, x2r, x2i;
X double val11, val12, val21, val22, offset;
X double sr, si, tr, ti, scale, t;
X double sumr, sumi;
X int is_imag;
X#ifdef DEBUG
X fprintf( stderr, "a1=%g, a2=%g, a3=%g\n", a1, a2, a3 );
X#endif
X
X if( fabs(a3) < EPSILON ) {
X#ifdef DEBUG
X fprintf( stderr, "Special case; 0 is one of the roots\n" );
X fprintf( stderr, "Now, solve the quadratic involving a1 and a2\n" );
X#endif
X *r0 = 0.0;
X if( fabs(a2) < EPSILON ) {
X#ifdef DEBUG
X fprintf( stderr, "Quadratic also has a root of zero.\n" );
X fprintf( stderr, "remaining root is -a1.\n" );
X *r1r = *r1i = 0.0;
X *r2r = -a1;
X *r2i = 0.0;
X#endif
X }
X else {
X#ifdef DEBUG
X fprintf( stderr, "Have to call quadratic solver\n" );
X#endif
X quadratic_roots(1.0,a1,a2,r1r,r1i,r2r,r2i,EPSILON);
X };
X#ifdef DEBUG
X fprintf( stderr, "Solutions are ( %g, %g ), ( %g, %g ), ( %g, %g )\n",
X *r0, 0.0, *r1r, *r1i, *r2r, *r2i );
X check_cubic_solution( a1, a2, a3, *r0, 0.0, &sumr, &sumi );
X check_cubic_solution( a1, a2, a3, *r1r, *r1i, &sumr, &sumi );
X check_cubic_solution( a1, a2, a3, *r2r, *r2i, &sumr, &sumi );
X#endif
X return;
X };
X
X p = a2 - (a1*a1)/3;
X q = (2*a1*a1*a1)/27 - a1*a2/3 + a3;
X#ifdef DEBUG
X fprintf( stderr, "p=%g, q=%g\n", p, q );
X#endif
X
X discr = (q*q/4.0) + (p*p*p/27.0);
X#ifdef DEBUG
X fprintf( stderr, "We have discr = %g\n", discr );
X#endif
X if (discr > 0) {
X discri = 0.0;
X discrr = sqrt(discr);
X }
X else {
X discrr = 0.0;
X discri = sqrt(fabs(discr));
X }
X#ifdef DEBUG
X fprintf( stderr, "And its sqrt is ( %g, %g )\n", discrr, discri );
X#endif
X
X u13r = (-q/2) + discrr;
X u13i = discri;
X v13r = (-q/2) - discrr;
X v13i = -discri;
X
X#ifdef DEBUG
X fprintf( stderr, "We have u1**3=( %g, %g ), v1**3=( %g, %g )\n", u13r, u13i, v13r, v13i );
X#endif
X
X cube_root(u13r, u13i, &u1r, &u1i);
X cube_root(v13r, v13i, &v1r, &v1i);
X
X#ifdef DEBUG
X fprintf( stderr, "u1=( %g, %g ), v1 = ( %g, %g )\n", u1r, u1i, v1r, v1i );
X#endif
X
X COMPLEXADD( y1r, y1i, u1r, u1i, v1r, v1i );
X
X#ifdef DEBUG
X fprintf( stderr, "y1 = ( %g, %g )\n", y1r, y1i );
X#endif
X
X w1r = -0.5;
X w1i = sqrt(3.0)/2.0;
X w2r = w1r;
X w2i = -w1i;
X
X COMPLEXMUL(t0r,t0i, u1r, u1i, w1r, w1i );
X COMPLEXMUL(t1r,t1i, v1r, v1i, w2r, w2i );
X COMPLEXADD(y2r, y2i, t0r,t0i, t1r,t1i );
X
X COMPLEXMUL(t0r,t0i, u1r, u1i, w2r, w2i );
X COMPLEXMUL(t1r,t1i, v1r, v1i, w1r, w1i );
X COMPLEXADD(y3r,y3i, t0r,t0i, t1r,t1i );
X
X#ifdef DEBUG
X fprintf( stderr, "y2=( %g, %g ), y3=( %g, %g )\n", y2r, y2i, y3r,y3i );
X#endif
X
X COMPLEXADD(x0r, x0i, y1r, y1i, (-a1/3), 0 );
X COMPLEXADD(x1r, x1i, y2r, y2i, (-a1/3), 0 );
X COMPLEXADD(x2r, x2i, y3r,y3i, (-a1/3), 0 );
X
X#ifdef DEBUG
X fprintf( stderr, "Solutions are ( %g, %g ), ( %g, %g ), ( %g, %g )\n",
X x0r, x0i, x1r, x1i, x2r, x2i );
X check_cubic_solution( a1, a2, a3, x0r, x0i, &sumr, &sumi );
X check_cubic_solution( a1, a2, a3, x1r, x1i, &sumr, &sumi );
X check_cubic_solution( a1, a2, a3, x2r, x2i, &sumr, &sumi );
X#endif
X
X /* put the root with the smallest absolute value imaginary (which should be
X the real root) as the first root. */
X
X if(fabs(x0i)>fabs(x1i)){
X t = x0i;
X x0i = x1i;
X x1i = t;
X t = x0r;
X x0r = x1r;
X x1r = t;
X };
X
X if(fabs(x0i)>fabs(x2i)) {
X t = x0i;
X x0i = x2i;
X x2i = t;
X t = x0r;
X x0r = x2r;
X x2r = t;
X };
X
X *r0 = x0r;
X *r1r = x1r;
X *r1i = x1i;
X *r2r = x2r;
X *r2i = x2i;
X#ifdef DEBUG
X fprintf( stderr, "Returning from cubic_roots\n" );
X#endif
X }
X
X/*
X| find cube root of a complex number. If the number is negative
X| real, the root is negative real, otherwise, it will fall between
X| (-pi/3 : pi/3)
X*/
X
X#define PI 3.1415926535897932384626433
X#define EPSILONSQ 1.0e-12
X
X#ifdef __STDC__
XVOID
Xcube_root( double rr, double ri, double *r0r, double *r0i )
X#else /* __STDC__ */
Xcube_root( rr, ri, r0r, r0i )
Xdouble rr, ri, *r0r, *r0i;
X#endif /* __STDC__ */
X{
X double theta, phi, radius, sinrot, cosrot, tempr, tempi;
X double rr2, ri2;
X
X#ifdef DEBUG
X fprintf( stderr, "Starting with ( %g, %g )\n", rr, ri );
X#endif
X
X if((rr == 0.0)&&(ri == 0.0)) {
X#ifdef DEBUG
X fprintf( stderr, "Is zero\n" );
X#endif
X *r0r = *r0i = 0.0;
X return;
X };
X
X if ( rr2 != 0.0 ) {
X rr2 = rr * rr;
X ri2 = ri * ri;
X if((ri2 / rr2) < EPSILONSQ) {
X ri = 0.0;
X }
X }
X if((ri == 0.0)&&(rr < 0.0)) {
X radius = -rr;
X radius = log(radius) / 3.0;
X radius = exp(radius);
X *r0r = -radius;
X *r0i = 0.0;
X return;
X }
X
X radius = sqrt(rr*rr + ri*ri);
X theta = atan2( ri, rr );
X
X#ifdef DEBUG
X fprintf( stderr, "Equivalent to ( %g, %g ) in polar coords\n", radius, theta );
X#endif
X
X radius = log(radius) / 3.0;
X radius = exp(radius);
X
X phi = theta / 3.0;
X
X#ifdef DEBUG
X fprintf( stderr, "Base cube root is ( %g, %g ) in polar coords\n", radius, phi );
X#endif
X
X if((ri==0.0)&&(rr < 0.0)) {
X phi += 2.0*PI/3.0;
X#ifdef DEBUG
X fprintf( stderr, "Rotate phi counterclockwise to %g\n", phi );
X#endif
X };
X *r0r = radius * cos( phi );
X *r0i = radius * sin( phi );
X
X#ifdef DEBUG
X fprintf( stderr, "Which is ( %g, %g ) in rectangular coords\n", *r0r, *r0i );
X#endif
X
X }
X#ifdef DEBUG
X/* procedure to use under DEBUG which checks the roots of the cubic */
X/* Only use the following with simple variables */
X
XVOID
Xcheck_cubic_solution( a1, a2, a3, rr, ri, sumr, sumi )
Xdouble a1, a2, a3; /* the eq is x**3 + a1*x**2 + a2+x + a3 = 0 */
Xdouble rr, ri; /* the real and imaginary parts of the proposed solution */
Xdouble *sumr, *sumi;
X{
X double rr_2, ri_2, rr_3, ri_3; /* the squared and cubed versions */
X
X COMPLEXMUL(rr_2, ri_2, rr, ri, rr, ri );
X COMPLEXMUL(rr_3, ri_3, rr, ri, rr_2, ri_2 );
X *sumr = rr_3 + a1 * rr_2 + a2 * rr + a3;
X *sumi = ri_3 + a1 * ri_2 + a2 * ri;
X fprintf( stderr, "For ( %g, %g ), x**3 + %g * x**2 + %g * x + %g = ( %g, %g )\n",
X rr, ri, a1, a2, a3, *sumr, *sumi );
X }
X#endif
X/* simple; given a quadratic equation of the form
X
X 2
X a * x + b * x + c = 0
X
X the last parameter is the error factor. If it is not 0,
X then if the ratio of the real to imaginary part is < error,
X then just set the imaginary part to 0.
X*/
X
XVOID
Xquadratic_roots(a,b,c,x1r,x1i,x2r,x2i,error)
Xdouble a,b,c;
Xdouble *x1r,*x1i,*x2r,*x2i;
Xdouble error;
X{
X double disc,p1,p2;
X#ifdef DEBUG
X fprintf( stderr, "QUADRATIC- a=%g, b=%g, c=%g\n", a, b, c );
X#endif
X if ( a == 0e0) {
X#ifdef DEBUG
X fprintf( stderr, "QUADRATIC- a == 0\n" );
X#endif
X if(b == 0e0) {
X#ifdef DEBUG
X fprintf( stderr, "QUADRATIC- b == 0\n" );
X#endif
X /* error : infinite number of solutions (or none) */
X pm_error( "QUADRATIC : cannot find solution." );
X };
X /* only one solution. We just set them both to the same */
X *x1r = *x2r = (-c) / b;
X *x1i = *x2i = 0e0;
X#ifdef DEBUG
X fprintf( stderr, "QUADRATIC- returning x1r = x2r = %g\n", *x1r );
X#endif
X return;
X };
X disc = b*b - 4 * a * c;
X#ifdef DEBUG
X fprintf( stderr, "QUADRATIC- disc = %g\n", disc );
X#endif
X if (disc >= 0e0) {
X p2 = sqrt(disc) / (2e0 * a);
X p1 = 0e0 - (b / ( 2e0 * a));
X#ifdef DEBUG
X fprintf( stderr, "QUADRATIC- p1 = %g, p2 = %g\n", p1, p2 );
X#endif
X *x1r = p1 + p2;
X *x2r = p1 - p2;
X *x1i = *x2i = 0e0;
X }
X else { /* conjugate complex roots */
X p2 = sqrt(fabs(disc)) / (2e0 * a);
X p1 = 0e0 - (b / (2e0 * a));
X#ifdef DEBUG
X fprintf( stderr, "QUADRATIC- Complex roots, p1 = %g, p2 = %g\n", p1, p2 );
X#endif
X *x1r = *x2r = p1;
X if(fabs( p2 ) < error) {
X /* with Newton-Raphson, we may be slightly off from correct
X root. If so, an imaginary part may creep in. We use the
X error parameter to tell if we expect this to happen */
X *x1i = *x2i = 0e0;
X }
X else {
X *x1i = p2;
X *x2i = (-(p2));
X }
X };
X }
X
XVOID
X#ifdef __STDC__
Xset_color_weights( double rval, double gval, double bval )
X#else /* __STDC__ */
Xset_color_weights( rval, gval, bval )
Xdouble rval, gval, bval;
X#endif /* __STDC__ */
X{
X red_weight = rval;
X green_weight = gval;
X blue_weight = bval;
X }
END_OF_FILE
if test 64461 -ne `wc -c <'inert_cut.c'`; then
echo shar: \"'inert_cut.c'\" unpacked with wrong size!
fi
# end of 'inert_cut.c'
fi
echo shar: End of archive 2 \(of 2\).
cp /dev/null ark2isdone