[gs-cvs] rev 7034 - in trunk/gs: . imdi

giles at ghostscript.com giles at ghostscript.com
Mon Sep 11 13:26:04 PDT 2006


Author: giles
Date: 2006-09-11 13:26:01 -0700 (Mon, 11 Sep 2006)
New Revision: 7034

Added:
   trunk/gs/imdi/
   trunk/gs/imdi/Jamfile
   trunk/gs/imdi/LICENSE
   trunk/gs/imdi/README
   trunk/gs/imdi/arch.h
   trunk/gs/imdi/cctiff.c
   trunk/gs/imdi/cgen.c
   trunk/gs/imdi/config.h
   trunk/gs/imdi/copyright.h
   trunk/gs/imdi/imdi.c
   trunk/gs/imdi/imdi.h
   trunk/gs/imdi/imdi_gen.c
   trunk/gs/imdi/imdi_gen.h
   trunk/gs/imdi/imdi_imp.h
   trunk/gs/imdi/imdi_k.c
   trunk/gs/imdi/imdi_k.h
   trunk/gs/imdi/imdi_tab.c
   trunk/gs/imdi/imdi_tab.h
Log:
Check in working files for the GPL imdi (integer multi-dimensional 
interpolation) library for color mapping. This is needed by the imdi
device. Port from the ghostpcl tree.


Added: trunk/gs/imdi/Jamfile
===================================================================
--- trunk/gs/imdi/Jamfile	2006-09-11 07:02:18 UTC (rev 7033)
+++ trunk/gs/imdi/Jamfile	2006-09-11 20:26:01 UTC (rev 7034)
@@ -0,0 +1,4 @@
+LINKLIBS += -lm ;
+Main imdi_gen : imdi_gen.c cgen.c ;
+GenFile imdi_k.h : imdi_gen ;
+# Library libimdi.lib : imdi.c imdi_tab.c ;

Added: trunk/gs/imdi/LICENSE
===================================================================
--- trunk/gs/imdi/LICENSE	2006-09-11 07:02:18 UTC (rev 7033)
+++ trunk/gs/imdi/LICENSE	2006-09-11 20:26:01 UTC (rev 7034)
@@ -0,0 +1,282 @@
+                   GNU GENERAL PUBLIC LICENSE
+                       Version 2, June 1991
+
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.
+                       59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+                            Preamble
+
+  The licenses for most software are designed to take away your
+freedom to share and change it.  By contrast, the GNU General Public
+License is intended to guarantee your freedom to share and change free
+software--to make sure the software is free for all its users.  This
+General Public License applies to most of the Free Software
+Foundation's software and to any other program whose authors commit to
+using it.  (Some other Free Software Foundation software is covered by
+the GNU Library General Public License instead.)  You can apply it to
+your programs, too.
+
+  When we speak of free software, we are referring to freedom, not
+price.  Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+this service if you wish), that you receive source code or can get it
+if you want it, that you can change the software or use pieces of it
+in new free programs; and that you know you can do these things.
+
+  To protect your rights, we need to make restrictions that forbid
+anyone to deny you these rights or to ask you to surrender the rights.
+These restrictions translate to certain responsibilities for you if you
+distribute copies of the software, or if you modify it.
+
+  For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must give the recipients all the rights that
+you have.  You must make sure that they, too, receive or can get the
+source code.  And you must show them these terms so they know their
+rights.
+
+  We protect your rights with two steps: (1) copyright the software, and
+(2) offer you this license which gives you legal permission to copy,
+distribute and/or modify the software.
+
+  Also, for each author's protection and ours, we want to make certain
+that everyone understands that there is no warranty for this free
+software.  If the software is modified by someone else and passed on, we
+want its recipients to know that what they have is not the original, so
+that any problems introduced by others will not reflect on the original
+authors' reputations.
+
+  Finally, any free program is threatened constantly by software
+patents.  We wish to avoid the danger that redistributors of a free
+program will individually obtain patent licenses, in effect making the
+program proprietary.  To prevent this, we have made it clear that any
+patent must be licensed for everyone's free use or not licensed at all.
+
+  The precise terms and conditions for copying, distribution and
+modification follow.
+
+                    GNU GENERAL PUBLIC LICENSE
+   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+  0. This License applies to any program or other work which contains
+a notice placed by the copyright holder saying it may be distributed
+under the terms of this General Public License.  The "Program", below,
+refers to any such program or work, and a "work based on the Program"
+means either the Program or any derivative work under copyright law:
+that is to say, a work containing the Program or a portion of it,
+either verbatim or with modifications and/or translated into another
+language.  (Hereinafter, translation is included without limitation in
+the term "modification".)  Each licensee is addressed as "you".
+
+Activities other than copying, distribution and modification are not
+covered by this License; they are outside its scope.  The act of
+running the Program is not restricted, and the output from the Program
+is covered only if its contents constitute a work based on the
+Program (independent of having been made by running the Program).
+Whether that is true depends on what the Program does.
+
+  1. You may copy and distribute verbatim copies of the Program's
+source code as you receive it, in any medium, provided that you
+conspicuously and appropriately publish on each copy an appropriate
+copyright notice and disclaimer of warranty; keep intact all the
+notices that refer to this License and to the absence of any warranty;
+and give any other recipients of the Program a copy of this License
+along with the Program.
+
+You may charge a fee for the physical act of transferring a copy, and
+you may at your option offer warranty protection in exchange for a fee.
+
+  2. You may modify your copy or copies of the Program or any portion
+of it, thus forming a work based on the Program, and copy and
+distribute such modifications or work under the terms of Section 1
+above, provided that you also meet all of these conditions:
+
+    a) You must cause the modified files to carry prominent notices
+    stating that you changed the files and the date of any change.
+
+    b) You must cause any work that you distribute or publish, that in
+    whole or in part contains or is derived from the Program or any
+    part thereof, to be licensed as a whole at no charge to all third
+    parties under the terms of this License.
+
+    c) If the modified program normally reads commands interactively
+    when run, you must cause it, when started running for such
+    interactive use in the most ordinary way, to print or display an
+    announcement including an appropriate copyright notice and a
+    notice that there is no warranty (or else, saying that you provide
+    a warranty) and that users may redistribute the program under
+    these conditions, and telling the user how to view a copy of this
+    License.  (Exception: if the Program itself is interactive but
+    does not normally print such an announcement, your work based on
+    the Program is not required to print an announcement.)
+
+These requirements apply to the modified work as a whole.  If
+identifiable sections of that work are not derived from the Program,
+and can be reasonably considered independent and separate works in
+themselves, then this License, and its terms, do not apply to those
+sections when you distribute them as separate works.  But when you
+distribute the same sections as part of a whole which is a work based
+on the Program, the distribution of the whole must be on the terms of
+this License, whose permissions for other licensees extend to the
+entire whole, and thus to each and every part regardless of who wrote it.
+
+Thus, it is not the intent of this section to claim rights or contest
+your rights to work written entirely by you; rather, the intent is to
+exercise the right to control the distribution of derivative or
+collective works based on the Program.
+
+In addition, mere aggregation of another work not based on the Program
+with the Program (or with a work based on the Program) on a volume of
+a storage or distribution medium does not bring the other work under
+the scope of this License.
+
+  3. You may copy and distribute the Program (or a work based on it,
+under Section 2) in object code or executable form under the terms of
+Sections 1 and 2 above provided that you also do one of the following:
+
+    a) Accompany it with the complete corresponding machine-readable
+    source code, which must be distributed under the terms of Sections
+    1 and 2 above on a medium customarily used for software interchange; or,
+
+    b) Accompany it with a written offer, valid for at least three
+    years, to give any third party, for a charge no more than your
+    cost of physically performing source distribution, a complete
+    machine-readable copy of the corresponding source code, to be
+    distributed under the terms of Sections 1 and 2 above on a medium
+    customarily used for software interchange; or,
+
+    c) Accompany it with the information you received as to the offer
+    to distribute corresponding source code.  (This alternative is
+    allowed only for noncommercial distribution and only if you
+    received the program in object code or executable form with such
+    an offer, in accord with Subsection b above.)
+
+The source code for a work means the preferred form of the work for
+making modifications to it.  For an executable work, complete source
+code means all the source code for all modules it contains, plus any
+associated interface definition files, plus the scripts used to
+control compilation and installation of the executable.  However, as a
+special exception, the source code distributed need not include
+anything that is normally distributed (in either source or binary
+form) with the major components (compiler, kernel, and so on) of the
+operating system on which the executable runs, unless that component
+itself accompanies the executable.
+
+If distribution of executable or object code is made by offering
+access to copy from a designated place, then offering equivalent
+access to copy the source code from the same place counts as
+distribution of the source code, even though third parties are not
+compelled to copy the source along with the object code.
+
+  4. You may not copy, modify, sublicense, or distribute the Program
+except as expressly provided under this License.  Any attempt
+otherwise to copy, modify, sublicense or distribute the Program is
+void, and will automatically terminate your rights under this License.
+However, parties who have received copies, or rights, from you under
+this License will not have their licenses terminated so long as such
+parties remain in full compliance.
+
+  5. You are not required to accept this License, since you have not
+signed it.  However, nothing else grants you permission to modify or
+distribute the Program or its derivative works.  These actions are
+prohibited by law if you do not accept this License.  Therefore, by
+modifying or distributing the Program (or any work based on the
+Program), you indicate your acceptance of this License to do so, and
+all its terms and conditions for copying, distributing or modifying
+the Program or works based on it.
+
+  6. Each time you redistribute the Program (or any work based on the
+Program), the recipient automatically receives a license from the
+original licensor to copy, distribute or modify the Program subject to
+these terms and conditions.  You may not impose any further
+restrictions on the recipients' exercise of the rights granted herein.
+You are not responsible for enforcing compliance by third parties to
+this License.
+
+  7. If, as a consequence of a court judgment or allegation of patent
+infringement or for any other reason (not limited to patent issues),
+conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License.  If you cannot
+distribute so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you
+may not distribute the Program at all.  For example, if a patent
+license would not permit royalty-free redistribution of the Program by
+all those who receive copies directly or indirectly through you, then
+the only way you could satisfy both it and this License would be to
+refrain entirely from distribution of the Program.
+
+If any portion of this section is held invalid or unenforceable under
+any particular circumstance, the balance of the section is intended to
+apply and the section as a whole is intended to apply in other
+circumstances.
+
+It is not the purpose of this section to induce you to infringe any
+patents or other property right claims or to contest validity of any
+such claims; this section has the sole purpose of protecting the
+integrity of the free software distribution system, which is
+implemented by public license practices.  Many people have made
+generous contributions to the wide range of software distributed
+through that system in reliance on consistent application of that
+system; it is up to the author/donor to decide if he or she is willing
+to distribute software through any other system and a licensee cannot
+impose that choice.
+
+This section is intended to make thoroughly clear what is believed to
+be a consequence of the rest of this License.
+
+  8. If the distribution and/or use of the Program is restricted in
+certain countries either by patents or by copyrighted interfaces, the
+original copyright holder who places the Program under this License
+may add an explicit geographical distribution limitation excluding
+those countries, so that distribution is permitted only in or among
+countries not thus excluded.  In such case, this License incorporates
+the limitation as if written in the body of this License.
+
+  9. The Free Software Foundation may publish revised and/or new versions
+of the General Public License from time to time.  Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+Each version is given a distinguishing version number.  If the Program
+specifies a version number of this License which applies to it and "any
+later version", you have the option of following the terms and conditions
+either of that version or of any later version published by the Free
+Software Foundation.  If the Program does not specify a version number of
+this License, you may choose any version ever published by the Free Software
+Foundation.
+
+  10. If you wish to incorporate parts of the Program into other free
+programs whose distribution conditions are different, write to the author
+to ask for permission.  For software which is copyrighted by the Free
+Software Foundation, write to the Free Software Foundation; we sometimes
+make exceptions for this.  Our decision will be guided by the two goals
+of preserving the free status of all derivatives of our free software and
+of promoting the sharing and reuse of software generally.
+
+                            NO WARRANTY
+
+  11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+  12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+
+                     END OF TERMS AND CONDITIONS
+
+

Added: trunk/gs/imdi/README
===================================================================
--- trunk/gs/imdi/README	2006-09-11 07:02:18 UTC (rev 7033)
+++ trunk/gs/imdi/README	2006-09-11 20:26:01 UTC (rev 7034)
@@ -0,0 +1,114 @@
+
+This is the development area for IMDI, the
+Interger Multi-Dimensional Interpolation routines.
+
+They provide a flexible and high performance
+system for applying color transforms to typical
+raster pixel data. Because they provide a means of
+applying arbitrary combination dependent mappings
+of multi-channel pixel data, there are many other
+possible uses for these sorts of routines as well,
+including high quality matting/compositing. For instance,
+one could create a smooth, proportional "chroma key"
+type of matt for matting one RGB image onto another
+by creating a 6 channel to 3 dimensional transform,
+that its applied to each pair of pixels from the
+source images and produces one combined output pixel.
+Additional input or output alpha channels are easy
+to add by simply adding more input and/or output
+dimensions. The matting calculatons can be almost
+arbitrarily complex, and the imdi will apply them
+to the pixel data at very high speed.
+
+The system has two parts, one that generates taylored,
+optimised source code for the transformation kernels,
+and the run time code that matches a transform request
+to a compiled kernel, and initialises the appropriate
+run time lookup tables.
+
+The kernel source generator is intended to accomodate
+various optimisations, such as assembly code, vector
+instruction set (ie. MMX, AltiVec etc.) versions, but
+at present only generates the more portable 'C' code
+kernels.
+
+Both 8 bit per component and 16 bit per component
+pixel data is handled, up to 8 input and output
+dimensions (but this limit could be trivially raised).
+
+imdi_gen.exe	is the module that triggers the generation of
+		optimised source code as configured for the color spaces
+		and pixel formats selected. By default creates
+		a single imdi_k.c and imdi_k.h file, but if
+		given the -i flag, creates a separate file
+		for each kernel variant.
+
+cgen.c		C code generator module.
+
+itest.c	regresion test routine.
+		Normally runs speed and accuracy tests for
+		all configured kernel variants.
+		The -q flag makes it run quicker,
+		but makes the benchmarking inacurate,
+		the -s flag will cause it to stop
+		if any routine has unexpectedly low
+		accuracy.
+
+cctiff.c	is the utility that takes an ICC device
+		profile link, and converts a TIFF file
+		from the input colorspace to the output
+		space. Both 8 bit and 16 bit TIFF files
+		are handled, as well as colorspaces up to
+		8 channels in and out.
+        This accepts either a device link ICC profile,
+        or links an input and output devce ICC profile
+        to define the color transform.
+
+
+greytiff.c  is a utility similar to cctiff, that
+        is an example of how to colorimetrically
+        convert an RGB file into a monochrome RGB file.
+
+
+
+Misc. Notes
+-----------
+
+   ITU-T Rec. T.42 specifies the ITULAB encoding in terms of a range
+    and offset for each component, which are related to the minimum and
+    maximum values as follows:
+
+        minimum = - (range x offset) / 2^n - 1
+        maximum = minimum + range
+
+    The Decode field default values depend on the color space. For the
+    ITULAB color space encoding, the default values correspond to the
+    base range and offset, as specified in ITU-T Rec. T.42 [T.42]. The
+    following table gives the base range and offset values for
+    BitsPerSample=8 and 12, and the corresponding default minimum and
+    maximum default values for the Decode field, calculated using the
+    equations above when PhotometricInterpetation=10.
+
+                       +-----------------------------------------------+
+                       | ITU-T Rec. T.42  |           Decode           |
+ +---------+-----------|   base values    |       default values       |
+ | BitsPer + Component +------------------+----------------------------+
+ | -Sample |           |  Range | Offset  |      Min     |     Max     |
+ +---------+-----------+--------+---------+--------------+-------------+
+ |    8    |    L*     |   100  |    0    |       0      |     100     |
+ |         +-----------+--------+---------+--------------+-------------+
+ |         |    a*     |   170  |   128   |  -21760/255  |  21590/255  |
+ |         +-----------+--------+---------+--------------+-------------+
+ |         |    b*     |   200  |    96   |  -19200/255  |  31800/255  |
+ +---------+-----------+--------+---------+--------------+-------------+
+ |   12    |    L*     |   100  |    0    |       0      |     100     |
+ |         +-----------+--------+---------+--------------+-------------+
+ |         |    a*     |   170  |  2048   | -348160/4095 | 347990/4095 |
+ |         +-----------+--------+---------+--------------+-------------+
+ |         |    b*     |   200  |  1536   | -307200/4095 | 511800/4095 |
+ +---------+-----------+--------+---------+--------------+-------------+
+
+   For example, when PhotometricInterpretation=10 and BitsPerSample=8,
+   the default value for Decode is (0, 100, -21760/255, 21590/255,
+   -19200/255, 31800/255).
+

Added: trunk/gs/imdi/arch.h
===================================================================
--- trunk/gs/imdi/arch.h	2006-09-11 07:02:18 UTC (rev 7033)
+++ trunk/gs/imdi/arch.h	2006-09-11 20:26:01 UTC (rev 7034)
@@ -0,0 +1,58 @@
+#ifndef ARCH_H
+#define ARCH_H
+
+/* Integer Multi-Dimensional Interpolation */
+/*
+ * Copyright 2000 Graeme W. Gill
+ *
+ * This material is licenced under the GNU GENERAL PUBLIC LICENCE :-
+ * see the Licence.txt file for licencing details.
+ */
+
+#define STR_DEF(def)  #def
+
+#ifdef ALLOW64
+
+/* Detect machine/compiler specifics here */
+#if defined(NT)
+#define longlong __int64
+#else	/* !NT, assume standard */
+#define longlong long long
+#endif	/* !NT */
+#define str_longlong STR_DEF(longlong)
+
+#endif /* ALLOW64 */
+
+
+
+/* Machine/Language architectural specifications */
+typedef struct {
+	int bits;		/* Bits in this data type */
+	char *name;		/* Name used to specify this type */
+	int align;		/* Non-zero if this type should be accessed aligned */
+} dtypes;
+
+#define MXDTYPES 6
+
+typedef struct {
+	int bigend;		/* Non-zero if this is a bigendian architecture */
+	int uwa;		/* Use wide memory access */
+
+	int    pbits;	/* Number of bits in a pointer */
+
+	int    nords;	/* Number of ord types */
+	dtypes ords[MXDTYPES];		/* Ordinal types, in size order */
+	int    natord;	/* Index of natural machine ordinal */
+
+	int    nints;	/* Number of int types */
+	dtypes ints[MXDTYPES];		/* Integer types, in size order */
+	int    natint;	/* Index of natural machine integer */
+
+	/* Optimisation settings */
+	int    shfm;	/* Non-zero to use shifts for masking */
+	int    oscale;	/* Maximum power of 2 scaled indexing mode, 0 for none. */
+	int    smmul;	/* Has fast small multiply for index scaling */
+
+} mach_arch;
+
+#endif /* ARCH_H */

Added: trunk/gs/imdi/cctiff.c
===================================================================
--- trunk/gs/imdi/cctiff.c	2006-09-11 07:02:18 UTC (rev 7033)
+++ trunk/gs/imdi/cctiff.c	2006-09-11 20:26:01 UTC (rev 7034)
@@ -0,0 +1,1190 @@
+
+/* 
+ * Color Correct a TIFF file, using an ICC Device link profile.
+ *
+ * Author:  Graeme W. Gill
+ * Date:    00/3/8
+ * Version: 1.30
+ *
+ * Copyright 2000 - 2004 Graeme W. Gill
+ * All rights reserved.
+ *
+ * This material is licenced under the GNU GENERAL PUBLIC LICENCE :-
+ * see the Licence.txt file for licencing details.
+ */
+
+/*
+ * Thanks to Neil Okamoto for the 16 bit TIFF mods.
+ */
+
+/* TTBD:
+ */
+
+/*
+	This program is a framework that exercises the
+	IMDI code, as well as a demonstration of simple
+    profile linking.  It can also do the conversion using the
+    floating point code in ICCLIB as a reference.
+
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdarg.h>
+#include <fcntl.h>
+#include <string.h>
+#include <math.h>
+#include "copyright.h"
+#include "config.h"
+#include "tiffio.h"
+#include "icc.h"
+#include "imdi.h"
+
+#undef TREAT_CMY_AS_RGB
+
+void error(char *fmt, ...), warning(char *fmt, ...);
+
+void usage(void) {
+	fprintf(stderr,"Color Correct a TIFF file using an ICC device link profile, V%s\n",ARGYLL_VERSION_STR);
+	fprintf(stderr,"Author: Graeme W. Gill, licensed under the GPL\n");
+	fprintf(stderr,"usage: cctiff [-options] devlinkprofile.icm infile.tif outfile.tif\n");
+	fprintf(stderr,"usage: cctiff [-options] -l inprofile.icm outprofile.icm infile.tif outfile.tif\n");
+	fprintf(stderr," -v            Verbose\n");
+	fprintf(stderr," -c            Combine linearisation curves into one transform\n");
+	fprintf(stderr," -p            Use slow precise correction\n");
+	fprintf(stderr," -k            Check fast result against precise, and report\n");
+	fprintf(stderr," -l            Link input and output profiles\n"); 
+	fprintf(stderr,"  -i in_intent  p = perceptual, r = relative colorimetric,\n");
+	fprintf(stderr,"                s = saturation, a = absolute colorimetric\n");
+	fprintf(stderr,"  -o out_intent p = perceptual, r = relative colorimetric,\n");
+	fprintf(stderr,"                s = saturation, a = absolute colorimetric\n");
+	exit(1);
+}
+
+/* Convert an ICC colorspace to the corresponding possible TIFF Photometric tags. */
+/* Return the number of matching tags, and 0 if there is no corresponding tag. */
+int
+ColorSpaceSignature2TiffPhotometric(
+uint16 tags[10],					/* Pointer to return array, up to 10 */
+icColorSpaceSignature cspace	/* Input ICC colorspace */
+) {
+	switch(cspace) {
+		case icSigGrayData:
+			tags[0] = PHOTOMETRIC_MINISBLACK;
+			return 1;
+		case icSigRgbData:
+#ifdef TREAT_CMY_AS_RGB
+		case icSigCmyData:
+#endif
+			tags[0] = PHOTOMETRIC_RGB;
+			return 1;
+#ifndef TREAT_CMY_AS_RGB
+		case icSigCmyData:
+#endif
+		case icSigCmykData:
+			tags[0] = PHOTOMETRIC_SEPARATED;
+			return 1;
+		case icSigYCbCrData:
+			tags[0] = PHOTOMETRIC_YCBCR;
+			return 1;
+		case icSigLabData:
+			tags[0] = PHOTOMETRIC_CIELAB;
+#ifdef PHOTOMETRIC_ICCLAB
+			tags[1] = PHOTOMETRIC_ICCLAB;
+			tags[2] = PHOTOMETRIC_ITULAB;
+#endif
+			return 3;
+
+		case icSigXYZData:
+		case icSigLuvData:
+		case icSigYxyData:
+		case icSigHsvData:
+		case icSigHlsData:
+			return 0;
+
+		case icSig2colorData:
+			tags[0] = PHOTOMETRIC_SEPARATED;
+			tags[1] = 2;			/* Cheat */
+			return 1;
+
+		case icSig3colorData:
+			tags[0] = PHOTOMETRIC_SEPARATED;
+			tags[1] = 3;			/* Cheat */
+			return 1;
+
+		case icSig4colorData:
+			tags[0] = PHOTOMETRIC_SEPARATED;
+			tags[1] = 4;			/* Cheat */
+			return 1;
+
+		case icSig5colorData:
+		case icSigMch5Data:
+			tags[0] = PHOTOMETRIC_SEPARATED;
+			tags[1] = 5;			/* Cheat */
+			return 1;
+
+		case icSig6colorData:
+		case icSigMch6Data:
+			tags[0] = PHOTOMETRIC_SEPARATED;
+			tags[1] = 6;			/* Cheat */
+			return 1;
+
+		case icSig7colorData:
+		case icSigMch7Data:
+			tags[0] = PHOTOMETRIC_SEPARATED;
+			tags[1] = 7;			/* Cheat */
+			return 1;
+
+		case icSig8colorData:
+		case icSigMch8Data:
+			tags[0] = PHOTOMETRIC_SEPARATED;
+			tags[1] = 8;			/* Cheat */
+			return 1;
+
+		case icSig9colorData:
+			tags[0] = PHOTOMETRIC_SEPARATED;
+			tags[1] = 9;			/* Cheat */
+			return 1;
+
+		case icSig10colorData:
+			tags[0] = PHOTOMETRIC_SEPARATED;
+			tags[1] = 10;			/* Cheat */
+			return 1;
+
+		case icSig11colorData:
+			tags[0] = PHOTOMETRIC_SEPARATED;
+			tags[1] = 11;			/* Cheat */
+			return 1;
+
+		case icSig12colorData:
+			tags[0] = PHOTOMETRIC_SEPARATED;
+			tags[1] = 12;			/* Cheat */
+			return 1;
+
+		case icSig13colorData:
+			tags[0] = PHOTOMETRIC_SEPARATED;
+			tags[1] = 13;			/* Cheat */
+			return 1;
+
+		case icSig14colorData:
+			tags[0] = PHOTOMETRIC_SEPARATED;
+			tags[1] = 14;			/* Cheat */
+			return 1;
+
+		case icSig15colorData:
+			tags[0] = PHOTOMETRIC_SEPARATED;
+			tags[1] = 15;			/* Cheat */
+			return 1;
+
+		default:
+			return 0;
+	}
+	return 0;
+}
+
+
+/* Compute the length of a double nul terminated string, including */
+/* the nuls. */
+static int zzstrlen(char *s) {
+	int i;
+	for (i = 0;; i++) {
+		if (s[i] == '\000' && s[i+1] == '\000')
+			return i+2;
+	}
+	return 0;
+}
+
+/* Convert an ICC colorspace to the corresponding TIFF Inkset tag */
+/* return 0xffff if not possible or applicable. */
+
+int
+ColorSpaceSignature2TiffInkset(
+icColorSpaceSignature cspace,
+int *len,						/* Return length of ASCII inknames */
+char **inknames					/* Return ASCII inknames if non NULL */
+) {
+	switch(cspace) {
+		case icSigCmyData:
+			return 0xffff;	// ~~9999
+			if (inknames != NULL) {
+				*inknames = "cyan\000magenta\000yellow\000\000";
+				*len = zzstrlen(*inknames);
+			}
+			return 0;				/* Not CMYK */
+		case icSigCmykData:
+			if (inknames != NULL) {
+				*inknames = NULL;	/* No inknames */
+				*len = 0;
+			}
+			return INKSET_CMYK;
+
+		case icSigGrayData:
+		case icSigRgbData:
+		case icSigYCbCrData:
+		case icSigLabData:
+		case icSigXYZData:
+		case icSigLuvData:
+		case icSigYxyData:
+		case icSigHsvData:
+		case icSigHlsData:
+		case icSig2colorData:
+		case icSig3colorData:
+		case icSig4colorData:
+		case icSig5colorData:
+		case icSigMch5Data:
+			return 0xffff;
+
+		case icSig6colorData:
+		case icSigMch6Data:
+				/* This is a cheat and a hack. Should really make sure that icclink */
+				/* transfers the right information from the destination */
+				/* profile, and then copies it to the device profile, */
+				/* allowing cctiff to read it. */
+			if (inknames != NULL) {
+				*inknames = "cyan\000magenta\000yellow\000black\000orange\000green\000\000";
+				*len = zzstrlen(*inknames);
+			}
+			return 0;				/* Not CMYK */
+
+		case icSig7colorData:
+		case icSigMch7Data:
+			return 0xffff;
+
+		case icSig8colorData:
+		case icSigMch8Data:
+				/* This is a cheat and a hack. Should really make sure that icclink */
+				/* transfers the right information from the destination */
+				/* profile, and then copies it to the device profile, */
+				/* allowing cctiff to read it. */
+			if (inknames != NULL) {
+				*inknames = "cyan\000magenta\000yellow\000black\000orange\000green\000lightcyan\000lightmagenta\000\000";
+				*len = zzstrlen(*inknames);
+			}
+			return 0;				/* Not CMYK */
+		case icSig9colorData:
+		case icSig10colorData:
+		case icSig11colorData:
+		case icSig12colorData:
+		case icSig13colorData:
+		case icSig14colorData:
+		case icSig15colorData:
+		default:
+			return 0xffff;
+	}
+	return 0xffff;
+}
+
+char *
+Photometric2str(
+int pmtc
+) {
+	static char buf[80];
+	switch (pmtc) {
+		case PHOTOMETRIC_MINISWHITE:
+			return "Subtractive Gray";
+		case PHOTOMETRIC_MINISBLACK:
+			return "Additive Gray";
+		case PHOTOMETRIC_RGB:
+			return "RGB";
+		case PHOTOMETRIC_PALETTE:
+			return "Indexed";
+		case PHOTOMETRIC_MASK:
+			return "Transparency Mask";
+		case PHOTOMETRIC_SEPARATED:
+			return "Separated";
+		case PHOTOMETRIC_YCBCR:
+			return "YCbCr";
+		case PHOTOMETRIC_CIELAB:
+			return "CIELab";
+#ifdef PHOTOMETRIC_ICCLAB
+		case PHOTOMETRIC_ICCLAB:
+			return "ICCLab";
+		case PHOTOMETRIC_ITULAB:
+			return "ITULab";
+#endif
+		case PHOTOMETRIC_LOGL:
+			return "CIELog2L";
+		case PHOTOMETRIC_LOGLUV:
+			return "CIELog2Luv";
+	}
+	sprintf(buf,"Unknonw Tag %d",pmtc);
+	return buf;
+}
+
+/* Callbacks used to initialise imdi */
+
+/* Information needed from a profile */
+struct _profinfo {
+	char name[100];
+	icmFile *fp;
+	icc *c;
+	icmHeader *h;
+	icRenderingIntent intent;
+	icmLuBase *luo;				/* Base Lookup type object */
+	icmLuAlgType alg;			/* Type of lookup algorithm */
+	int chan;					/* Device channels */
+}; typedef struct _profinfo profinfo;
+
+/* Context for imdi setup callbacks */
+typedef struct {
+	/* Overall parameters */
+	int verb;			/* Non-zero if verbose */
+	icColorSpaceSignature ins, outs;	/* Input/Output spaces */
+	int id, od;			/* Input/Output dimensions */
+	int icombine;		/* Non-zero if input curves are to be combined */
+	int ocombine;		/* Non-zero if output curves are to be combined */
+	int link;			/* Non-zero if input and output profiles are to be linked */
+
+	profinfo dev;		/* Device link profile */
+	profinfo in;		/* Device to PCS profile */
+	profinfo out;		/* PCS to Device profile */
+} sucntx;
+
+/* Input curve function */
+double input_curve(
+	void *cntx,
+	int ch,
+	double in_val
+) {
+	sucntx *rx = (sucntx *)cntx;
+	int i;
+	double vals[MAX_CHAN];
+
+	if (rx->icombine)
+		return in_val;
+
+	if (rx->link) {
+
+		for (i = 0; i < rx->id; i++)
+			vals[i] = 0.0;
+		vals[ch] = in_val;
+
+		switch(rx->in.alg) {
+		    case icmMonoFwdType: {
+				icmLuMono *lu = (icmLuMono *)rx->in.luo;	/* Safe to coerce */
+				lu->fwd_curve(lu, vals, vals);
+				break;
+			}
+		    case icmMatrixFwdType: {
+				icmLuMatrix *lu = (icmLuMatrix *)rx->in.luo;	/* Safe to coerce */
+				lu->fwd_curve(lu, vals, vals);
+				break;
+			}
+		    case icmLutType: {
+				icmLuLut *lu = (icmLuLut *)rx->in.luo;	/* Safe to coerce */
+				/* Since not PCS, in_abs and matrix cannot be valid, */
+				/* so input curve on own is ok to use. */
+				lu->input(lu, vals, vals);
+				break;
+			}
+			default:
+				error("Unexpected algorithm type in input curve");
+		}
+	} else {
+		icmLuLut *lu = (icmLuLut *)rx->dev.luo;		/* Safe to coerce */
+	
+		for (i = 0; i < rx->id; i++)
+			vals[i] = 0.0;
+		vals[ch] = in_val;
+	
+		/* Since input not PCS, in_abs and matrix cannot be valid, */
+		/* so input curve on own is ok to use. */
+		lu->input(lu, vals, vals);
+	
+	}
+	return vals[ch];
+}
+
+/* Multi-dim table function */
+void md_table(
+void *cntx,
+double *out_vals,
+double *in_vals
+) {
+	sucntx *rx = (sucntx *)cntx;
+	double pcsv[3];
+	int i;
+	
+	if (rx->link) {
+		double vals[MAX_CHAN];
+
+		switch(rx->in.alg) {
+		    case icmMonoFwdType: {
+				icmLuMono *lu = (icmLuMono *)rx->in.luo;	/* Safe to coerce */
+				if (rx->icombine) {
+					lu->fwd_curve(lu, vals, in_vals);
+					lu->fwd_map(lu, pcsv, vals);
+				} else {
+					lu->fwd_map(lu, pcsv, in_vals);
+				}
+				lu->fwd_abs(lu, pcsv, pcsv);
+				break;
+			}
+		    case icmMatrixFwdType: {
+				icmLuMatrix *lu = (icmLuMatrix *)rx->in.luo;	/* Safe to coerce */
+				if (rx->icombine) {
+					lu->fwd_curve(lu, vals, in_vals);
+					lu->fwd_matrix(lu, pcsv, vals);
+				} else {
+					lu->fwd_matrix(lu, pcsv, in_vals);
+				}
+				lu->fwd_abs(lu, pcsv, pcsv);
+				break;
+			}
+		    case icmLutType: {
+				icmLuLut *lu = (icmLuLut *)rx->in.luo;	/* Safe to coerce */
+				if (rx->icombine) {
+					/* Since not PCS, in_abs and matrix cannot be valid, */
+					/* so input curve on own is ok to use. */
+					lu->input(lu, vals, in_vals);
+					lu->clut(lu, pcsv, vals);
+				} else {
+					lu->clut(lu, pcsv, in_vals);
+				}
+				lu->output(lu, pcsv, pcsv);
+				lu->out_abs(lu, pcsv, pcsv);
+				break;
+			}
+			default:
+				error("Unexpected algorithm type in clut lookup");
+		}
+
+		switch(rx->out.alg) {
+		    case icmMonoBwdType: {
+				icmLuMono *lu = (icmLuMono *)rx->out.luo;	/* Safe to coerce */
+				lu->bwd_abs(lu, pcsv, pcsv);
+				lu->bwd_map(lu, out_vals, pcsv);
+				if (rx->ocombine) {
+					lu->bwd_curve(lu, out_vals, out_vals);
+				}
+				break;
+			}
+		    case icmMatrixBwdType: {
+				icmLuMatrix *lu = (icmLuMatrix *)rx->out.luo;	/* Safe to coerce */
+				lu->bwd_abs(lu, pcsv, pcsv);
+				lu->bwd_matrix(lu, out_vals, pcsv);
+				if (rx->ocombine) {
+					lu->bwd_curve(lu, out_vals, out_vals);
+				}
+				break;
+			}
+		    case icmLutType: {
+				icmLuLut *lu = (icmLuLut *)rx->out.luo;	/* Safe to coerce */
+				lu->in_abs(lu, pcsv, pcsv);
+				lu->matrix(lu, pcsv, pcsv);
+				lu->input(lu, pcsv, pcsv);
+				lu->clut(lu, out_vals, pcsv);
+				if (rx->ocombine) {
+					lu->output(lu, out_vals, out_vals);
+					/* Since not PCS, out_abs is never used */
+				}
+				break;
+			}
+
+			default:
+				error("Unexpected algorithm type in clut lookup");
+		}
+	} else {
+		icmLuLut *lu = (icmLuLut *)rx->dev.luo;		/* Safe to coerce */
+
+		if (rx->icombine && rx->ocombine) {
+			lu->lookup((icmLuBase *)lu, out_vals, in_vals);		/* Do everything here */
+		} else {
+			lu->clut(lu, out_vals, in_vals);
+		}
+	}
+}
+
+/* Output curve function */
+double output_curve(
+void *cntx,
+int ch,
+double in_val
+) {
+	sucntx *rx = (sucntx *)cntx;
+	int i;
+	double vals[MAX_CHAN];
+	
+	if (rx->ocombine)
+		return in_val;
+
+	if (rx->link) {
+		for (i = 0; i < rx->od; i++)
+			vals[i] = 0.0;
+		vals[ch] = in_val;
+	
+		switch(rx->out.alg) {
+		    case icmMonoBwdType: {
+				icmLuMono *lu = (icmLuMono *)rx->out.luo;	/* Safe to coerce */
+				lu->bwd_curve(lu, vals, vals);
+				break;
+			}
+		    case icmMatrixBwdType: {
+				icmLuMatrix *lu = (icmLuMatrix *)rx->out.luo;	/* Safe to coerce */
+				lu->bwd_curve(lu, vals, vals);
+				break;
+			}
+		    case icmLutType: {
+				icmLuLut *lu = (icmLuLut *)rx->out.luo;	/* Safe to coerce */
+				lu->output(lu, vals, vals);
+				/* Since not PCS, out_abs is never used */
+				break;
+			}
+			default:
+				error("Unexpected algorithm type in devop_devo()");
+		}
+
+	} else {
+		icmLuLut *lu = (icmLuLut *)rx->dev.luo;		/* Safe to coerce */
+	
+		for (i = 0; i < rx->od; i++)
+			vals[i] = 0.0;
+		vals[ch] = in_val;
+	
+		/* Since output not PCS, out_abs cannot be valid, */
+		lu->output(lu, vals, vals);
+	
+	}
+	return vals[ch];
+}
+
+
+int
+main(int argc, char *argv[]) {
+	int fa,nfa;				/* argument we're looking at */
+	char in_name[100];		/* Raster file name */
+	char out_name[100];		/* Raster file name */
+	int slow = 0;
+	int check = 0;
+	int i, rv = 0;
+
+	TIFF *rh = NULL, *wh = NULL;
+	int x, y, width, height;					/* Size of image */
+	uint16 samplesperpixel, bitspersample;
+	int no_pmtc;								/* Number of input photometrics */
+	uint16 photometric, pmtc[10];				/* Photometrics of input file, and input profile */
+	uint16 pconfig;								/* Planar configuration */
+	uint16 resunits;
+	float resx, resy;
+	tdata_t *inbuf, *outbuf, *checkbuf;
+
+	/* IMDI */
+	imdi *s = NULL;
+	sucntx su;		/* Setup context */
+	unsigned char *inp[MAX_CHAN];
+	unsigned char *outp[MAX_CHAN];
+	int clutres = 33;
+
+	/* Error check */
+	int mxerr = 0;
+	double avgerr = 0.0;
+	double avgcount = 0.0;
+
+	if (argc < 2)
+		usage();
+
+	su.verb = 0;
+	su.icombine = 0;
+	su.ocombine = 0;
+	su.link = 0;
+	su.in.intent = icmDefaultIntent;
+	su.out.intent = icmDefaultIntent;
+
+	/* Process the arguments */
+	for(fa = 1;fa < argc;fa++) {
+		nfa = fa;					/* skip to nfa if next argument is used */
+		if (argv[fa][0] == '-')	{	/* Look for any flags */
+			char *na = NULL;		/* next argument after flag, null if none */
+
+			if (argv[fa][2] != '\000')
+				na = &argv[fa][2];		/* next is directly after flag */
+			else {
+				if ((fa+1) < argc) {
+					if (argv[fa+1][0] != '-') {
+						nfa = fa + 1;
+						na = argv[nfa];		/* next is seperate non-flag argument */
+					}
+				}
+			}
+
+			if (argv[fa][1] == '?')
+				usage();
+
+			/* Slow, Precise */
+			else if (argv[fa][1] == 'p' || argv[fa][1] == 'P') {
+				slow = 1;
+			}
+
+			/* Combine per channel curves */
+			else if (argv[fa][1] == 'c' || argv[fa][1] == 'C') {
+				su.icombine = 1;
+				su.ocombine = 1;
+			}
+
+			/* Check curves */
+			else if (argv[fa][1] == 'k' || argv[fa][1] == 'K') {
+				check = 1;
+			}
+
+			/* Link profiles */
+			else if (argv[fa][1] == 'l' || argv[fa][1] == 'L') {
+				su.link = 1;
+			}
+
+			/* Input profile Intent */
+			else if (argv[fa][1] == 'i' || argv[fa][1] == 'I') {
+				fa = nfa;
+				if (na == NULL) usage();
+    			switch (na[0]) {
+					case 'p':
+					case 'P':
+						su.in.intent = icPerceptual;
+						break;
+					case 'r':
+					case 'R':
+						su.in.intent = icRelativeColorimetric;
+						break;
+					case 's':
+					case 'S':
+						su.in.intent = icSaturation;
+						break;
+					case 'a':
+					case 'A':
+						su.in.intent = icAbsoluteColorimetric;
+						break;
+					default:
+						usage();
+				}
+			}
+
+			/* Output profile Intent */
+			else if (argv[fa][1] == 'o' || argv[fa][1] == 'O') {
+				fa = nfa;
+				if (na == NULL) usage();
+    			switch (na[0]) {
+					case 'p':
+					case 'P':
+						su.out.intent = icPerceptual;
+						break;
+					case 'r':
+					case 'R':
+						su.out.intent = icRelativeColorimetric;
+						break;
+					case 's':
+					case 'S':
+						su.out.intent = icSaturation;
+						break;
+					case 'a':
+					case 'A':
+						su.out.intent = icAbsoluteColorimetric;
+						break;
+					default:
+						usage();
+				}
+			}
+
+			/* Verbosity */
+			else if (argv[fa][1] == 'v' || argv[fa][1] == 'V') {
+				su.verb = 1;
+			}
+
+			else 
+				usage();
+		} else
+			break;
+	}
+
+	if (su.link) {
+		if (fa >= argc || argv[fa][0] == '-') usage();
+		strcpy(su.in.name,argv[fa++]);
+
+		if (fa >= argc || argv[fa][0] == '-') usage();
+		strcpy(su.out.name,argv[fa++]);
+	} else {
+		if (fa >= argc || argv[fa][0] == '-') usage();
+		strcpy(su.dev.name,argv[fa++]);
+	}
+
+	if (fa >= argc || argv[fa][0] == '-') usage();
+	strcpy(in_name,argv[fa++]);
+
+	if (fa >= argc || argv[fa][0] == '-') usage();
+	strcpy(out_name,argv[fa++]);
+
+	/* - - - - - - - - - - - - - - - - */
+
+	if (su.link) {
+		icColorSpaceSignature natpcs;
+
+		/* Open up the input device profile for reading */
+		if ((su.in.fp = new_icmFileStd_name(su.in.name,"r")) == NULL)
+			error ("Can't open file '%s'",su.in.name);
+		
+		if ((su.in.c = new_icc()) == NULL)
+			error ("Creation of Input profile ICC object failed");
+	
+		/* Read header etc. */
+		if ((rv = su.in.c->read(su.in.c,su.in.fp,0)) != 0)
+			error ("%d, %s on file '%s'",rv,su.in.c->err,su.in.name);
+		su.in.h = su.in.c->header;
+	
+		/* Check that it is a suitable device input icc */
+		if (su.in.h->deviceClass != icSigInputClass
+		 && su.in.h->deviceClass != icSigDisplayClass
+		 && su.in.h->deviceClass != icSigOutputClass
+		 && su.in.h->deviceClass != icSigColorSpaceClass)	/* For sRGB etc. */
+			error("Input profile isn't a device profile");
+	
+		/* Get a conversion object */
+		if ((su.in.luo = su.in.c->get_luobj(su.in.c, icmFwd, su.in.intent,
+		                                             icSigLabData, icmLuOrdNorm)) == NULL)
+			error ("%d, %s for profile '%s'",su.in.c->errc, su.in.c->err, su.in.name);
+	
+		/* Get details of conversion (Arguments may be NULL if info not needed) */
+		su.in.luo->spaces(su.in.luo, &su.ins, &su.id, NULL, NULL, &su.in.alg, NULL, NULL, NULL);
+
+		/* Get native PCS space */
+		su.in.luo->lutspaces(su.in.luo, NULL, NULL, NULL, NULL, &natpcs);
+	
+		if (natpcs == icSigXYZData) {
+			su.icombine = 1;			/* XYZ is to non-linear to be a benefit */
+		}
+
+		/* Open up the output device profile for reading */
+		if ((su.out.fp = new_icmFileStd_name(su.out.name,"r")) == NULL)
+			error ("Can't open file '%s'",su.out.name);
+		
+		if ((su.out.c = new_icc()) == NULL)
+			error ("Creation of Output profile ICC object failed");
+	
+		/* Read header etc. */
+		if ((rv = su.out.c->read(su.out.c,su.out.fp,0)) != 0)
+			error ("%d, %s on file '%s'",rv,su.out.c->err,su.out.name);
+		su.out.h = su.out.c->header;
+	
+		/* Check that it is a suitable device output icc */
+		if (su.out.h->deviceClass != icSigInputClass
+		 && su.out.h->deviceClass != icSigDisplayClass
+		 && su.out.h->deviceClass != icSigOutputClass
+		 && su.out.h->deviceClass != icSigColorSpaceClass)	/* For sRGB etc. */
+			error("Output profile isn't a device profile");
+	
+		/* Get a conversion object */
+		if ((su.out.luo = su.out.c->get_luobj(su.out.c, icmBwd, su.out.intent,
+		                                             icSigLabData, icmLuOrdNorm)) == NULL)
+			error ("%d, %s for profile '%s'",su.out.c->errc, su.out.c->err, su.out.name);
+	
+		/* Get details of conversion (Arguments may be NULL if info not needed) */
+		su.out.luo->spaces(su.out.luo, NULL, NULL, &su.outs, &su.od, &su.out.alg, NULL, NULL, NULL);
+
+		/* Get native PCS space */
+		su.out.luo->lutspaces(su.out.luo, NULL, NULL, NULL, NULL, &natpcs);
+	
+		if (natpcs == icSigXYZData) {
+			su.ocombine = 1;			/* XYZ is to non-linear to be a benefit */
+		}
+
+		/* See discussion in imdi/imdi_gen.c for ideal numbers */
+		/* Use "high quality" resolution numbers */
+		switch (su.id) {
+			case 0:
+				error ("Illegal number of input chanels");
+			case 1:
+	  		  	clutres = 256;
+				break;
+			case 2:
+	  		  	clutres = 256;
+				break;
+			case 3:
+	  		  	clutres = 33;
+				break;
+			case 4:
+  		  		clutres = 18;
+				break;
+			case 5:
+  		  		clutres = 16;
+				break;
+			case 6:
+  		  		clutres = 9;
+				break;
+			case 7:
+  		  		clutres = 7;
+				break;
+			case 8:
+  		  		clutres = 6;
+				break;
+			deault: /* > 8 chan */
+				clutres = 3;
+				break;
+		}	
+
+	} else {
+		icmLut *lut;						/* ICC LUT table */
+		icmLuLut *luluo;					/* LUT lookup object */
+
+		/* Open up the device link profile for reading */
+		if ((su.dev.fp = new_icmFileStd_name(su.dev.name,"r")) == NULL)
+			error ("Can't open file '%s'",su.dev.name);
+	
+		if ((su.dev.c = new_icc()) == NULL)
+			error ("Creation of ICC object failed");
+	
+		if ((rv = su.dev.c->read(su.dev.c, su.dev.fp, 0)) != 0)
+			error ("%d, %s",rv,su.dev.c->err);
+		su.dev.h = su.dev.c->header;
+
+		if (su.verb) {
+			icmFile *op;
+			if ((op = new_icmFileStd_fp(stdout)) == NULL)
+				error ("Can't open stdout");
+			su.dev.h->dump(su.dev.h, op, 1);
+			op->del(op);
+		}
+
+		/* Check that the profile is appropriate */
+		if (su.dev.h->deviceClass != icSigLinkClass)
+			error("Profile isn't a device link profile");
+
+		/* Get a conversion object */
+		if ((su.dev.luo = su.dev.c->get_luobj(su.dev.c, icmFwd, icmDefaultIntent,
+		                                             icmSigDefaultData, icmLuOrdNorm)) == NULL)
+			error ("%d, %s",su.dev.c->errc, su.dev.c->err);
+	
+		/* Get details of conversion (Arguments may be NULL if info not needed) */
+		su.dev.luo->spaces(su.dev.luo, &su.ins, &su.id, &su.outs, &su.od, &su.dev.alg, NULL, NULL, NULL);
+
+		if (su.dev.alg != icmLutType)
+			error ("DeviceLink profile doesn't have Lut !");
+
+		luluo = (icmLuLut *)su.dev.luo;		/* Safe to coerce */
+		luluo->get_info(luluo, &lut, NULL, NULL, NULL);	/* Get some details */
+		clutres = lut->clutPoints;			/* Desired table resolution */
+	}
+
+	/* - - - - - - - - - - - - - - - */
+	/* Open up input tiff file ready for reading */
+	/* Got arguments, so setup to process the file */
+	if ((rh = TIFFOpen(in_name, "r")) == NULL)
+		error("error opening read file '%s'",in_name);
+
+	TIFFGetField(rh, TIFFTAG_IMAGEWIDTH,  &width);
+	TIFFGetField(rh, TIFFTAG_IMAGELENGTH, &height);
+
+	TIFFGetField(rh, TIFFTAG_BITSPERSAMPLE, &bitspersample);
+	if (bitspersample != 8 && bitspersample != 16) {
+		error("TIFF Input file must be 8 or 16 bit/channel");
+	}
+
+	TIFFGetField(rh, TIFFTAG_PHOTOMETRIC, &photometric);
+	if  ((no_pmtc = ColorSpaceSignature2TiffPhotometric(pmtc, su.ins)) == 0)
+		error("ICC  input colorspace '%s' can't be handled by a TIFF file!",
+		      icm2str(icmColorSpaceSignature, su.ins));
+	for (i = 0; i < no_pmtc; i++) {
+		if (pmtc[i] == photometric)
+			break;				/* Matches */
+	}
+	if (i >= no_pmtc) {
+		switch (no_pmtc) {
+			case 1:
+				error("ICC input colorspace '%s' doesn't match TIFF photometric '%s'!",
+			      icm2str(icmColorSpaceSignature, su.ins), Photometric2str(pmtc[0]));
+			case 2:
+				error("ICC input colorspace '%s' doesn't match TIFF photometric '%s' or '%s'!",
+			      icm2str(icmColorSpaceSignature, su.ins), Photometric2str(pmtc[0]),
+			      Photometric2str(pmtc[1]));
+			default:
+				error("ICC input colorspace '%s' doesn't match TIFF photometric '%s', '%s' or '%s'!",
+			      icm2str(icmColorSpaceSignature, su.ins), Photometric2str(pmtc[0]),
+			      Photometric2str(pmtc[1]), Photometric2str(pmtc[2]));
+		}
+	}
+
+	TIFFGetField(rh, TIFFTAG_SAMPLESPERPIXEL, &samplesperpixel);
+	if (su.id != samplesperpixel)
+		error ("TIFF Input file has %d input channels mismatched to colorspace '%s'",
+		       samplesperpixel, icm2str(icmColorSpaceSignature, su.ins));
+
+	TIFFGetField(rh, TIFFTAG_PLANARCONFIG, &pconfig);
+	if (pconfig != PLANARCONFIG_CONTIG)
+		error ("TIFF Input file must be planar");
+
+	TIFFGetField(rh, TIFFTAG_RESOLUTIONUNIT, &resunits);
+	TIFFGetField(rh, TIFFTAG_XRESOLUTION, &resx);
+	TIFFGetField(rh, TIFFTAG_YRESOLUTION, &resy);
+
+	/* - - - - - - - - - - - - - - - */
+	if ((wh = TIFFOpen(out_name, "w")) == NULL)
+		error("Can\'t create TIFF file '%s'!",out_name);
+	
+	TIFFSetField(wh, TIFFTAG_IMAGEWIDTH,  width);
+	TIFFSetField(wh, TIFFTAG_IMAGELENGTH, height);
+	TIFFSetField(wh, TIFFTAG_ORIENTATION, ORIENTATION_TOPLEFT);
+	TIFFSetField(wh, TIFFTAG_SAMPLESPERPIXEL, su.od);
+	TIFFSetField(wh, TIFFTAG_BITSPERSAMPLE, bitspersample);
+	TIFFSetField(wh, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
+	if  ((no_pmtc = ColorSpaceSignature2TiffPhotometric(pmtc, su.outs)) == 0)
+		error("TIFF file can't handle output colorspace '%s'!",
+		      icm2str(icmColorSpaceSignature, su.outs));
+	TIFFSetField(wh, TIFFTAG_PHOTOMETRIC, pmtc[0]);	/* Use first returned */
+	if (pmtc[0] == PHOTOMETRIC_SEPARATED) {
+		int iset;
+		int inlen;
+		char *inames;
+		iset = ColorSpaceSignature2TiffInkset(su.outs, &inlen, &inames);
+		if (iset != 0xffff && inlen > 0 && inames != NULL) {
+			TIFFSetField(wh, TIFFTAG_INKSET, iset);
+			if (inames != NULL) {
+				TIFFSetField(wh, TIFFTAG_INKNAMES, inlen, inames);
+			}
+		}
+	}
+	TIFFSetField(wh, TIFFTAG_COMPRESSION, COMPRESSION_NONE);
+	if (resunits) {
+		TIFFSetField(wh, TIFFTAG_RESOLUTIONUNIT, resunits);
+		TIFFSetField(wh, TIFFTAG_XRESOLUTION, resx);
+		TIFFSetField(wh, TIFFTAG_YRESOLUTION, resy);
+	}
+	TIFFSetField(wh, TIFFTAG_IMAGEDESCRIPTION, "Color corrected by Argyll");
+
+	/* - - - - - - - - - - - - - - - */
+	/* Setup the imdi */
+
+	if (!slow) {
+		s = new_imdi(
+			su.id,			/* Number of input dimensions */
+			su.od,			/* Number of output dimensions */
+							/* Input pixel representation */
+			bitspersample == 8 ? pixint8 : pixint16,
+							/* Output pixel representation */
+			0x0,			/* Treat every channel as unsigned */
+			bitspersample == 8 ? pixint8 : pixint16,
+			0x0,			/* Treat every channel as unsigned */
+			clutres,		/* Desired table resolution */
+			input_curve,	/* Callback functions */
+			md_table,
+			output_curve,
+			(void *)&su		/* Context to callbacks */
+		);
+	
+		if (s == NULL)
+			error("new_imdi failed");
+	}
+
+	/* - - - - - - - - - - - - - - - */
+	/* Process colors to translate */
+	/* (Should fix this to process a group of lines at a time ?) */
+
+	inbuf  = _TIFFmalloc(TIFFScanlineSize(rh));
+	outbuf = _TIFFmalloc(TIFFScanlineSize(wh));
+	if (check)
+		checkbuf = _TIFFmalloc(TIFFScanlineSize(wh));
+
+	inp[0] = (unsigned char *)inbuf;
+	outp[0] = (unsigned char *)outbuf;
+
+	if (!slow) {		/* Fast */
+		for (y = 0; y < height; y++) {
+
+			/* Read in the next line */
+			if (TIFFReadScanline(rh, inbuf, y, 0) < 0)
+				error ("Failed to read TIFF line %d",y);
+
+			/* Do fast conversion */
+			s->interp(s, (void **)outp, (void **)inp, width);
+			
+			if (check) {
+				/* Do floating point conversion */
+				for (x = 0; x < width; x++) {
+					int i;
+					double in[MAX_CHAN], out[MAX_CHAN];
+					
+					if (bitspersample == 8)
+						for (i = 0; i < su.id; i++)
+							in[i] = ((unsigned char *)inbuf)[x * su.id + i]/255.0;
+					else
+						for (i = 0; i < su.id; i++)
+							in[i] = ((unsigned short *)inbuf)[x * su.id + i]/65535.0;
+					
+					if (su.link) {
+						if ((rv = su.in.luo->lookup(su.in.luo, out, in)) > 1)
+							error ("%d, %s",su.dev.c->errc,su.dev.c->err);
+						if ((rv = su.out.luo->lookup(su.out.luo, out, out)) > 1)
+							error ("%d, %s",su.dev.c->errc,su.dev.c->err);
+					} else {
+						if ((rv = su.dev.luo->lookup(su.dev.luo, out, in)) > 1)
+							error ("%d, %s",su.dev.c->errc,su.dev.c->err);
+					}
+	
+					if (bitspersample == 8)
+						for (i = 0; i < su.od; i++)
+							((unsigned char *)checkbuf)[x * su.od + i] = (int)(out[i] * 255.0 + 0.5);
+					else
+						for (i = 0; i < su.od; i++)
+							((unsigned short *)checkbuf)[x * su.od + i] = (int)(out[i] * 65535.0 + 0.5);
+				}
+				/* Compute the errors */
+				for (x = 0; x < (width * su.od); x++) {
+					int err;
+					if (bitspersample == 8)
+						err = ((unsigned char *)outbuf)[x] - ((unsigned char *)checkbuf)[x];
+					else
+						err = ((unsigned short *)outbuf)[x] - ((unsigned short *)checkbuf)[x];
+					if (err < 0)
+						err = -err;
+					if (err > mxerr)
+						mxerr = err;
+					avgerr += (double)err;
+					avgcount++;
+				}
+			}
+				
+			if (TIFFWriteScanline(wh, outbuf, y, 0) < 0)
+				error ("Failed to write TIFF line %d",y);
+
+		}
+
+	} else {	/* Slow but precise */
+		if (bitspersample == 8) {
+			for (y = 0; y < height; y++) {
+
+				/* Read in the next line */
+				if (TIFFReadScanline(rh, inbuf, y, 0) < 0)
+					error ("Failed to read TIFF line %d",y);
+
+				/* Do floating point conversion */
+				for (x = 0; x < width; x++) {
+					int i;
+					double in[MAX_CHAN], out[MAX_CHAN];
+					
+					for (i = 0; i < su.id; i++) {
+						in[i] = ((unsigned char *)inbuf)[x * su.id + i]/255.0;
+					}
+					
+					if (su.link) {
+						if ((rv = su.in.luo->lookup(su.in.luo, out, in)) > 1)
+							error ("%d, %s",su.dev.c->errc,su.dev.c->err);
+						if ((rv = su.out.luo->lookup(su.out.luo, out, out)) > 1)
+							error ("%d, %s",su.dev.c->errc,su.dev.c->err);
+					} else {
+						if ((rv = su.dev.luo->lookup(su.dev.luo, out, in)) > 1)
+							error ("%d, %s",su.dev.c->errc,su.dev.c->err);
+					}
+
+					for (i = 0; i < su.od; i++) {
+						double outi = out[i];
+						if (outi < 0.0)			/* Protect against sillies */
+							outi = 0.0;
+						else if (outi > 1.0)
+							outi = 1.0;
+						((unsigned char *)outbuf)[x * su.od + i] = (int)(outi * 255.0 + 0.5);
+					}
+				}
+				if (TIFFWriteScanline(wh, outbuf, y, 0) < 0)
+					error ("Failed to write TIFF line %d",y);
+			}
+		} else if (bitspersample == 16) {
+			for (y = 0; y < height; y++) {
+
+				/* Read in the next line */
+				if (TIFFReadScanline(rh, inbuf, y, 0) < 0)
+					error ("Failed to read TIFF line %d",y);
+
+				/* Do floating point conversion */
+				for (x = 0; x < width; x++) {
+					int i;
+					double in[MAX_CHAN], out[MAX_CHAN];
+					
+					for (i = 0; i < su.id; i++) {
+						in[i] = ((unsigned short *)inbuf)[x * su.id + i]/65535.0;
+					}
+					
+					if (su.link) {
+						if ((rv = su.in.luo->lookup(su.in.luo, out, in)) > 1)
+							error ("%d, %s",su.dev.c->errc,su.dev.c->err);
+						if ((rv = su.out.luo->lookup(su.out.luo, out, out)) > 1)
+							error ("%d, %s",su.dev.c->errc,su.dev.c->err);
+					} else {
+						if ((rv = su.dev.luo->lookup(su.dev.luo, out, in)) > 1)
+							error ("%d, %s",su.dev.c->errc,su.dev.c->err);
+					}
+
+					for (i = 0; i < su.od; i++) {
+						double outi = out[i];
+						if (outi < 0.0)			/* Protect against sillies */
+							outi = 0.0;
+						else if (outi > 1.0)
+							outi = 1.0;
+						((unsigned short *)outbuf)[x * su.od + i] = (int)(outi * 65535.0 + 0.5);
+					}
+				}
+				if (TIFFWriteScanline(wh, outbuf, y, 0) < 0)
+					error ("Failed to write TIFF line %d",y);
+			}
+		}
+	}
+
+	if (check) {
+		printf("Worst error = %d bits, average error = %f bits\n", mxerr, avgerr/avgcount);
+		if (bitspersample == 8)
+			printf("Worst error = %f%%, average error = %f%%\n",
+			       mxerr/2.55, avgerr/(2.55 * avgcount));
+		else
+			printf("Worst error = %f%%, average error = %f%%\n",
+			       mxerr/655.35, avgerr/(655.35 * avgcount));
+	}
+
+	/* Done with lookup object */
+	if (s != NULL)
+		s->done(s);
+
+	if (su.link) {
+		su.in.luo->del(su.in.luo);
+		su.in.c->del(su.in.c);
+		su.in.fp->del(su.in.fp);
+		su.out.luo->del(su.out.luo);
+		su.out.c->del(su.out.c);
+		su.out.fp->del(su.out.fp);
+	} else {
+		su.dev.luo->del(su.dev.luo);
+		su.dev.c->del(su.dev.c);
+		su.dev.fp->del(su.dev.fp);
+	}
+
+	_TIFFfree(inbuf);
+	_TIFFfree(outbuf);
+	if (check)
+		_TIFFfree(checkbuf);
+
+	TIFFClose(rh);		/* Close Input file */
+	TIFFClose(wh);		/* Close Output file */
+
+	return 0;
+}
+
+
+/* Basic printf type error() and warning() routines */
+
+void
+error(char *fmt, ...)
+{
+	va_list args;
+
+	fprintf(stderr,"cctiff: Error - ");
+	va_start(args, fmt);
+	vfprintf(stderr, fmt, args);
+	va_end(args);
+	fprintf(stderr, "\n");
+	exit (-1);
+}
+
+void
+warning(char *fmt, ...)
+{
+	va_list args;
+
+	fprintf(stderr,"cctiff: Warning - ");
+	va_start(args, fmt);
+	vfprintf(stderr, fmt, args);
+	va_end(args);
+	fprintf(stderr, "\n");
+}

Added: trunk/gs/imdi/cgen.c
===================================================================
--- trunk/gs/imdi/cgen.c	2006-09-11 07:02:18 UTC (rev 7033)
+++ trunk/gs/imdi/cgen.c	2006-09-11 20:26:01 UTC (rev 7034)
@@ -0,0 +1,1861 @@
+
+/* Integer Multi-Dimensional Interpolation */
+
+/*
+ * Copyright 2000 - 2002 Graeme W. Gill
+ * All rights reserved.
+ *
+ * This material is licenced under the GNU GENERAL PUBLIC LICENCE :-
+ * see the Licence.txt file for licencing details.
+ */
+
+/* 'C' code color transform kernel code generator. */
+
+/*
+   This module generates C code routines which implement
+   an integer multi-channel transform. The input values
+   are read, passed through per channel lookup tables,
+   a multi-dimentional interpolation table, and then
+   a per channel output lookup table, before being written.
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <stdarg.h>
+#include <string.h>
+
+#include "imdi_imp.h"
+#include "imdi_gen.h"
+#include "imdi_tab.h"
+
+#undef VERBOSE
+#undef FORCESORT		/* Use sort algorithm allways */
+
+/*
+ * TTBD:
+ *		Need to implement g->dir
+ *      Haven't used t->it_map[] or t->im_map[].
+ *
+ *
+ */
+
+/* ------------------------------------ */
+/* Context */
+typedef struct {
+	FILE *of;			/* Output file */
+	int indt;			/* Indent */
+
+	/* Other info */
+	genspec *g;			/* Generation specifications */
+	tabspec *t;			/* Table setup data */
+	mach_arch *a;		/* Machine architecture and tuning data */
+
+	/* Code generation information */
+	/* if() conditions are for entry usage */
+
+	/* Pixel read information */
+	int ipt[IXDI];		/* Input pointer types */
+	int nip;			/* Actual number of input pointers, accounting for pint */
+	int chv_bits;		/* Bits in chv temp variable ?? */
+
+	/* Input table entry */
+	int itet;			/* Input table entry type */
+	int itvt;			/* Input table variable type */
+	int itmnb;			/* Input table minimum bits (actual is it_ab) */
+
+	/* Interpolation index */
+	int ixet;			/* Interpolation index entry type */
+	int ixvt;			/* Interpolation index variable type */
+	int ixmnb;			/* Interpolation index minimum bits (actual is ix_ab) */
+	int ixmxres;		/* Interpolation table maximum resolution */
+
+	/* Simplex index: if(!sort && it_xs) */
+	int sxet;			/* Simplex index entry type  */
+	int sxvt;			/* Simplex index variable type */
+	int sxmnb;			/* Simplex index bits minimum (actual is sx_ab) */
+	int sxmxres;		/* Simplex table maximum resolution (0 if sort) */
+
+	/* Combination Weighting + Vertex offset values: if(it_xs && !wo_xs) */
+	int woet;			/* Weighting+offset entry type  */
+	int wovt;			/* Weighting+offset variable type */
+	int womnb;			/* Weighting+offset index bits minimum (actual is wo_ab) */
+
+	/* Weighting value: if(it_xs && wo_xs) */
+	int weet;			/* Weighting entry type  */
+	int wevt;			/* Weighting variable type */
+	int wemnb;			/* Weighting index bits minimum (actual is we_ab) */
+
+	/* Vertex offset value: if(it_xs && wo_xs) */
+	int voet;			/* Vertex offset entry type  */
+	int vovt;			/* Vertex offset variable type */
+	int vomnb;			/* Vertex offset index bits minimum (actual is vo_ab) */
+
+	/* Interpolation table entry: */
+	int imovb;			/* Interpolation output value bits per channel required */
+	int imfvt;			/* Interpolation full entry & variable type */
+	int impvt;			/* Interpolation partial entry variable type */
+
+	/* Interpolation accumulators: */
+	int iaovb;			/* Interpolation output value bits per channel required */
+	int iafvt;			/* Interpolation full entry & variable type */
+	int iapvt;			/* Interpolation partial entry variable type */
+	int ian;			/* Total number of accumulators */
+
+	/* Output table lookup */
+	int otit;			/* Output table index type */
+	int otvt;			/* Output table value type (size is ot_ts bytes) */
+
+	/* Write information */
+	int opt[IXDO];		/* Output pointer types */
+	int nop;			/* Actual number of output pointers, accounting for pint */
+
+} fileo;
+
+void line(fileo *f, char *fmt, ...);	/* Output one line */
+void sline(fileo *f, char *fmt, ...);	/* Output start of line line */
+void mline(fileo *f, char *fmt, ...);	/* Output middle of line */
+void eline(fileo *f, char *fmt, ...);	/* Output end of line */
+void cr(fileo *f) { line(f,""); }		/* Output a blank line */
+void inc(fileo *f) { f->indt++; }		/* Increment the indent level */
+void dec(fileo *f) { f->indt--; }		/* Decrement the indent level */
+/* ------------------------------------ */
+
+int findord(fileo *f, int bits);		/* Find ordinal with bits or more */
+int nord(fileo *f, int ov);				/* Round ordinal type up to natural size */
+int findnord(fileo *f, int bits);		/* Find ordinal with bits, or natural larger */
+int findint(fileo *f, int bits);		/* Find integer with bits or more */
+int nint(fileo *f, int iv);				/* Round integer type up to natural size */
+int findnint(fileo *f, int bits);		/* Find integer with bits, or natural larger */
+static void doheader(fileo *f);
+
+static int calc_bits(int dim, int res);
+static int calc_res(int dim, int bits);
+static int calc_obits(int dim, int res, int esize);
+static int calc_ores(int dim, int bits, int esize);
+
+
+/* return a hexadecimal mask string */
+/* take care of the case when bits >= 32 */
+char *hmask(int bits) {
+	static char buf[20];
+
+	if (bits < 32) {
+		sprintf(buf, "0x%x",(1 << bits)-1);
+	} else if (bits == 32) {
+		return "0xffffffff";
+	} else if (bits == 64) {
+		return "0xffffffffffffffff";
+	} else {	/* Bits > 32 */
+		sprintf(buf, "0x%xffffffff",(1 << (bits-32))-1);
+	}
+	return buf;
+}
+
+/* Generate a source file to implement the specified */
+/* interpolation kernel. Fill in return values and return 0 if OK. */
+/* Return non-zero on error. */
+int gen_c_kernel(
+	genspec *g,				/* Specification of what to generate */
+	mach_arch *a,
+	FILE *fp,				/* File to write to */
+	int index				/* Identification index, 1 = first */
+) {
+	unsigned char kk[] = { 0x43, 0x6F, 0x70, 0x79, 0x72, 0x69, 0x67, 0x68,
+	                       0x74, 0x20, 0x32, 0x30, 0x30, 0x34, 0x20, 0x47,
+	                       0x72, 0x61, 0x65, 0x6D, 0x65, 0x20, 0x57, 0x2E,
+	                       0x20, 0x47, 0x69, 0x6C, 0x6C, 0x00 };
+	fileo f[1];
+	int e, i;
+	tabspec tabsp, *t = &tabsp;
+	int timp = 0;		/* Flag to use temporary imp pointer. */
+						/* Seem to make x86 MSVC++ slower */
+						/* Has no effect on x86 IBMCC */
+
+	sprintf(g->kname, "imdi_k%d",index); /* Kernel routine base name */
+	strcpy(g->kkeys, kk);				 /* Kernel keys for this session */
+	
+	/* Setup the file output context */
+	f->of = fp;
+	f->indt = 0;			/* Start with no indentation */
+	f->g = g;
+	f->t = t;
+	f->a = a;
+
+	if (g->prec == 8) {
+		if (g->id <= 4)
+			t->sort = 0;		/* Implicit sort using simplex table lookup */
+		else
+			t->sort = 1;		/* Explicit sort */
+
+	} else  if (g->prec == 16) {
+		t->sort = 1;			/* Explit sort, no simplex table */
+
+	} else {
+		fprintf(stderr,"Can't cope with requested precision of %d bits\n",g->prec);
+		exit(-1);
+	}
+#ifdef FORCESORT
+	t->sort = 1;
+#endif
+
+	/* Compute input read and input table lookup stuff */
+
+	/* Compute number of input pointers */
+	if (g->in.pint != 0)	/* Pixel interleaved */
+		f->nip = 1;
+	else
+		f->nip = g->id;
+
+	/* Figure out the input pointer types */
+	for (e = 0; e < f->nip; e++) {
+		if ((f->ipt[e] = findord(f, g->in.bpch[e])) < 0) {
+			fprintf(stderr,"Input channel size can't be handled\n");
+			exit(-1);
+		}
+	}
+
+	/* Set a default input channel mapping */
+	for (e = 0; e < g->id; e++)
+		t->it_map[e] = e;
+
+	/* Do the rest of the input table size calculations after figuring */
+	/* out simplex and interpolation table sizes. */
+
+
+	/* Figure out the interpolation multi-dimentional table structure */
+	/* and output accumulation variable sizes. */
+
+	if (g->prec == 8
+	 || g->prec == 16 && a->ords[a->nords-1].bits >= (g->prec * 4)) {
+		int tiby;		/* Total interpolation bytes needed */
+
+		/* We assume that we can normally compute more than one */
+		/* output value at a time, so we need to hold the interpolation */
+		/* output data in the expanded fixed point format in both the */
+		/* table and accumulator. */
+		t->im_cd = 1;
+		f->imovb = g->prec * 2;		/* 16 bits needed for 8 bit precision, */
+		f->iaovb = g->prec * 2;		/* 32 bits needed for 16 bit precision */
+		f->imfvt = a->nords-1;		/* Full variable entry type is biggest available */
+		f->iafvt = a->nords-1;		/* Full variable accum. type is same */
+
+		if (a->ords[f->imfvt].bits < f->imovb) {
+			fprintf(stderr,"Interpolation table entry size can't be handled\n");
+			exit(-1);
+		}
+
+		/* Compute details of table entry sizes, number */
+		tiby = (f->imovb * g->od)/8;				/* Total table bytes needed */
+		t->im_fs = a->ords[f->imfvt].bits/8;		/* Full entry bytes */
+		t->im_fv = (t->im_fs * 8)/f->imovb;			/* output values per full entry . */
+		t->im_fn = tiby/t->im_fs;					/* Number of full entries (may be 0) */
+		t->im_ts = t->im_fn * t->im_fs;				/* Structure size so far */
+		tiby -= t->im_fn * t->im_fs;				/* Remaining bytes */
+
+		if (tiby <= 0) {
+			t->im_pn = 0;		/* No partials */
+			t->im_ps = 0;
+			t->im_pv = 0;
+			f->impvt = 0;
+			f->iapvt = 0;
+
+		} else {
+			t->im_pn = 1;					/* Must be just 1 partial */
+			t->im_pv = (tiby * 8)/f->imovb;	/* Partial holds remaining entries */ 
+
+			if ((f->impvt = findnord(f, tiby * 8)) < 0) {
+				fprintf(stderr,"Can't find partial interp table entry variable size\n");
+				exit(-1);
+			}
+			f->iapvt = f->impvt;
+			t->im_ps = a->ords[f->impvt].bits/8;/* Partial entry bytes */
+
+			if (a->ords[f->imfvt].align)		/* If full entry's need to be aligned */
+				t->im_ts += t->im_fs;			/* Round out struct size by full entry */
+			else
+				t->im_ts += t->im_ps;			/* Round out to natural size */
+		}
+	
+	} else {
+		/* One 16 bit output value per entry + 32 bit accumulator. */
+		/* We can conserve table space by not holding the table data in expanded */
+		/* fixed point format, but expanding it when it is read. */
+		/* Without resorting to compicated code, this restricts us */
+		/* to only computing one output value per accumulator. */
+		t->im_cd = 0;
+		f->imovb = g->prec;			/* Table holds 16 bit entries with no fractions */
+		f->iaovb = g->prec * 2;		/* 32 bits needed for 16 bit precision in comp. */
+
+		if ((f->imfvt = findord(f, f->imovb)) < 0) {
+			fprintf(stderr,"Interpolation table entry size can't be handled\n");
+			exit(-1);
+		}
+		if ((f->iafvt = findord(f, f->iaovb)) < 0) {
+			fprintf(stderr,"Interpolation accumulator size can't be handled\n");
+			exit(-1);
+		}
+
+		/* Compute details of table entry sizes, number */
+		t->im_fs = a->ords[f->imfvt].bits/8;		/* Full entry bytes */
+		t->im_fv = 1;								/* output values per full entry . */
+		t->im_fn = g->od;							/* Number of full entries */
+		t->im_ts = t->im_fn * t->im_fs;				/* Total structure size */
+
+		t->im_pn = 0;		/* No partials */
+		t->im_ps = 0;
+		t->im_pv = 0;
+		f->impvt = 0;
+		f->iapvt = 0;
+	}
+	f->ian = t->im_fn + t->im_pn;			/* Total number of output accumulators */
+
+	/* Figure out how much of the interpolation entry offset to put in the */
+	/* vertex offset value, and how much to make explicit in accessing the */
+	/* interpolation table enty. */
+	if (a->oscale > 0) {		/* We have a scaled index mode */
+		/* Use as much of the scaled index mode as possible */
+		/* and then do the balance by scaling the simplex index entry. */
+		for (t->im_oc = a->oscale; ; t->im_oc >>= 1) {
+			t->vo_om = t->im_ts/t->im_oc;		/* Simplex index multiplier */
+			if ((t->vo_om * t->im_oc) == t->im_ts)
+				break;				/* Got appropriate offset scale */
+		}
+	} else if (a->smmul) {		/* Architecure supports fast small multiply */
+		t->im_oc = t->im_ts;	/* Do scale by structure size explicitly */
+		t->vo_om = 1;			/* Do none in the Simplex index */
+	} else {					/* We have no fast tricks */
+		t->im_oc = 1;			/* Do none explicitly */
+		t->vo_om = t->im_ts;	/* Do all in Simplex index */
+	}
+
+	/* Compute the number of bits needed to hold an index into */
+	/* the interpolation table (index is in terms of table entry size). */
+	/* This value is used to figure out the room needed in the input */
+	/* table to accumulate the interpolation cube base offset value. (IM_O macro) */
+	f->ixmnb = calc_bits(g->id, g->itres);
+
+	/* Set a default output channel mapping */
+	for (e = 0; e < g->od; e++)
+		t->im_map[e] = e;
+
+#ifdef VERBOSE
+	/* Summarise the interpolation table arrangements */
+	printf("\n");
+	printf("Interpolation table structure:\n");
+	printf("  Minimum bits needed to index table %d\n", f->ixmnb);
+	printf("  Entry total size %d bytes\n", t->im_ts);
+	printf("  Simplex entry offset scale %d\n", t->vo_om);
+	printf("  Explicit entry offset scale %d\n", t->im_oc);
+	printf("  %d full entries, size %d bytes\n", t->im_fn, t->im_fs);
+	printf("  %d partial entries, size %d bytes\n", t->im_pn, t->im_ps);
+	printf("  to hold %d output values of %d bits\n", g->od, f->imovb);
+
+#endif /* VERBOSE */
+
+	/* Number of bits needed for the weighting value */
+	f->wemnb = g->prec+1;	/* Need to hold a weighting factor of 0 - 256 for 8 bits */
+							/* Need to hold a weighting factor of 0 - 65536 for 16 bits */
+
+	/* Variable that would be used to hold it */
+	if ((f->wevt = findnord(f, f->wemnb)) < 0) {
+		fprintf(stderr,"Can't find entry size to hold weighting variable\n");
+		exit(-1);
+	}
+
+	/* Number of bits needed for vertex offset value */
+	f->vomnb = calc_obits(g->id, g->itres, t->vo_om);
+
+	/* Variable that would be used to hold it */
+	if ((f->vovt = findnord(f, f->vomnb)) < 0) {
+		fprintf(stderr,"Can't find entry size to hold vertex offset variable\n");
+		exit(-1);
+	}
+
+	if (t->sort) {
+		/* If we are using an explicit sort, we need to figure how many */
+		/* separate entries we need to use to hold the interpolation index, */
+		/* weighting factor and vertex offset values in the input table. */
+
+		/* First try all three in one entry */
+		if ((f->itet = findord(f, f->ixmnb + f->wemnb + f->vomnb)) >= 0) {/* size to read */
+			int rem;						/* Remainder bits */
+
+			t->it_xs = 0;					/* Combined interp+weight+offset */
+			t->wo_xs = 0;
+			t->it_ab = a->ords[f->itet].bits;	/* Bits in combined input entry */
+			rem = t->it_ab - f->ixmnb - f->wemnb - f->vomnb; /* Spair bits */
+			t->we_ab = f->wemnb;				/* Get minimum weight bits */
+			t->vo_ab = f->vomnb + rem/2;		/* vertex offset index bits actually available */
+			t->ix_ab = t->it_ab - t->vo_ab - t->we_ab;	/* interp index bits actually available */
+			t->wo_ab = t->we_ab + t->vo_ab;		/* Weight & offset total bits */
+			t->it_ts = a->ords[f->itet].bits/8;	/* total size in bytes */
+			f->itvt = nord(f, f->itet);			/* Variable type */
+
+			if ((f->wovt = findnord(f, t->we_ab + t->vo_ab)) < 0) {
+				fprintf(stderr,"Can't find variable size to hold weight/offset\n");
+				exit(-1);
+			}
+			if ((f->wevt = findnord(f, t->we_ab)) < 0) {
+				fprintf(stderr,"Can't find variable size to hold weighting factor\n");
+				exit(-1);
+			}
+			if ((f->vovt = findnord(f, t->vo_ab)) < 0) {
+				fprintf(stderr,"Can't find variable size to hold vertex offset index\n");
+				exit(-1);
+			}
+			if ((f->ixvt = findnord(f, t->ix_ab)) < 0) {
+				fprintf(stderr,"Interp index variable size can't be handled\n");
+				exit(-1);
+			}
+		} else {	/* Interp index will be a separate entry */
+			int wit, oft, bigt;		/* weighting type, offset type, biggest type */
+			int combt;				/* Combined type */
+			int sepbits, combits;	/* Total separate, combined bits */
+
+			t->it_xs = 1;				/* Separate interp index and weighting+offset */
+			if ((f->ixet = findord(f, f->ixmnb)) < 0) {
+				fprintf(stderr,"Interp index entry size can't be handled\n");
+				exit(-1);
+			}
+			f->ixvt = nord(f, f->ixet);		/* Variable type */
+			t->ix_ab = a->ords[f->ixet].bits;
+			t->ix_es = t->ix_ab/8;
+			t->ix_eo = 0;
+			t->it_ts = t->ix_es;			/* Input table size so far */
+
+			/* Now figure weighting and vertex offset */
+
+			/* See if we can fit them into separately readable entries, or whether */
+			/* they should be combined to minimise overall table size. */
+
+			if ((wit = findord(f, f->wemnb)) < 0) {
+				fprintf(stderr,"Can't find entry size to hold weighting factor\n");
+				exit(-1);
+			}
+			if ((oft = findord(f, f->vomnb)) < 0) {
+				fprintf(stderr,"Can't find entry size to hold vertex offset index\n");
+				exit(-1);
+			}
+			bigt = wit > oft ? wit : oft;			/* Bigest separate type */
+
+			if ((combt = findord(f, f->wemnb + f->vomnb)) < 0) {/* Combined isn't possible */
+				sepbits = 2 * a->ords[bigt].bits;		/* Total separate bits */
+				combits = sepbits;						/* Force separate entries */
+			} else {
+				sepbits = 2 * a->ords[bigt].bits;		/* Total separate bits */
+				combits = a->ords[combt].bits;			/* Total combined bits */
+			}
+
+			if (sepbits <= combits) {				/* We will use separate entries */
+				t->wo_xs = 1;
+				t->we_es = a->ords[bigt].bits/8;	/* size in bytes for weighting entry */
+				t->we_ab = a->ords[bigt].bits;		/* bits available for weighting */
+				t->we_eo = t->ix_es;				/* Entry offset in input table */
+				t->vo_es = a->ords[bigt].bits/8;	/* size in bytes for vertex offset entry */
+				t->vo_ab = a->ords[bigt].bits;		/* bits available for vertex offset */
+				t->vo_eo = t->ix_es + t->we_es;		/* Entry offset in input table */
+				t->wo_es = t->we_es + t->vo_es;		/* Total entry size for each vertex */
+				t->it_ts += t->we_es + t->vo_es;	/* Total input entry size in bytes */
+
+				f->weet = bigt;				/* Variable type for accessing weighting entry */
+				f->voet = bigt;				/* Variable type for accessing vertex offset entry */
+				f->wevt = nord(f, wit);		/* Variable type for holding weight value */
+				f->vovt = nord(f, oft);		/* Variable type for holding offset value */
+
+			} else {								/* We will combine the two entries */
+				t->wo_xs = 0;
+				t->wo_es = a->ords[combt].bits/8;	/* entry size in bytes for each entry */
+				t->wo_ab = a->ords[combt].bits;		/* bits in weightig + offset */
+				t->we_ab = f->wemnb;				/* bits available for weighting */
+				t->vo_ab = t->wo_ab - t->we_ab;		/* Allow all spare bits to vertex offset */
+				t->wo_eo = t->ix_es;				/* entry offset in input table */
+				t->it_ts += t->wo_es;				/* Final input table size */
+
+				f->woet = combt;			/* Variable type for accessing combined entry */
+				f->wovt = nord(f, combt);	/* Variable type holding weight/offset read value */
+
+				if ((f->wevt = findnord(f, t->we_ab)) < 0) {
+					fprintf(stderr,"Can't find variable size to hold weighting factor\n");
+					exit(-1);
+				}
+				if ((f->vovt = findnord(f, t->vo_ab)) < 0) {
+					fprintf(stderr,"Can't find variable size to hold vertex offset index\n");
+					exit(-1);
+				}
+			}
+		}
+#ifdef VERBOSE
+		/* Summarise the input table arrangements */
+		printf("\n");
+		printf("Input table structure:\n");
+		printf("  Input value re-ordering is [");
+		for (e = 0; e < g->id; e++)
+			printf("%s%d",e > 0 ? " " : "", t->it_map[e]);
+		printf("]\n");
+		printf("  Input table entry size = %d bytes\n",t->it_ts);
+		if (t->it_ix) {
+			printf("  Input table extracts value from read values\n");
+			if (t->wo_xs) {
+				printf("  Separate Interp., Weighting and Offset values\n");
+				printf("  Interp. index is at offset %d, size %d bytes\n",t->ix_eo, t->ix_es);
+				printf("  Weighting is at offset %d, size %d bytes\n",t->we_eo, t->we_es);
+				printf("  Vertex offset is at offset %d, size %d bytes\n",t->vo_eo, t->vo_es);
+			} else {
+				printf("  Separate Interp. index and Weightint+Offset value\n");
+				printf("  Interp. index is at offset %d, size %d bytes\n",t->ix_eo, t->ix_es);
+				printf("  Weighting+Offset is at offset %d, size %d bytes\n",t->wo_eo, t->wo_es);
+				printf("  Weighting     = %d bits\n",t->we_ab);
+				printf("  Vertex offset = %d bits\n",t->vo_ab);
+			}
+		} else {
+			printf("  Combined InterpIndex+Weighting+Voffset values\n");
+			printf("  Values are stored in size %d bytes\n",t->it_ts);
+			printf("  Interp. index = %d bits\n",t->ix_ab);
+			printf("  Weighting     = %d bits\n",t->we_ab);
+			printf("  Vertex offset = %d bits\n",t->vo_ab);
+		}
+#endif /* VERBOSE */
+
+	} else {	/* Simplex table */
+		/* If we are going to use a simplex table, figure out how we */
+		/* will store the weighting value and vertex offset values in it, */
+		/* as well as the size of index we'll need to address it. */
+		int wit, oft, bigt;		/* weighting type, offset type, biggest type */
+		int combt;				/* Combined type */
+		int sepbits, combits;	/* Total separate, combined bits */
+
+		/* See if we can fit them into separately readable entries, or whether */
+		/* they should be combined to minimise overall table size. */
+
+		if ((wit = findord(f, f->wemnb)) < 0) {
+			fprintf(stderr,"Can't find entry size to hold weighting factor\n");
+			exit(-1);
+		}
+		if ((oft = findord(f, f->vomnb)) < 0) {
+			fprintf(stderr,"Can't find entry size to hold vertex offset index\n");
+			exit(-1);
+		}
+		bigt = wit > oft ? wit : oft;			/* Bigest separate type */
+
+		if ((combt = findord(f, f->wemnb + f->vomnb)) < 0) {/* Combined isn't possible */
+			sepbits = 2 * a->ords[bigt].bits;		/* Total separate bits */
+			combits = sepbits;						/* Force separate entries */
+		} else {
+			sepbits = 2 * a->ords[bigt].bits;		/* Total separate bits */
+			combits = a->ords[combt].bits;			/* Total combined bits */
+		}
+
+		if (sepbits <= combits) {				/* We will use separate entries */
+			t->wo_xs = 1;
+			t->we_es = a->ords[bigt].bits/8;	/* size in bytes for weighting entry */
+			t->we_ab = a->ords[bigt].bits;		/* bits available for weighting */
+			t->we_eo = 0;						/* Entry offset in simplex table */
+			t->vo_es = a->ords[bigt].bits/8;	/* size in bytes for vertex offset entry */
+			t->vo_ab = a->ords[bigt].bits;		/* bits available for vertex offset */
+			t->vo_eo = t->we_es;				/* Entry offset in simplex table */
+			t->wo_es = t->we_es + t->vo_es;		/* Total entry size for each vertex */
+			t->sm_ts = (g->id + 1) * (t->we_es + t->vo_es) ;	/* Total size in bytes */
+
+			f->weet = bigt;				/* Variable type for accessing weighting entry */
+			f->voet = bigt;				/* Variable type for accessing vertex offset entry */
+			f->wevt = nord(f, wit);		/* Variable type for holding weight value */
+			f->vovt = nord(f, oft);		/* Variable type for holding offset value */
+
+		} else {								/* We will combine the two entries */
+			t->wo_xs = 0;
+			t->wo_es = a->ords[combt].bits/8;	/* entry size in bytes for each entry */
+			t->wo_ab = a->ords[combt].bits;		/* bits in weightig + offset */
+			t->we_ab = f->wemnb;				/* bits available for weighting */
+			t->vo_ab = t->wo_ab - t->we_ab;		/* Allow all spare bits to vertex offset */
+			t->wo_eo = 0;						/* entry offset in simplex table */
+			t->sm_ts = (g->id + 1) * t->wo_es;	/* Total size in bytes */
+
+			f->woet = combt;			/* Variable type for accessing combined entry */
+			f->wovt = nord(f, combt);	/* Variable type holding weight/offset read value */
+
+			if ((f->wevt = findnord(f, t->we_ab)) < 0) {
+				fprintf(stderr,"Can't find variable size to hold weighting factor\n");
+				exit(-1);
+			}
+			if ((f->vovt = findnord(f, t->vo_ab)) < 0) {
+				fprintf(stderr,"Can't find variable size to hold vertex offset index\n");
+				exit(-1);
+			}
+		}
+
+		/* Compute the number of bits needed to hold an index into */
+		/* the simplex table (index is in terms of table entry size). */
+		/* This value is used to figure out the room needed in the input */
+		/* table to accumulate the simplex cube base offset value. (SW_O macro) */
+		f->sxmnb = calc_bits(g->id, g->stres);
+
+#ifdef VERBOSE
+		/* Summarise the simplex table arrangements */
+		printf("\n");
+		printf("Simplex table structure:\n");
+		printf("  Minimum bits needed to index table %d\n", f->sxmnb);
+		printf("  Total simplex entry size %d bytes to hold %d entries\n",t->sm_ts, g->id+1);
+		if (t->wo_xs) {
+			printf("  Separate entries for offset and weight\n");
+			printf("  Weighting entry size %d bytes\n",t->we_es);
+			printf("  Offset entry size %d bytes\n",t->vo_es);
+		} else {
+			printf("  Combined offset and weight entries in %d bytes\n",t->wo_es);
+			printf("  Weighting entry size %d bits\n",t->we_ab);
+			printf("  Offset entry size %d bits\n",t->vo_ab);
+		}
+		printf("  Vertex offset scale factor %d\n", t->vo_om);
+#endif /* VERBOSE */
+
+		/* We known how big the interpolation and simplex */
+		/* tables indexes are going to be, so complete figuring out */
+		/* how big the input table entries have to be. */
+		if ((f->itet = findord(f, f->sxmnb + f->ixmnb)) >= 0) {/* size to read */
+			int rem;						/* Remainder bits */
+
+			t->it_xs = 0;					/* Combined simplex+interp index */
+
+			t->it_ab = a->ords[f->itet].bits;	/* Bits in combined input entry */
+			rem = t->it_ab - f->sxmnb - f->ixmnb;
+			t->sx_ab = f->sxmnb + rem/2;		/* simplex index bits actually available */
+			t->ix_ab = t->it_ab - t->sx_ab;		/* interp index bits actually available */
+			t->it_ts = a->ords[f->itet].bits/8;	/* total size in bytes */
+			f->itvt = nord(f, f->itet);			/* Variable type */
+
+			if ((f->sxvt = findnord(f, t->sx_ab)) < 0) {
+				fprintf(stderr,"Simplex index variable size can't be handled\n");
+				exit(-1);
+			}
+			if ((f->ixvt = findnord(f, t->ix_ab)) < 0) {
+				fprintf(stderr,"Interp index variable size can't be handled\n");
+				exit(-1);
+			}
+		} else {						/* Separate entries */
+			int bbits;					/* Largest number of bits needed */
+
+			t->it_xs = 1;				/* Separate simplex+interp indexes */
+			bbits = f->sxmnb > f->ixmnb ? f->sxmnb : f->ixmnb;
+
+			/* Allocate same size for both so that total structure size is power of 2 */
+			if ((f->sxet = f->ixet = findord(f, bbits)) < 0) {
+				fprintf(stderr,"Interp/Simplex index entry size can't be handled\n");
+				exit(-1);
+			}
+
+			t->sx_ab = a->ords[f->sxet].bits;		/* Actual bits available */
+			t->sx_es = t->sx_ab/8;					/* Entry size in bytes */
+			t->ix_ab = a->ords[f->ixet].bits;
+			t->ix_es = t->sx_ab/8;
+			t->it_ts = t->sx_es + t->ix_es;		/* total size in bytes */
+			t->sx_eo = 0;						/* simplex index offset in bytes */
+			t->ix_eo = t->sx_es;				/* interp. index offset in bytes */
+			f->sxvt = nord(f, f->sxet);			/* Variable type */
+			f->ixvt = nord(f, f->ixet);			/* Variable type */
+		}
+
+#ifdef VERBOSE
+		/* Summarise the input table arrangements */
+		printf("\n");
+		printf("Input table structure:\n");
+		if (t->it_ix) {
+			printf("  Input table extracts value from read values\n");
+		} else {
+			printf("  Value extraction read values is explicit\n");
+		}
+		printf("  Input value re-ordering is [");
+		for (e = 0; e < g->id; e++)
+			printf("%s%d",e > 0 ? " " : "", t->it_map[e]);
+		printf("]\n");
+		printf("  Input table entry size = %d bytes\n",t->it_ts);
+		if (t->it_xs) {
+			printf("  Separate Interp. and Simplex index values\n");
+			printf("  Interp. index is at offset %d, size %d bytes\n",t->ix_eo, t->ix_es);
+			printf("  Simplex index is at offset %d, size %d bytes\n",t->sx_eo, t->sx_es);
+		} else {
+			printf("  Combined Interp. and Simplex index values\n");
+			printf("  Values are size %d bytes\n",t->it_ts);
+			printf("  Interp. index = %d bits\n",t->ix_ab);
+			printf("  Simplex index = %d bits\n",t->sx_ab);
+		}
+#endif /* VERBOSE */
+	}
+
+	/* Figure out output table stuff */
+	{
+		/* A variable to hold the index into an output table */
+		if ((f->otit = findord(f, g->prec)) < 0) {
+			fprintf(stderr,"Can't find output table index size\n");
+			exit(-1);
+		}
+		f->otit = nord(f,f->otit);				/* Make temp variable natural size */
+
+		if (g->out.pint != 0)	/* Pixel interleaved */
+			f->nop = 1;
+		else
+			f->nop = g->od;
+	
+		/* Figure out the output pointer types */
+		f->otvt = 0;			/* Output table value type */
+		for (e = 0; e < f->nop; e++) {
+			if ((f->opt[e] = findord(f, g->out.bpch[e])) < 0) {
+				fprintf(stderr,"Output channel size can't be handled\n");
+				exit(-1);
+			}
+			if (f->opt[e] > f->otvt)
+				f->otvt = f->opt[e];	/* Make value type big enough for any channel size */
+		}
+		t->ot_ts = a->ords[f->otvt].bits/8;	/* Output table entry size in bytes */
+
+		/* Setup information on data placement in output table enries */
+		for (e = 0; e < g->od; e++) {
+			t->ot_off[e] = g->out.bov[e];		/* Transfer info from generation spec. */
+			t->ot_bits[e] = g->out.bpv[e];
+		}
+	}
+
+#ifdef VERBOSE
+	/* Summarise the output table arrangements */
+	printf("  Output value re-ordering is [");
+	for (e = 0; e < g->od; e++)
+		printf("%s%d",e > 0 ? " " : "", t->im_map[e]);
+	printf("]\n");
+	printf("\n");
+
+	printf("Output table structure:\n");
+	printf("  Entry size = %d bytes\n",t->ot_ts);
+	printf("  Output value placement within each enry is:\n");
+	for (e = 0; e < f->nop; e++) {
+		printf("    %d: Offset %d bits, size %d bits\n", e, t->ot_off[e], t->ot_bits[e]);
+	}
+#endif /* VERBOSE */
+
+	/* Compute the maximum interpolation table resolution we will be able to handle */
+	{
+		int res, ores;
+
+		res = calc_res(g->id, t->ix_ab);
+		ores = calc_ores(g->id, t->vo_ab, t->vo_om);
+		f->ixmxres = res < ores ? res : ores;
+	}
+
+	/* Compute the maximum simplex table resolution we will be able to handle */
+	if (t->sort) {
+		f->sxmxres = 0;
+	} else {
+		f->sxmxres = calc_res(g->id, t->sx_ab);
+	}
+
+#ifdef VERBOSE
+	printf("Emitting introductory code\n"); fflush(stdout);
+#endif /* VERBOSE */
+
+	/* Start of code generation */
+	doheader(f);			/* Output the header comments */
+
+	/* We need an include file */
+	line(f,"#ifndef  IMDI_INCLUDED");
+	line(f,"#include <memory.h>");
+	line(f,"#include \"imdi_imp.h\"");
+	line(f,"#define  IMDI_INCLUDED");
+	line(f,"#endif  /* IMDI_INCLUDED */");
+	cr(f);
+
+	/* Declare our explicit pointer type */
+	line(f,"#ifndef DEFINED_pointer");
+	line(f,"#define DEFINED_pointer");
+	line(f,"typedef unsigned char * pointer;");
+	line(f,"#endif");
+	cr(f);
+
+	/* Declare our explicit structure access macros */
+
+#ifdef VERBOSE
+	printf("Declaring macros\n"); fflush(stdout);
+#endif /* VERBOSE */
+
+	/* Macros for accessing input table entries */
+	if (t->sort) {
+		if (t->it_xs) {
+			line(f,"/* Input table interp. index */");
+			line(f,"#define IT_IX(p, off) *((%s *)((p) + %d + (off) * %d))",
+			     a->ords[f->ixet].name, t->ix_eo, t->it_ts);
+			cr(f);
+			if (t->wo_xs) {
+				line(f,"/* Input table input weighting enty */");
+				line(f,"#define IT_WE(p, off) *((%s *)((p) + %d + (off) * %d))",
+				     a->ords[f->weet].name, t->we_eo, t->it_ts);
+				cr(f);
+				line(f,"/* Input table input offset value enty */");
+				line(f,"#define IT_VO(p, off) *((%s *)((p) + %d + (off) * %d))",
+				     a->ords[f->voet].name, t->vo_eo, t->it_ts);
+				cr(f);
+			} else {
+				line(f,"/* Input table input weighting/offset value enty */");
+				line(f,"#define IT_WO(p, off) *((%s *)((p) + %d + (off) * %d))",
+				     a->ords[f->woet].name, t->wo_eo, t->it_ts);
+				cr(f);
+			}
+		} else {
+			line(f,"/* Input table interp index, weighting and vertex offset */");
+			line(f,"#define IT_IT(p, off) *((%s *)((p) + %d + (off) * %d))",
+			     a->ords[f->itet].name, 0, t->it_ts);
+			cr(f);
+		}
+
+		/* Macro to conditionally exchange two varibles */
+		/* Doing this in place using an xor seems to be fastest */
+		/* on most architectures. */
+		line(f,"/* Conditional exchange for sorting */");
+		if (t->wo_xs) {
+			line(f,"#define CEX(A, AA, B, BB) if (A < B) { \\");
+			line(f,"            A ^= B; B ^= A; A ^= B; AA ^= BB; BB ^= AA; AA ^= BB; }");
+		} else
+			line(f,"#define CEX(A, B) if (A < B) { A ^= B; B ^= A; A ^= B; }");
+		cr(f);
+
+	} else {	/* Simplex table */
+		if (t->it_xs) {
+			line(f,"/* Input table interp. index */");
+			line(f,"#define IT_IX(p, off) *((%s *)((p) + %d + (off) * %d))",
+			     a->ords[f->ixet].name, t->ix_eo, t->it_ts);
+			cr(f);
+			line(f,"/* Input table simplex index enty */");
+			line(f,"#define IT_SX(p, off) *((%s *)((p) + %d + (off) * %d))",
+			     a->ords[f->sxet].name, t->sx_eo, t->it_ts);
+			cr(f);
+		} else {
+			line(f,"/* Input table inter & simplex indexes */");
+			line(f,"#define IT_IT(p, off) *((%s *)((p) + %d + (off) * %d))",
+			     a->ords[f->itet].name, 0, t->it_ts);
+			cr(f);
+		}
+	}
+
+	if (!t->sort) {
+		/* Macro for computing a simplex table entry */
+		line(f,"/* Simplex weighting table access */");
+		line(f,"#define SW_O(off) ((off) * %d)", t->sm_ts);
+		cr(f);
+
+		/* Macros for accessing the contents of the simplex table */
+		if (t->wo_xs) { 			/* If separate */
+			line(f,"/* Simplex table - get weighting value */");
+			line(f,"#define SX_WE(p, v) *((%s *)((p) + (v) * %d + %d))",
+			     a->ords[f->weet].name, t->wo_es, t->we_eo);
+			cr(f);
+
+			line(f,"/* Simplex table - get offset value */");
+			line(f,"#define SX_VO(p, v) *((%s *)((p) + (v) * %d + %d))",
+			     a->ords[f->voet].name, t->wo_es, t->vo_eo);
+			cr(f);
+
+		} else {	/* Combined */
+			line(f,"/* Simplex table - get weighting/offset value */");
+			line(f,"#define SX_WO(p, v) *((%s *)((p) + (v) * %d))",
+			     a->ords[f->woet].name, t->wo_es);
+			cr(f);
+		}
+	}
+
+	/* Macro for computing an interpolation table entry */
+	line(f,"/* Interpolation multi-dim. table access */");
+	line(f,"#define IM_O(off) ((off) * %d)", t->im_ts);
+	cr(f);
+
+	/* Macro for accessing an entry in the interpolation table */
+	line(f,"/* Interpolation table - get vertex values */");
+
+	if (t->im_fn > 0) {
+		/* Arguments to macro are cell base address, vertex offset, data offset */
+
+		if (f->imfvt == f->iafvt) {	/* Table and accumulator are the same size */
+			if (!timp || t->im_fn == 1) 
+				line(f,"#define IM_FE(p, v, c) *((%s *)((p) + (v) * %d + (c) * %d))",
+				     a->ords[f->imfvt].name, t->im_oc, t->im_fs);
+			else {
+				line(f,"#define IM_TP(p, v) ((p) + (v) * %d)", t->im_oc);
+				line(f,"#define IM_FE(p, c) *((%s *)((p) + (c) * %d))",
+				     a->ords[f->imfvt].name, t->im_fs);
+			}
+		} else {					/* Expand single table entry to accumulator size */
+			if (!timp || t->im_fn == 1) 
+				line(f,"#define IM_FE(p, v, c) ((%s)*((%s *)((p) + (v) * %d + (c) * %d)))",
+				     a->ords[f->iafvt].name,
+				     a->ords[f->imfvt].name, t->im_oc, t->im_fs);
+			else {
+				line(f,"#define IM_TP(p, v) ((p) + (v) * %d)", t->im_oc);
+				line(f,"#define IM_FE(p, c) ((%s)*((%s *)((p) + (c) * %d)))",
+				     a->ords[f->iafvt].name,
+				     a->ords[f->imfvt].name, t->im_fs);
+			}
+		}
+	}
+	if (t->im_pn > 0) {
+		/* Arguments to macro are cell base address, vertex offset */
+		/* There is no data offset since there can be only be one partial entry */
+
+		if (f->imfvt == f->iafvt)	/* Table and accumulator are the same size */
+			line(f,"#define IM_PE(p, v) *((%s *)((p) + %d + (v) * %d))",
+			     a->ords[f->impvt].name, t->im_fn * t->im_fs, t->im_oc);
+		else						/* Expand single table entry to accumulator size */
+			line(f,"#define IM_PE(p, v) ((%s)*((%s *)((p) + %d + (v) * %d)))",
+			     a->ords[f->iafvt].name,
+			     a->ords[f->impvt].name, t->im_fn * t->im_fs, t->im_oc);
+	}
+	cr(f);
+
+	/* Macro for accessing an output table entry */
+	line(f,"/* Output table indexes */");
+	line(f,"#define OT_E(p, off) *((%s *)((p) + (off) * %d))",
+	     a->ords[f->otvt].name, t->ot_ts);
+	cr(f);
+
+	/* =============================================== */
+
+#ifdef VERBOSE
+	printf("Starting interpolation function\n"); fflush(stdout);
+#endif /* VERBOSE */
+
+	/* Declare the function */
+	line(f,"void");
+	line(f, "imdi_k%d(",index);
+	line(f, "imdi *s,			/* imdi context */");
+	line(f, "void **outp,		/* pointer to output pointers */");
+	line(f, "void **inp,		/* pointer to input pointers */");
+	line(f, "unsigned int npix	/* Number of pixels to process */");
+	line(f, ") {");
+	inc(f);
+
+	/* We need access to the imdi_imp */
+	line(f, "imdi_imp *p = (imdi_imp *)(s->impl);");
+
+	/* Declare the input pointers and init them */
+	for (e = 0; e < f->nip; e++) {
+		line(f, "%s *ip%d = (%s *)inp[%d];",
+		     a->ords[f->ipt[e]].name, e, a->ords[f->ipt[e]].name, e);
+	}
+
+	/* Declare the output pointers and init them */
+	for (e = 0; e < f->nop; e++) {
+		line(f, "%s *op%d = (%s *)outp[%d];",
+		     a->ords[f->opt[e]].name, e, a->ords[f->opt[e]].name, e);
+	}
+
+	/* Declare and intialise the end pointer */
+	line(f, "%s *ep = ip0 + npix * %d ;",
+	        a->ords[f->ipt[0]].name, g->in.chi[0]);
+
+	/* Declare and initialise the input table pointers */
+	for (e = 0; e < g->id; e++)
+		line(f,"pointer it%d = (pointer)p->in_tables[%d];",e,e);
+
+	/* Declare and initialise the output table pointers */
+	for (e = 0; e < g->od; e++)
+		line(f,"pointer ot%d = (pointer)p->out_tables[%d];",e,e);
+
+	if (!t->sort) {
+		/* Declare and initialise the Simplex weighting base pointer */
+		line(f,"pointer sw_base = (pointer)p->sw_table;");
+	}
+
+	/* Declare and initialise the Interpolation multidim base pointer */
+	line(f,"pointer im_base = (pointer)p->im_table;");
+
+	/* Figure out whether input channel reads can be used directly as table offsets */
+	t->it_ix = 1;				/* Default use input table lookup to extract value */
+
+	if (g->in.packed != 0)
+		t->it_ix = 0;				/* Extract will be done explicitly */
+
+	for (e = 0; e < g->id; e++) {
+		int ee = (g->in.pint != 0) ? 0 : e;		/* bpch index */
+
+		if ((g->in.bov[e] + g->in.bpv[e]) <= 12)
+			continue;							/* Table can do extract */
+
+		if (g->in.bov[e] != 0 || g->in.bpv[e] != g->in.bpch[ee]) {
+			t->it_ix = 0;						/* Extract will be done explicitly */
+			break;
+		}
+	}
+	
+	/* ------------------------------- */
+#ifdef VERBOSE
+	printf("Starting pixel processing loop\n"); fflush(stdout);
+#endif /* VERBOSE */
+
+	/* Start the pixel processing loop */
+	cr(f);
+	sline(f, "for(;ip0 < ep;");
+	for (e = 0; e < f->nip; e++)
+		mline(f, " ip%d += %d,", e, g->in.chi[e]);
+	for (e = 0; e < f->nop; e++)
+		mline(f, " op%d += %d%s", e, g->out.chi[e], ((e+1) < f->nop) ? "," : "");
+	eline(f, ") {");
+	inc(f);
+
+	/* Declare output value accumulator(s) */
+	for (i = 0; i < t->im_fn; i++) {
+		line(f,"%s ova%d;	/* Output value accumulator */",a->ords[f->iafvt].name,i);
+	}
+	for (; i < f->ian; i++) {
+		line(f,"%s ova%d;	/* Output value partial accumulator */",a->ords[f->iapvt].name,i);
+	}
+
+	/* Context around interp/Simplex table lookup */
+	line(f, "{");
+	inc(f);
+
+	if (!t->sort)
+		line(f,"pointer swp;");		/* Declare Simplex weighting pointer */
+	line(f,"pointer imp;");			/* Declare Interpolation multidim pointer */
+
+	/* Declare the input weighting/vertex offset variables */
+	if (t->sort) {
+		for (e = 0; e < g->id; e++) {
+			if (t->wo_xs) {
+				line(f,"%s we%d;	/* Weighting value variable */",
+				       a->ords[f->wevt].name, e);
+				line(f,"%s vo%d;	/* Vertex offset variable */",
+				       a->ords[f->vovt].name, e);
+			} else {
+				line(f,"%s wo%d;	/* Weighting value and vertex offset variable */",
+				       a->ords[f->wovt].name, e);
+			}
+		}
+	}
+
+	/* Context around input table processing */
+	line(f, "{");
+	inc(f);
+
+	/* Declare the table index variables/input weighting/vertex offset variables */
+	if (t->sort) {
+		if (!t->it_xs)
+			line(f,"%s ti;		/* Input table entry variable */",a->ords[f->itvt].name);
+		line(f,"%s ti_i;	/* Interpolation index variable */",a->ords[f->ixvt].name);
+	} else {
+		if (t->it_xs) {
+			line(f,"%s ti_s;	/* Simplex index variable */",a->ords[f->sxvt].name);
+			line(f,"%s ti_i;	/* Interpolation index variable */",a->ords[f->ixvt].name);
+		} else {
+			line(f,"%s ti;	/* Simplex+Interpolation index variable */",a->ords[f->itvt].name);
+		}
+	}
+
+	if (g->in.packed != 0)	/* We need to unpack from a single read */
+		line(f,"%s rdv;		/* Read value */",a->ords[f->ipt[0]].name);
+
+	if (t->it_ix == 0) {
+		int bv = 0;
+		for (e = 0; e < f->nip; e++) {	/* Find largest input type */
+			if (f->ipt[e] > bv)
+				bv = f->ipt[e];
+		}
+		bv = nord(f, bv);
+		line(f,"%s chv;	/* Channel value */",a->ords[bv].name);
+		f->chv_bits = a->ords[bv].bits;
+	}
+	cr(f);
+
+#ifdef VERBOSE
+	printf("Read code\n"); fflush(stdout);
+#endif /* VERBOSE */
+
+	/* For all the input channels */
+	for (e = 0; e < g->id; e++) {
+		char rde[50];		/* Read expression */
+		char toff[50];		/* Table offset expression */
+		int ee = (g->in.pint != 0) ? 0 : e;		/* bpch index */
+
+		if (g->in.pint != 0) 	/* Pixel interleaved */
+			sprintf(rde,"ip0[%d]",e);	/* Offset from single pointer */
+		else
+			sprintf(rde,"*ip%d",e);		/* Pointer per channel */
+
+		if (g->in.packed != 0) {
+			if (e == 0)
+				line(f,"rdv = %s;",rde);	/* Do single read */
+			sprintf(rde,"rdv");				/* Use read value for extraction */
+		}
+
+		if (t->it_ix == 0) {
+			if (g->in.bov[e] == 0 ) {				/* No offset */
+				if (g->in.bpv[e] == g->in.bpch[ee])	/* No mask */
+					line(f,"chv = %s;",rde);
+				else								/* Just mask  */
+					line(f,"chv = (%s & %s);",rde, hmask(g->in.bpv[e]));
+			} else {								/* Offset */
+				if ((g->in.bov[e] + g->in.bpv[e]) == g->in.bpch[ee])
+					line(f,"chv = (%s >> %d);",rde, g->in.bov[e]);
+				else {								/* Offset and mask */
+					if (a->shfm || g->in.bpv[e] > 32) {
+						/* Extract using just shifts */
+						line(f,"chv = ((%s << %d) >> %d);", rde,
+						        f->chv_bits - g->in.bpv[e] - g->in.bov[e],
+						        f->chv_bits - g->in.bpv[e]);
+					} else {
+						/* Extract using shift and mask */
+						line(f,"chv = ((%s >> %d) & %s);",
+						        rde, g->in.bov[e], hmask(g->in.bpv[e]));
+					}
+				}
+			}
+			sprintf(toff,"chv");
+		} else {									/* No extraction */
+			sprintf(toff,"%s",rde);
+		}
+
+		if (t->sort) {
+			if (t->it_xs) {
+				line(f,"ti_i %s= IT_IX(it%d, %s);", e ? "+" : " ", e, toff);
+				if (t->wo_xs) {
+					line(f,"we%d   = IT_WE(it%d, %s);", e, e, toff);
+					line(f,"vo%d   = IT_VO(it%d, %s);", e, e, toff);
+				} else {
+					line(f,"wo%d   = IT_WO(it%d, %s);", e, e, toff);
+				}
+			} else {	/* All three combined */
+				line(f,"ti = IT_IT(it%d, %s);", e, toff);
+				if (a->shfm || t->wo_ab > 32) {
+					/* Extract using just shifts */
+					line(f,"wo%d   = ((ti << %d) >> %d);	"
+					     "/* Extract weighting/vertex offset value */",
+					     e, a->ords[f->wovt].bits - t->wo_ab, a->ords[f->wovt].bits - t->wo_ab);
+					line(f,"ti_i %s= (ti >> %d);	"
+					     "/* Extract interpolation table value */",
+					     e ? "+" : " ", t->wo_ab);
+				} else {
+					/* Extract using shift and mask */
+					line(f,"wo%d   = (ti & %s);	"
+					     "/* Extract weighting/vertex offset value */",
+					     e, hmask(t->wo_ab));
+					line(f,"ti_i %s= (ti >> %d);	"
+					     "/* Extract interpolation table value */",
+					     e ? "+" : " ", t->wo_ab);
+				}
+			}
+
+		} else {	/* Simplex */
+			if (t->it_xs) {
+				/* ~~~~ should toff be forced to be a temp variable ?? */
+				/* (ie. force use of rde (above) if t->it_xs is nonz) */ 
+				line(f,"ti_i %s= IT_IX(it%d, %s);", e ? "+" : " ", e, toff);
+				line(f,"ti_s %s= IT_SX(it%d, %s);", e ? "+" : " ", e, toff);
+			} else {
+				line(f,"ti %s= IT_IT(it%d, %s);", e ? "+" : " ", e, toff);
+			}
+		}
+	}
+
+#ifdef VERBOSE
+	printf("Index extraction code\n"); fflush(stdout);
+#endif /* VERBOSE */
+
+	cr(f);
+
+	if (t->sort) {
+		/* Extract Simplex and Interpolation indexes from accumulator */
+		line(f,"imp = im_base + IM_O(ti_i);		/* Compute interp. table entry pointer */");
+	} else {
+		if (t->it_xs) {		/* Extract Simplex and Interpolation indexes from accumulator */
+			line(f,"swp = sw_base + SW_O(ti_s);		/* Compute simplex table entry pointer */");
+			line(f,"imp = im_base + IM_O(ti_i);		/* Compute interp. table entry pointer */");
+		} else {
+			line(f,"imp = im_base + IM_O(ti >> %d);		"
+			     "/* Extract interp. index and comp. entry */",
+			     t->sx_ab);
+			if (a->shfm || t->sx_ab > 32) {
+				/* Extract using just shifts */
+				line(f,"swp = sw_base + SW_O((ti << %d) >> %d);	"
+				     "/* Extract simplex index & comp. entry */",
+				     a->ords[f->itvt].bits - t->sx_ab, a->ords[f->itvt].bits - t->sx_ab);
+			} else {
+				/* Extract using shift and mask */
+				line(f,"swp = sw_base + SW_O(ti & %s);	"
+				     "/* Extract simplex index and comp. entry */",
+				     hmask(t->sx_ab));
+			}
+		}
+	}
+
+	/* Do the explicit sort now */
+	if (t->sort) {
+		cr(f);
+		/* Sort from largest to smallest using a selection sort */
+		/* Use simple sequence for the moment. More elaborate sequence */
+		/* may allow other optimisations. */
+		line(f,"/* Sort weighting values and vertex offset values */");
+		for (i = 0; i < (g->id-1); i++) {
+			for (e = i+1; e < g->id; e++) {
+				if (t->wo_xs)
+					line(f,"CEX(we%d, vo%d, we%d, vo%d);",i,i,e,e);
+				else
+					line(f,"CEX(wo%d, wo%d);",i,e);
+			}
+		}
+	}
+
+	/* End of input table processing context */
+	dec(f);
+	line(f,"}");
+
+	line(f,"{");	/* Context around vertex lookup and accumulation */
+	inc(f);
+
+	/* Declare vertex offset and weight variables */
+	if (t->sort && t->wo_xs == 0) {
+		line(f,"%s nvof;	/* Next vertex offset value */",a->ords[f->vovt].name);
+	} else {
+		if (!t->wo_xs)	/* If combined in table */
+			line(f,"%s vowr;	/* Vertex offset/weight value */",a->ords[f->wovt].name);
+	}
+	line(f,"%s vof;	/* Vertex offset value */",a->ords[f->vovt].name);
+	line(f,"%s vwe;	/* Vertex weighting */",a->ords[f->wevt].name);
+	if (timp && t->im_fn > 1) 
+		line(f,"pointer timp;		/* Temporary interpolation table pointer */");
+	cr(f);
+
+#ifdef VERBOSE
+	printf("Vertex offset and weight code\n"); fflush(stdout);
+#endif /* VERBOSE */
+
+	/* For each vertex in the simplex */
+	for (e = 0; e < (g->id +1); e++) {
+
+		if (t->sort) {
+
+			if (e == 0) {
+				line(f,"vof = 0;				/* First vertex offset is 0 */");
+			} else {
+				if (t->wo_xs)
+					line(f,"vof += vo%d;			/* Move to next vertex */",e-1);
+				else
+					line(f,"vof += nvof;			/* Move to next vertex */");
+			}
+
+			/* Extract the vertex offset and weight values from the sorted input values */
+			if (e < g->id && !t->wo_xs) {
+				if (a->shfm || t->vo_ab > 32) {
+					/* Extract using just shifts */
+					line(f,"nvof = ((wo%d << %d) >> %d);	"
+					     "/* Extract offset value */",
+					     e, a->ords[f->vovt].bits - t->vo_ab, a->ords[f->vovt].bits - t->vo_ab);
+					line(f,"wo%d = (wo%d >> %d);	"
+					     "	/* Extract weighting value */",
+					     e, e, t->vo_ab);
+				} else {
+					/* Extract using shift and mask */
+					line(f,"nvof = (wo%d & %s);	"
+					     "/* Extract offset value */",
+					     e, hmask(t->vo_ab));
+					line(f,"wo%d = (wo%d >> %d);	"
+					     "	/* Extract weighting value */",
+					     e, e, t->vo_ab);
+				}
+			}
+			/* Compute the weighting value */
+			if (!t->wo_xs) {
+				if (e == 0) {
+					line(f,"vwe = %d - wo%d;		/* Baricentric weighting */", 1 << g->prec, e);
+				} else if (e < g->id) {
+					line(f,"vwe = wo%d - wo%d;		/* Baricentric weighting */", e-1, e);
+				} else {
+					line(f,"vwe = wo%d;				/* Baricentric weighting */", e-1);
+				}
+			} else {
+				if (e == 0) {
+					line(f,"vwe = %d - we%d;		/* Baricentric weighting */", 1 << g->prec, e);
+				} else if (e < g->id) {
+					line(f,"vwe = we%d - we%d;		/* Baricentric weighting */", e-1, e);
+				} else {
+					line(f,"vwe = we%d;				/* Baricentric weighting */", e-1);
+				}
+			}
+
+		} else {	/* Not sort */
+			/* Read the vertex offset and weight values from the simplex table */
+			if (t->wo_xs) { 			/* If separate */
+				line(f,"vof = SX_VO(swp, %d);	/* Read vertex offset value */", e);
+				line(f,"vwe = SX_WE(swp, %d);	/* Read vertex weighting value */", e);
+			} else { 			/* If combined in table */
+				line(f,"vowr = SX_WO(swp, %d);	/* Read vertex offset+weighting values */", e);
+				if (a->shfm || t->vo_ab > 32) {
+					/* Extract using just shifts */
+					line(f,"vof = ((vowr << %d) >> %d);	"
+					     "/* Extract offset value */",
+					     a->ords[f->vovt].bits - t->vo_ab, a->ords[f->vovt].bits - t->vo_ab);
+					line(f,"vwe = (vowr >> %d);	"
+					     "/* Extract weighting value */",
+					     t->vo_ab);
+				} else {
+					/* Extract using shift and mask */
+					line(f,"vof = (vowr & %s);	"
+					     "/* Extract offset value */",
+					     hmask(t->vo_ab));
+					line(f,"vwe = (vowr >> %d);	"
+					     "/* Extract weighting value */",
+					     t->vo_ab);
+				}
+			}
+		}
+
+		/* Lookup the vertex value, weight it, and accumulate it into output value */
+		if (timp && t->im_fn > 1) 
+			line(f,"timp = IM_TP(imp, vof);	/* Vertex address */");
+		for (i = 0; i < f->ian; i++) {		/* For each output accumulation chunk */
+			if (i < t->im_fn) { 	/* Full entry */
+				if (!timp || t->im_fn == 1) 
+					line(f,"ova%d %s= IM_FE(imp, vof, %d) * vwe;	"
+					     "/* Accumulate weighted output values */",
+					     i, e ? "+" : " ", i);
+				else
+					line(f,"ova%d %s= IM_FE(timp, %d) * vwe;	"
+					     "/* Accumulate weighted output values */",
+					     i, e ? "+" : " ", i);
+			} else				/* One partial entry */
+				line(f,"ova%d %s= IM_PE(imp, vof) * vwe;	"
+				     "/* Accumulate last weighted output values */",
+				     i, e ? "+" : " ");
+		}
+	}
+
+	dec(f);
+	line(f, "}"); 	/* End of output value lookup context */
+
+	dec(f);
+	line(f, "}"); 	/* End of output value accumulation context */
+
+	/* Start of output lookup and write */
+	line(f,"{");
+	inc(f);
+
+#ifdef VERBOSE
+	printf("Output table code\n"); fflush(stdout);
+#endif /* VERBOSE */
+
+	{
+		char wre[50];		/* Write destination expression */
+
+		if (g->out.packed != 0)	/* We need to pack results into a single write */
+			line(f,"%s wrv;		/* Write value */",a->ords[f->ipt[0]].name);
+	
+		/* Declare temporary to hold index into output lookup table */
+		line(f,"%s oti;	/* Vertex offset value */",a->ords[f->otit].name);
+	
+		/* For each accumulator value */
+		/* (Assume they are in output order for the moment ?) */
+		for (e = i = 0; i < f->ian; i++) {		/* For each output accumulation chunk */
+			int vpa = i < t->im_fn ? t->im_fv : t->im_pv;		/* Chanel values per accumulator */
+			int oat = i < t->im_fn ? f->iafvt : f->iapvt;		/* Output accumulator type */
+			int ee;		/* Relative e to this accumulator */
+	
+			/* For each output value in this accumulator */
+			for (ee = 0; ee < vpa && e < g->od; ee++, e++) {
+				int off, size;		/* Bits to be extracted */
+		
+				/* Extract wanted 8 bits from the 8.8 bit result in accumulator */
+				off = ee * f->iaovb + (f->iaovb - g->prec);	
+				size = g->prec;
+	
+				if (e == 0 || g->out.packed == 0) {
+					if (g->out.pint != 0) 	/* Pixel interleaved */
+						sprintf(wre,"op0[%d]",e);	/* Offset from single pointer */
+					else
+						sprintf(wre,"*op%d",e);		/* Pointer per channel */
+				}
+	
+				if (a->shfm || size > 32) {
+					/* Extract using just shifts */
+					line(f,"oti = ((ova%d << %d) >> %d);	"
+					     "/* Extract integer part of result */",
+					     i, a->ords[oat].bits - off - size, a->ords[oat].bits - size);
+				} else {
+					/* Extract using shift and mask */
+					line(f,"oti = ((ova%d >> %d) & %s);	"
+					     "/* Extract integer part of result */",
+					     i, off, hmask(size));
+				}
+	
+				/* Lookup in output table and write to destination */
+				if (g->out.packed != 0) {
+					line(f,"wrv %s= OT_E(ot%d, oti);", e ? "+" : "", e);
+				} else {
+					line(f,"%s = OT_E(ot%d, oti);	/* Write result */", wre, e);
+				}
+			}
+		}
+	
+		if (g->out.packed != 0) {	/* Write out the accumulated value */
+			line(f,"%s = wrv;	/* Write result */", wre);
+		}
+	}
+
+	/* The end of the output lookup and write */
+	dec(f);
+	line(f, "}");
+
+	/* The end of the pixel processing loop */
+	dec(f);
+	line(f, "}");
+
+	/* The end of the function */
+	dec(f);
+	line(f, "}");
+
+	/* Undefine all the macros */
+	if (t->sort) {
+		if (t->it_xs) {
+			if (t->wo_xs) {
+				line(f,"#undef IT_WE");
+				line(f,"#undef IT_VO");
+			} else
+				line(f,"#undef IT_WO");
+			line(f,"#undef IT_IX");
+		} else {
+			line(f,"#undef IT_IT");
+		}
+		line(f,"#undef CEX");
+	} else {
+		if (t->it_xs) {
+			line(f,"#undef IT_IX");
+			line(f,"#undef IT_SX");
+		} else {
+			line(f,"#undef IT_IT");
+		}
+
+		line(f,"#undef SW_O");
+		if (t->wo_xs) {
+			line(f,"#undef SX_WE");
+			line(f,"#undef SX_VO");
+		} else {
+			line(f,"#undef SX_WO");
+		}
+	}
+	line(f,"#undef IM_O");
+	if (t->im_fn > 0) {
+		if (timp && t->im_fn > 1) 
+			line(f,"#undef IM_TP");
+		line(f,"#undef IM_FE");
+	}
+	if (t->im_pn > 0) {
+		line(f,"#undef IM_PE");
+	}
+	line(f,"#undef OT_E");
+
+	/* =============================================== */
+#ifdef VERBOSE
+	printf("Done interpolation code\n"); fflush(stdout);
+#endif /* VERBOSE */
+
+	/* =============================================== */
+
+	{
+		int sog = sizeof(genspec);	/* Size of genspec structure */
+		unsigned char *dp = (unsigned char *)g;
+
+		int s_stres, s_itres;	/* Save values */
+
+		s_stres = g->stres;
+		s_itres = g->itres;
+		g->stres = f->sxmxres;		/* Set maximum values */
+		g->itres = f->ixmxres;
+		
+		/* Declare the generation structure data function */
+		cr(f);
+		line(f,"void");
+		line(f, "imdi_k%d_gen(",index);
+		line(f, "genspec *g			/* structure to be initialised */");
+		line(f, ") {");
+		inc(f);
+	
+		/* Declare the genspec initialisation data */
+		line(f, "static unsigned char data[] = {");
+		inc(f);
+		for (i = 0; i < sog; i++) {
+			if ((i & 7) == 0)
+				sline(f,"");
+			mline(f, "0x%02x%s ", dp[i], (i+1) < sog ? "," : "", dp[i]);
+			if ((i & 7) == 7 || (i+1) == sog)
+				eline(f,"");
+		}
+		dec(f);
+		line(f, "};	/* Structure image */");
+
+		cr(f);
+		line(f, "memcpy(g, data, sizeof(data));	/* Initialise the structure */");
+		/* The end of the function */
+		dec(f);
+		line(f, "}");
+
+		g->stres = s_stres;		/* Restore entry values */
+		g->itres = s_itres;
+	}
+
+	/* =============================================== */
+
+	{
+		int sot = sizeof(tabspec);	/* Size of tabspec structure */
+		unsigned char *dp = (unsigned char *)t;
+
+		/* Declare the generation structure data function */
+		cr(f);
+		line(f,"void");
+		line(f, "imdi_k%d_tab(",index);
+		line(f, "tabspec *t			/* structure to be initialised */");
+		line(f, ") {");
+		inc(f);
+	
+		/* Declare the genspec initialisation data */
+		line(f, "static unsigned char data[] = {");
+		inc(f);
+		for (i = 0; i < sot; i++) {
+			if ((i & 7) == 0)
+				sline(f,"");
+			mline(f, "0x%02x%s ", dp[i], (i+1) < sot ? "," : "", dp[i]);
+			if ((i & 7) == 7 || (i+1) == sot)
+				eline(f,"");
+		}
+		dec(f);
+		line(f, "};	/* Structure image */");
+
+		cr(f);
+		line(f, "memcpy(t, data, sizeof(data));	/* Initialise the structure */");
+		/* The end of the function */
+		dec(f);
+		line(f, "}");
+	}
+
+	/* =============================================== */
+
+	cr(f); cr(f); cr(f); cr(f); cr(f); cr(f);
+
+	return 0;
+}
+
+
+/* Return bits needed to store index into table of */
+/* given resolution and dimensionality. */
+static int
+calc_bits(
+int dim,
+int res) {
+
+	return ceil(log((double)res) * (double)dim/log(2.0) - 1e-14);
+}
+
+/* Return maximum resolution possible given dimensionality */
+/* and number of index bits. */
+static int
+calc_res(
+int dim,
+int bits) {
+	double fres;
+
+	fres = log(2.0) * (double)bits/(double)dim;
+	if (fres > 12 || (fres = exp(fres)) > 65536.0)
+		fres = 65536.0;		/* Limit to a sane value */
+	return (int)(fres + 1e-14);
+}
+
+/* Return bits needed to store a relative offset of 1, */
+/* into a table of given resolution, dimensionality , and */
+/* entry size. */
+static int
+calc_obits(
+int dim,
+int res,
+int esize) {
+	double off;		/* Maximum diagonal offset value */
+	int bits;
+
+	if (res == 0 || res == 1)
+		return 0;
+	if (dim == 1)
+		off = esize;
+	else {
+		off = (double)esize * floor(exp(log((double)res) * dim - log(res-1.0)));
+	}
+
+	bits = ceil(log(off)/log(2.0) - 1e-14);
+	return bits;
+}
+
+/* Return maximum resolution possible given dimensionality */
+/* number of index bits, and entry size */
+static int
+calc_ores(
+int dim,
+int bits,
+int esize) {
+	int res;
+
+	/* Find resolution. Stop at arbitrary 65536 */
+	for (res = 1; res < 65537; res++) {
+		int bn;
+		bn = calc_obits(dim, res, esize);
+		if (bn > bits) {
+			return res-1;
+		}
+	}
+	return res-1;
+}
+
+
+
+/* Output the introductory comments */
+static void
+doheader(
+	fileo *f
+) {
+	genspec *g = f->g;
+	tabspec *t = f->t;
+	mach_arch *a  = f->a;
+	int e;
+
+	/* - - - - - - - - - - - - */
+	/* Output file title block */
+	line(f,"/* Integer Multi-Dimensional Interpolation */");
+	line(f,"/* Interpolation Kernel Code */");
+	line(f,"/* Generated by cgen */");
+	line(f,"/* Copyright 2000 - 2002 Graeme W. Gill */");
+	line(f,"/* This material is licenced under the GNU GENERAL PUBLIC LICENCE :- */\n");
+	line(f,"/* see the Licence.txt file for licencing details.*/\n");
+	cr(f);
+
+	/* - - - - - - - - - - - - */
+	/* Output the specification */
+	line(f,"/*");
+	line(f,"   Interpolation kernel specs:");
+	cr(f);
+	line(f,"   Input channels per pixel = %d",g->id);
+	for (e = 0; e < g->id; e++) {
+		line(f,"   Input channel %d bits = %d",e, g->in.bpch[e]);
+		line(f,"   Input channel %d increment = %d",e, g->in.chi[e]);
+	}
+	if (g->in.pint != 0)
+		line(f,"   Input is channel interleaved");
+	else
+		line(f,"   Input is plane interleaved");
+
+	if (g->in.packed != 0)
+		line(f,"   Input channels are packed into one word");
+	else
+		line(f,"   Input channels are separate words");
+
+	if (t->it_ix)
+		line(f,"   Input value extraction is done in input table lookup");
+	cr(f);
+	
+	line(f,"   Output channels per pixel = %d",g->od);
+	for (e = 0; e < g->od; e++) {
+		line(f,"   Output channel %d bits = %d",e, g->out.bpch[e]);
+		line(f,"   Output channel %d increment = %d",e, g->out.chi[e]);
+	}
+	if (g->out.pint != 0)
+		line(f,"   Output is channel interleaved");
+	else
+		line(f,"   Output is plane interleaved");
+	cr(f);
+	if (g->out.packed != 0)
+		line(f,"   Output channels are packed into one word");
+	else
+		line(f,"   Output channels are separate words");
+
+	
+	if (t->sort)
+		line(f,"   Weight+voffset bits       = %d",t->sx_ab);
+	else
+		line(f,"   Simplex table index bits       = %d",t->sx_ab);
+	line(f,"   Interpolation table index bits = %d",t->ix_ab);
+	if (!t->sort)
+		line(f,"   Simplex table max resolution = %d",f->sxmxres);
+	line(f,"   Interpolation table max resolution = %d",f->ixmxres);
+	line(f," */");
+	cr(f);
+
+	/* - - - - - - - - - - - - */
+	line(f,"/*");
+	line(f,"   Machine architecture specs:");
+	cr(f);
+	if (a->bigend != 0)
+		line(f,"   Big Endian");
+	else
+		line(f,"   Little endian");
+
+	if (a->uwa != 0)
+		line(f,"   Using maximum sized memory accesses where possible");
+	else
+		line(f,"   Reading and writing pixel values separately");
+
+	line(f,"   Pointer size = %d bits",a->pbits);
+	cr(f);
+
+	for (e = 0; e < a->nords; e++) {
+		line(f,"   Ordinal size %2d bits is known as '%s'",
+		        a->ords[e].bits,a->ords[e].name);
+	}
+	line(f,"   Natural ordinal is '%s'", a->ords[a->natord].name);
+	cr(f);
+
+	for (e = 0; e < a->nints; e++) {
+		line(f,"   Integer size %2d bits is known as '%s'",
+		        a->ints[e].bits,a->ints[e].name);
+	}
+	line(f,"   Natural integer is '%s'", a->ints[a->natint].name);
+	cr(f);
+
+	line(f," */");
+	cr(f);
+}
+
+
+/* ---------------------------------------- */
+/* Architecture support */
+/* Find an ordinal with at least bits size */
+/* Return -1 if failed */
+int findord(
+fileo *f,
+int bits
+) {
+	mach_arch *a  = f->a;
+	int i;
+
+	for (i = 0; i < a->nords; i++) {
+		if (a->ords[i].bits >= bits)
+			return i;
+	}
+	return -1;
+}
+
+/* Round ordinal type up to natural size */
+int nord(
+	fileo *f,
+	int ov
+) {
+	if (ov >= 0 && ov < f->a->natord)
+		ov = f->a->natord;
+	return ov;
+}
+
+/* Find an ordinal with at least bits size, */
+/* or natural size, whichever is greater. */
+/* Return -1 if failed */
+int findnord(
+	fileo *f,
+	int bits
+) {
+	int ov;
+
+	ov = findord(f, bits);
+	ov = nord(f, ov);
+	return ov;
+}
+
+/* Find an integer with at least bits size */
+/* Return -1 if failed */
+int findint(
+	fileo *f,
+	int bits
+) {
+	mach_arch *a  = f->a;
+	int i;
+
+	for (i = 0; i < a->nints; i++) {
+		if (a->ints[i].bits >= bits)
+			return i;
+	}
+	return -1;
+}
+
+/* Round integer type up to natural size */
+int nint(
+	fileo *f,
+	int iv
+) {
+	if (iv >= 0 && iv < f->a->natint)
+		iv = f->a->natint;
+	return iv;
+}
+
+/* Find an interger with at least bits size, */
+/* or natural size, whichever is greater. */
+/* Return -1 if failed */
+int findnint(
+	fileo *f,
+	int bits
+) {
+	int iv;
+
+	iv = findint(f, bits);
+	iv = nint(f, iv);
+	return iv;
+}
+
+
+/* ------------------------------------ */
+/* File output support */
+
+/* Output a line to the file (including trailing \n) */
+void
+line(fileo *f, char *fmt, ...)
+{
+	int i;
+	va_list args;
+
+	/* Indent to the correct level */
+	for (i = 0; i < f->indt; i++)
+		fprintf(f->of,"	");
+	
+	va_start(args, fmt);
+	vfprintf(f->of, fmt, args);
+	va_end(args);
+	fprintf(f->of, "\n");
+}
+
+/* Output the start of a line to the file) */
+void
+sline(fileo *f, char *fmt, ...)
+{
+	int i;
+	va_list args;
+
+	/* Indent to the correct level */
+	for (i = 0; i < f->indt; i++)
+		fprintf(f->of,"	");
+	
+	va_start(args, fmt);
+	vfprintf(f->of, fmt, args);
+	va_end(args);
+}
+
+/* Output the middle of a line to the file) */
+void
+mline(fileo *f, char *fmt, ...)
+{
+	int i;
+	va_list args;
+
+	va_start(args, fmt);
+	vfprintf(f->of, fmt, args);
+	va_end(args);
+}
+
+/* Output the end of a line to the file (including trailing \n) */
+void
+eline(fileo *f, char *fmt, ...)
+{
+	int i;
+	va_list args;
+
+	va_start(args, fmt);
+	vfprintf(f->of, fmt, args);
+	va_end(args);
+	fprintf(f->of, "\n");
+}
+/* ------------------------------------ */
+
+
+
+

Added: trunk/gs/imdi/config.h
===================================================================
--- trunk/gs/imdi/config.h	2006-09-11 07:02:18 UTC (rev 7033)
+++ trunk/gs/imdi/config.h	2006-09-11 20:26:01 UTC (rev 7034)
@@ -0,0 +1,16 @@
+
+#ifndef __CONFIG_H__
+#define __CONFIG_H__
+
+/* General project wide configuration */
+
+
+/* Version of Argyll release */
+
+#define ARGYLL_VERSION 0x000503
+#define ARGYLL_VERSION_STR "0.53"
+
+/* Maximum file path length */
+#define MAXNAMEL 512
+
+#endif /* __CONFIG_H__ */

Added: trunk/gs/imd