summaryrefslogtreecommitdiff
path: root/scanin/scanrd.c
diff options
context:
space:
mode:
authorJörg Frings-Fürst <debian@jff-webhosting.net>2014-09-01 13:56:46 +0200
committerJörg Frings-Fürst <debian@jff-webhosting.net>2014-09-01 13:56:46 +0200
commit22f703cab05b7cd368f4de9e03991b7664dc5022 (patch)
tree6f4d50beaa42328e24b1c6b56b6ec059e4ef21a5 /scanin/scanrd.c
Initial import of argyll version 1.5.1-8debian/1.5.1-8
Diffstat (limited to 'scanin/scanrd.c')
-rw-r--r--scanin/scanrd.c4659
1 files changed, 4659 insertions, 0 deletions
diff --git a/scanin/scanrd.c b/scanin/scanrd.c
new file mode 100644
index 0000000..a915e44
--- /dev/null
+++ b/scanin/scanrd.c
@@ -0,0 +1,4659 @@
+
+/*
+ * Raster Color Target Scan Input module
+ * This is the core chart recognition code.
+ *
+ * Author: Graeme Gill
+ *
+ * Copyright 1995 - 2008 Graeme W. Gill, All right reserved.
+ * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
+ * see the License.txt file for licencing details.
+ */
+
+/*
+ * To Do:
+ * Add option to output a raster file made from the .cht and example values.
+ *
+ * Fix sboxes parameters/digitization to fix "droop" in box areas.
+ * Scale parameters with image size.
+ * To handle high res, introduce automatic sub-sampler.
+ * Change reference parser to make it more forgiving - use cgats parser ?
+ */
+
+#undef DEBUG
+
+#define VERSION "1.0"
+
+/* Behaviour defines */
+#undef DIAGN /* Allow diagonal connectivity of groups */
+#define AA_LINES /* Plot diagnostics using anti-aliased lines */
+
+#define MATCHCC 0.3 /* Match correlation threshold - reject any match under this */
+ /* (Might want to be able to override this in command line) */
+
+#define ALT_ROT_TH 0.7 /* Correlation threshold of alternate rotations to be greater than this */
+
+#define TH (20.0 * 20.0) /* Initial color change threshhold */
+
+#undef DBG
+#define dbgo stdout
+#define DBG(aaa) fprintf aaa, fflush(dbgo)
+
+#include <stdio.h>
+/* #include <fcntl.h> */ /* In case DOS binary stuff is needed */
+#include <string.h>
+#include <math.h>
+
+#include <stdlib.h>
+#include <sys/stat.h>
+/* #include <fname.h> */
+
+#include "numlib.h"
+#include "scanrd_.h"
+
+/* ------------------------------------------------- */
+/* Implementations of public functions */
+static void free_scanrd(scanrd *s);
+static int scanrd_reset(scanrd *s);
+static int scanrd_read(scanrd *ps, char *id, double *P, double *mP,
+ double *sdP, int *cnt);
+static unsigned int scanrd_error(scanrd *s, char **errm);
+
+/* Forward internal function declaration */
+static scanrd_ *new_scanrd(int flags, int verb, double gammav,
+ int (*write_line)(void *ddata, int y, char *src), void *ddata,
+ int w, int h, int d, int td, int p,
+ int (*read_line)(void *fdata, int y, char *dst), void *fdata,
+ char *refname);
+static int read_input(scanrd_ *s);
+static int calc_lines(scanrd_ *s);
+static int show_lines(scanrd_ *s);
+static int calc_perspective(scanrd_ *s);
+static int calc_rotation(scanrd_ *s);
+static int calc_elists(scanrd_ *s, int ref);
+static int write_elists(scanrd_ *s);
+static int read_relists(scanrd_ *s);
+static int do_match(scanrd_ *s);
+static int compute_ptrans(scanrd_ *s);
+static int compute_man_ptrans(scanrd_ *s, double *sfids);
+static int improve_match(scanrd_ *s);
+static int setup_sboxes(scanrd_ *s);
+static int do_value_scan(scanrd_ *s);
+static int compute_xcc(scanrd_ *s);
+//static int restore_best(scanrd_ *s);
+static int show_sbox(scanrd_ *s);
+static int show_groups(scanrd_ *s);
+static int scanrd_write_diag(scanrd_ *s);
+static void toRGB(unsigned char *dst, unsigned char *src, int depth, int bpp);
+static void XYZ2Lab(double *out, double *in);
+static void pval2Lab(double *out, double *in, int depth);
+/* ------------------------------------------------- */
+
+/* Read in a chart, and either create a reference or make values available, */
+/* by using reset() and read() to get values read */
+scanrd *do_scanrd(
+int flags, /* option flags */
+int verb, /* verbosity level */
+
+double gammav, /* Apprimate gamma encoding of image (0.0 = default 2.2) */
+double *sfid, /* Specified four fiducials x1, y1 .. x4, y4, NULL if auto recognition */
+ /* Typical clockwise from top left */
+
+int w, int h, /* Width and Height of input raster in pixels */
+int d, int td, int p, /* Useful plane depth, Total depth, Bit presision of input pixels */
+int (*read_line)(void *fdata, int y, char *dst), /* Read RGB line of source file */
+void *fdata, /* Opaque data for read_line */
+
+char *refname, /* reference file name */
+
+int (*write_line)(void *ddata, int y, char *src), /* Write RGB line of diag file */
+void *ddata /* Opaque data for write_line */
+) {
+ scanrd_ *s;
+
+ /* allocate the basic object */
+ if (verb >= 2)
+ DBG((dbgo,"About to allocate scanrd_ object\n"));
+ if ((s = new_scanrd(flags, verb, gammav, write_line, ddata, w, h, d, td, p, read_line, fdata, refname)) == NULL)
+ return NULL;
+
+ if (s->errv != 0) /* Some other error from new_scanrd() */
+ return (scanrd *)s;
+
+ if (s->verb >= 2)
+ DBG((dbgo,"About to read input tiff file and discover groups\n"));
+ if (read_input(s))
+ goto sierr; /* Error */
+
+ if (s->flags & SI_SHOW_GROUPS)
+ if (show_groups(s))
+ goto sierr; /* Error */
+
+ if (s->verb >= 2)
+ DBG((dbgo,"About to calculate edge lines\n"));
+ if (calc_lines(s))
+ goto sierr; /* Error */
+ if (s->verb >= 2)
+ DBG((dbgo,"%d useful edges out of %d\n",s->novlines, s->noslines));
+
+ if (s->flags & SI_PERSPECTIVE) {
+ if (s->verb >= 2)
+ DBG((dbgo,"About to calculate perspective correction\n"));
+ if (calc_perspective(s)) {
+ if (s->flags & SI_SHOW_LINES) {
+ s->flags &= ~SI_SHOW_PERS; /* Calc perspective failed! */
+ s->flags &= ~SI_SHOW_ROT; /* Calc rotation not done! */
+ show_lines(s);
+ }
+ goto sierr; /* Error */
+ }
+ }
+
+ if (s->verb >= 2)
+ DBG((dbgo,"About to calculate rotation\n"));
+ if (calc_rotation(s)) {
+ if (s->flags & SI_SHOW_LINES) {
+ s->flags &= ~SI_SHOW_ROT; /* Calc rotation failed! */
+ show_lines(s);
+ }
+ goto sierr; /* Error */
+ }
+
+ if (s->flags & SI_BUILD_REF) { /* If generating a chart reference file */
+ /* Calculate the edge lists and write it to the file */
+ if (s->verb >= 2)
+ DBG((dbgo,"About to build feature information\n"));
+ if (calc_elists(s, 1)) /* reference */
+ goto sierr; /* Error */
+
+ if (s->verb >= 2)
+ DBG((dbgo,"About to write feature reference information\n"));
+ if (write_elists(s))
+ goto sierr; /* Error */
+ } else {
+ /* If we are matching to the reference and generating an output data file */
+ int rv;
+
+ /* Calculate the edge lists read for a match */
+ if (s->verb >= 2)
+ DBG((dbgo,"About to calculate feature information\n"));
+ if (calc_elists(s, 0)) /* match */
+ goto sierr; /* Error */
+
+ if (s->verb >= 2)
+ DBG((dbgo,"About to read reference feature information\n"));
+ if (read_relists(s))
+ goto sierr; /* Error */
+ if (s->verb >= 2)
+ DBG((dbgo,"Read of chart reference file succeeded\n"));
+
+ if (sfid != NULL) { /* Manual matching */
+ if (s->verb >= 2)
+ DBG((dbgo,"Using manual matching\n"));
+
+ if (s->havefids == 0) {
+ s->errv = SI_NO_FIDUCIALS_ERR;
+ sprintf(s->errm,"Chart recognition definition file doesn't contain fiducials");
+ goto sierr; /* Error */
+ }
+ if (compute_man_ptrans(s, sfid))
+ goto sierr;
+
+ /* Do the actual scan given out manual transformation matrix */
+ if (s->verb >= 2)
+ DBG((dbgo,"About to setup value scanrdg boxes\n"));
+ if (setup_sboxes(s))
+ goto sierr;
+ if (s->verb >= 2)
+ DBG((dbgo,"About to read raster values\n"));
+ if (do_value_scan(s))
+ goto sierr;
+
+ } else { /* Automatic matching */
+
+ /* Attempt to match input file with reference */
+ if (s->verb >= 2)
+ DBG((dbgo,"About to match features\n"));
+ if ((rv = do_match(s)) != 0) {
+ if (rv == 1) { /* No reasonable rotation found */
+ s->errv = SI_POOR_MATCH;
+ sprintf(s->errm,"Pattern match wasn't good enough");
+ }
+ goto sierr;
+ }
+
+ /* If there is patch matching data and more than one */
+ /* feasible matching rotation, try and discriminate between them. */
+ if (s->xpt && s->norots > 1) {
+ int i, j;
+ int flags = s->flags;
+
+ s->flags &= ~SI_SHOW_SAMPLED_AREA; /* Don't show areas for trials */
+
+ /* For each candidate rotation, scan in the pixel values */
+ for (s->crot = 0; s->crot < s->norots; s->crot++) {
+
+ /* Compute transformation from reference to input file */
+ if (s->verb >= 2)
+ DBG((dbgo,"About to compute match transform for rotation %f deg.\n",
+ DEG(s->rots[s->crot].irot)));
+ if (compute_ptrans(s))
+ goto sierr;
+
+ /* Setup the input boxes ready for scanning in the input values */
+ if (s->verb >= 2)
+ DBG((dbgo,"About to setup value scanrdg boxes\n"));
+ if (setup_sboxes(s))
+ goto sierr;
+
+ /* Scan in the pixel values */
+ if (s->verb >= 2)
+ DBG((dbgo,"About to read raster values\n"));
+ if (do_value_scan(s))
+ goto sierr;
+
+ /* Copy to this rotation values so that the best can be restored */
+ if (s->xpt != 0) { /* Got expected patch values to compare with */
+ if (s->verb >= 2)
+ DBG((dbgo,"About to compute expected value correlation\n"));
+ if (compute_xcc(s))
+ goto sierr;
+ }
+ }
+
+ /* Pick the best from the candidate rotation */
+ if (s->verb >= 2) {
+ DBG((dbgo,"Expected value distance values are:\n"));
+ for (i = 0; i < s->norots; i++) {
+ DBG((dbgo,"%d, rot %f: %f\n", i, DEG(s->rots[i].irot), s->rots[i].xcc));
+ }
+ }
+
+ for (j = 0, i = 1; i < s->norots; i++) {
+ if (s->rots[i].xcc < s->rots[j].xcc)
+ j = i;
+ }
+
+ if (s->verb >= 2)
+ DBG((dbgo,"Chosen rotation %f deg. as best\n",DEG(s->rots[j].irot)));
+
+ s->crot = j;
+ s->flags = flags; /* Restore flags */
+ }
+
+ /* Setup transformation to be that for chosen rotation for diagnostics */
+ if (s->verb >= 2)
+ DBG((dbgo,"About to compute final match transform\n"));
+ if (compute_ptrans(s))
+ goto sierr;
+
+ if (s->verb >= 2)
+ DBG((dbgo,"Improve match\n"));
+ if (improve_match(s))
+ goto sierr;
+
+ /* After choosing rotation of improving the fit, rescan the values */
+ if (s->verb >= 2)
+ DBG((dbgo,"About to setup value scanrdg boxes\n"));
+ if (setup_sboxes(s))
+ goto sierr;
+ if (s->verb >= 2)
+ DBG((dbgo,"About to read raster values\n"));
+ if (do_value_scan(s))
+ goto sierr;
+ }
+
+ if (s->flags & SI_SHOW_SBOX) {
+ show_sbox(s); /* Draw sample box outlines on diagnostic raster */
+ }
+ }
+ if (s->flags & SI_SHOW_LINES)
+ if(show_lines(s))
+ goto sierr; /* Error */
+
+sierr:;
+ if (s->verb >= 2)
+ DBG((dbgo,"About to write diag file\n"));
+ if (scanrd_write_diag(s))
+ return (scanrd *)s; /* Error */
+
+ return (scanrd *)s;
+}
+
+
+/********************************************************************************/
+
+/* Allocate the basic scanrd object */
+/* Return NULL on failure to allocate */
+/* Need to check errv for other problems */
+static scanrd_
+*new_scanrd(
+ int flags, /* option flags */
+ int verb, /* verbosity level */
+ double gammav, /* Approximate gamma encoding of image (0.0 = default 2.2) */
+ int (*write_line)(void *ddata, int y, char *src), /* Write RGB line of diag file */
+ void *ddata, /* Opaque data for write_line() */
+ int w, int h, /* Width and Height of input raster in pixels */
+ int d, int td, int p, /* Useful plane Depth, Total depth, Bit presision of input pixels */
+ int (*read_line)(void *fdata, int y, char *dst), /* Read RGB line of source file */
+ void *fdata, /* Opaque data for read_line() */
+
+ char *refname /* reference file name */
+) {
+ scanrd_ *s;
+
+ if ((s = (scanrd_ *)calloc(1, sizeof(scanrd_))) == NULL)
+ return NULL;
+
+ /* Public functions */
+ s->public.reset = scanrd_reset;
+ s->public.read = scanrd_read;
+ s->public.error = scanrd_error;
+ s->public.free = free_scanrd;
+
+ if (flags & (SI_SHOW_ROT | SI_SHOW_PERS | SI_SHOW_IMPL | SI_SHOW_ALL_LINES))
+ flags |= SI_SHOW_LINES; /* Key all line stuff off SI_SHOW_LINES */
+
+ if (flags & (SI_SHOW_SBOX_OUTLINES | SI_SHOW_SBOX_NAMES | SI_SHOW_SBOX_AREAS))
+ flags |= SI_SHOW_SBOX;; /* Key all sample box stuff off SI_SHOW_SBOX */
+
+ if (write_line == NULL)
+ flags &= ~SI_SHOW_FLAGS; /* If no diag file, turn off show flags */
+
+ s->flags = flags;
+ s->verb = verb;
+
+ s->errv = 0;
+ s->errm[0] = '\0';
+
+ if (gammav <= 0.0)
+ gammav = 2.2; /* default */
+ s->gammav = gammav;
+ s->width = w;
+ s->height = h;
+ s->depth = d;
+ s->tdepth = td;
+ s->bpp = p;
+
+ if (d > MXDE) {
+ s->errv = SI_PIX_DEPTH_ERR;
+ sprintf(s->errm,"scanrd: Pixel depth is too large");
+ return s;
+ }
+
+ if (p != 8 && p != 16) {
+ s->errv = SI_BIT_DEPTH_ERR;
+ sprintf(s->errm,"scanrd: Pixel bits/pixel is not 8 or 16");
+ return s;
+ }
+ if (p == 8)
+ s->bypp = 1;
+ else
+ s->bypp = 2;
+
+ if (verb >= 2)
+ DBG((dbgo,"Verbosity = %d, flags = 0x%x\n",verb, flags));
+
+ /* RGB Diagnostic output raster array requested */
+ if ((flags & SI_SHOW_FLAGS) && write_line != NULL) {
+ if ((s->out = malloc(3 * w * h)) == NULL) {
+ s->errv = SI_MALLOC_DIAG_RAST;
+ sprintf(s->errm,"scanrd: Diagnostic output raster array malloc failed");
+ return s;
+ }
+ }
+
+ s->noslines = 0;
+ s->novlines = 0;
+ s->gdone = NULL;
+ s->irot = 0.0;
+ s->norots = 0;
+
+ s->ppc[0] = 0.0;
+ s->ppc[1] = 0.0;
+ s->ppc[2] = 0.0;
+ s->ppc[3] = 0.0;
+
+ /* Set overall perspective transform to null */
+ s->ptrans[0] = 1.0;
+ s->ptrans[1] = 0.0;
+ s->ptrans[2] = 0.0;
+ s->ptrans[3] = 0.0;
+ s->ptrans[4] = 1.0;
+ s->ptrans[5] = 0.0;
+ s->ptrans[6] = 0.0;
+ s->ptrans[7] = 0.0;
+
+ INIT_ELIST(s->xelist);
+ INIT_ELIST(s->yelist);
+ INIT_ELIST(s->ixelist);
+ INIT_ELIST(s->iyelist);
+ INIT_ELIST(s->rxelist);
+ INIT_ELIST(s->ryelist);
+ s->rbox_shrink = 0.9;
+ s->xpt = 0;
+
+ s->nsbox = 0;
+ s->sboxes = NULL;
+ s->sbstart = NULL;
+ s->sbend = NULL;
+ s->csi = 0;
+ s->cei = 0;
+ s->alist = NULL;
+
+ s->next_read = 0;
+
+ s->refname = refname;
+
+ s->inited = 0;
+ s->vrego = s->vregn = NULL;
+ s->no_vo = s->no_vn = 0;
+ s->hrego = s->hregn = NULL;
+ s->no_ho = s->no_hn = 0;
+ s->th = TH;
+ s->divval = 0.25;
+ s->adivval = 0.0;
+ s->divc = 0;
+
+ /* aa line init */
+ s->aa_inited = 0; /* Let line init do the rest */
+ s->coverage = NULL;
+
+ /* Callbacks */
+ s->read_line = read_line;
+ s->fdata = fdata;
+
+ s->write_line = write_line;
+ s->ddata = ddata;
+
+ return s;
+}
+
+static void free_elist_array(elist *el);
+
+/* Free the object up */
+static void
+free_scanrd(
+scanrd *ps
+) {
+ scanrd_ *s = (scanrd_ *)ps; /* Cast public to private */
+ points *tp;
+
+ free_elist_array(&s->xelist);
+ free_elist_array(&s->yelist);
+ free_elist_array(&s->ixelist);
+ free_elist_array(&s->iyelist);
+ free_elist_array(&s->rxelist);
+ free_elist_array(&s->ryelist);
+
+ if (s->sboxes != NULL)
+ free(s->sboxes);
+ if (s->sbstart != NULL)
+ free(s->sbstart);
+ if (s->sbend != NULL)
+ free(s->sbend);
+ s->alist = NULL;
+
+ /* Free up done line list */
+ tp = s->gdone;
+ FOR_ALL_ITEMS(points, tp)
+ if (tp->r != NULL)
+ free(tp->r);
+ free(tp);
+ END_FOR_ALL_ITEMS(tp);
+ s->gdone = NULL;
+
+ /* Points were deleted with gdone ??? */
+ if(s->vrego != NULL)
+ free(s->vrego);
+ if(s->vregn)
+ free(s->vregn);
+ if(s->hrego != NULL)
+ free(s->hrego);
+ if(s->hregn != NULL)
+ free(s->hregn);
+ s->inited = 1;
+
+ /* Free up output diag array */
+ if (s->out != NULL)
+ free(s->out);
+
+ /* Free up aa line array */
+ if (s->coverage != NULL)
+ free(s->coverage);
+ free(s);
+}
+
+
+/* Return the error flag, and set the message pointer */
+static unsigned int
+scanrd_error(scanrd *ps, char **errm) {
+ scanrd_ *s = (scanrd_ *)ps; /* Cast public to private */
+ *errm = s->errm;
+ return s->errv;
+}
+
+/********************************************************************************/
+static int analize(scanrd_ *s, unsigned char *inp[6], int y);
+
+/* Read in and process the input file */
+/* Return non-zero on error */
+static int
+read_input(scanrd_ *s) {
+ unsigned char *in[6]; /* Pointer to six input buffers */
+ int w = s->width; /* Raster width */
+ int h = s->height; /* Raster height */
+ int i, y;
+
+ /* Allocate input line buffers */
+ for (i = 0; i < 6; i++) {
+ if ((in[i] = malloc(s->tdepth * w * s->bypp)) == NULL) {
+ s->errv = SI_MALLOC_INPUT_BUF;
+ sprintf(s->errm,"scanrd: Failed to malloc input line buffers");
+ return 1;
+ }
+ }
+
+ /* Prime the input buffers with 5 lines */
+ for (y = 0; y < 5; y++) {
+ if (s->read_line(s->fdata, y, (char *)in[y])) {
+ s->errv = SI_RAST_READ_ERR;
+ sprintf(s->errm,"scanrd: read_line() returned error");
+ return 1;
+ }
+ }
+ /* Process the tiff file line by line (Assume at least 6 lines in total raster) */
+ for (; y < h; ++y) {
+ unsigned char *tt;
+ if (s->read_line(s->fdata, y, (char *)in[5])) {
+ s->errv = SI_RAST_READ_ERR;
+ sprintf(s->errm,"scanrd: read_line() returned error");
+ return 1;
+ }
+
+ if (analize(s, in, y)) {
+ return 1;
+ }
+
+ tt = in[0]; /* Shuffle buffers about */
+ in[0] = in[1];
+ in[1] = in[2];
+ in[2] = in[3];
+ in[3] = in[4];
+ in[4] = in[5];
+ in[5] = tt;
+
+ }
+ s->adivval /= (double)s->divc; /* Average divider value, 1.0 = 0 degrees, 0.0 = 45 degrees */
+ if (s->adivval < 0.0)
+ s->adivval = 0.0;
+ else if (s->adivval > 1.0)
+ s->adivval = 1.0;
+
+ if (s->verb >= 2)
+ DBG((dbgo,"adivval = %f\n",s->adivval));
+
+ /* Free the input line buffers */
+ for (i = 0; i < 6; i++)
+ free(in[i]);
+
+ return 0;
+}
+
+/********************************************************************************/
+
+#ifdef NEVER /* Before 22/5/2004 */
+#define THRN 1.0 /* Threshold above average ratio - numerator */
+#define THRD 2.0 /* Threshold above average ratio - denominator */
+
+#define THAWF 1.0 /* Threshold average adaptation filter weight, fixed value (TH) */
+#define THAWP 4.0 /* Threshold average adaptation filter weight, previous value */
+#define THAWN 1.0 /* Threshold average adaptation filter weight, new value */
+
+#else /* Current values */
+
+#define THRN 1.0 /* Threshold above average ratio - numerator */
+#define THRD 1.5 /* Threshold above average ratio - denominator */
+
+#define THAWF 1.0 /* Threshold average adaptation filter weight, fixed value (TH) */
+#define THAWP 5.0 /* Threshold average adaptation filter weight, previous value */
+#define THAWN 1.0 /* Threshold average adaptation filter weight, new value */
+
+#endif
+
+/* ~~~ minimum raster size needs to be specified/checked ~~~~ */
+#define MIN_NO_LINES 16 /* Minimum number of valid fitted lines to estimate rotation */
+
+/* Criteria for accepting lines for angle calculation (valid lines) */
+#define MAX_MWID_TO_LEN 0.1
+#define MIN_POINT_TO_AREA 0.9 /* Minimum point desity over the lines area */
+#define SD_WINDOW 1.5 /* Allow += 1.5 of a standard deviation for robust angle calc. */
+#define ELISTCDIST 800 /* 1/ELISTCDIST = portion of refence edge list legth to coalesce over */
+
+/* Criteria for accepting lines for improring final fit */
+#define IMP_MATCH 0.10 /* Proportion of average tick spacing */
+
+/* The following should be scaled to the resolution of the image ? */
+#define MIN_POINTS 10 /* Minimum points to calculate line */
+#define MIN_LINE_LENGTH 10.0
+#define CUT_CHUNKS 128 /* cut groups along diagonals - must be power of 2 */
+
+static int add_region(scanrd_ *s, region *rego, int no_o, region *regn, int no_n, int y);
+
+/* Process a line of the TIFF file */
+/* return non-zero on error */
+static int
+analize(
+scanrd_ *s,
+unsigned char *inp[6], /* current and previous 5 lines */
+int y /* Current line y */
+) {
+ int w = s->width;
+ int stride = s->tdepth * s->width; /* In pixels */
+ unsigned short *gamma = s->gamma;
+ int x,i;
+ unsigned short *inp2[6]; /* current and previous 5 lines (16bpp) equivalent of inp[] */
+ unsigned char *in[6]; /* six input lines (8bpp) */
+ unsigned short *in2[6]; /* six input lines (16bpp) */
+ region *tr;
+ double tdh,tdv; /* Horizontal/virtical detect levels */
+ double tdmag;
+ double atdmag = 0.0; /* Average magnitude over a line */
+ int atdmagc = 0; /* Average magnitude over a line count */
+ double linedv = 0.0; /* Lines average divider value */
+ int linedc = 0; /* Lines average count */
+ int xo3 = s->tdepth * 3; /* Xoffset by 3 pixels */
+ int xo2 = s->tdepth * 2; /* Xoffset by 2 pixels */
+ int xo1 = s->tdepth * 1; /* Xoffset by 1 pixels */
+
+ for (x = 0; x < 6; x++) /* Create 16 bpp version of line pointers */
+ inp2[x] = (unsigned short *)inp[x];
+
+ if (s->inited == 0) {
+ /* Init gamma conversion lookup and region tracking. */
+ /* The assumption is that a typical chart has an approx. visually */
+ /* uniform distribution of samples, so that a typically gamma */
+ /* encoded scan image will have an average pixel value of 50%. */
+ /* If a the chart has a different gamma encoding (ie. linear), */
+ /* then we convert it to gamma 2.2 encoded to (hopefuly) enhance */
+ /* the patch contrast. */
+ if (s->bpp == 8)
+ for (i = 0; i < 256; i++) {
+ int byteb1;
+
+ byteb1 = (int)(0.5 + 255 * pow( i / 255.0, s->gammav/2.2 ));
+ gamma[i] = byteb1;
+ }
+ else
+ for (i = 0; i < 65536; i++) {
+ int byteb1;
+
+ byteb1 = (int)(0.5 + 65535 * pow( i / 65535.0, s->gammav/2.2 ));
+ gamma[i] = byteb1;
+ }
+
+ if ((s->vrego = (region *) malloc(sizeof(region) * (w+1)/2)) == NULL) {
+ s->errv = SI_MALLOC_VREGION;
+ sprintf(s->errm,"vreg malloc failed");
+ return 1;
+ }
+ s->no_vo = 0;
+ if ((s->vregn = (region *) malloc(sizeof(region) * (w+1)/2)) == NULL) {
+ s->errv = SI_MALLOC_VREGION;
+ sprintf(s->errm,"vreg malloc failed");
+ return 1;
+ }
+ s->no_vn = 0;
+ if ((s->hrego = (region *) malloc(sizeof(region) * (w+1)/2)) == NULL) {
+ s->errv = SI_MALLOC_VREGION;
+ sprintf(s->errm,"vreg malloc failed");
+ return 1;
+ }
+ s->no_ho = 0;
+ if ((s->hregn = (region *) malloc(sizeof(region) * (w+1)/2)) == NULL) {
+ s->errv = SI_MALLOC_VREGION;
+ sprintf(s->errm,"vreg malloc failed");
+ return 1;
+ }
+ s->no_hn = 0;
+ INIT_LIST(s->gdone);
+ s->inited = 1;
+ }
+
+ /* Un-gamma correct the latest input line */
+ if (s->bpp == 8)
+ for (x = 0; x < stride; x++)
+ inp[5][x] = (unsigned char)gamma[inp[5][x]];
+ else
+ for (x = 0; x < stride; x++)
+ inp2[5][x] = gamma[inp2[5][x]];
+
+ /* Compute difference output for line y-3 */
+ atdmagc = w - 5; /* Magnitude count (to compute average) */
+ for (x = 3; x < (w-2); x++) { /* Allow for -3 to +2 from x */
+ unsigned char *out = s->out;
+ int e;
+ int ss;
+ int idx = ((y-2) * w + x) * 3; /* Output raster index in bytes */
+
+ if (s->bpp == 8)
+ for (i = 0; i < 6; i++)
+ in[i] = inp[i] + x * s->tdepth; /* Strength reduce */
+ else
+ for (i = 0; i < 6; i++) {
+ in2[i] = inp2[i] + x * s->tdepth; /* Strength reduce */
+ in[i] = (unsigned char *)in2[i]; /* track 8bpp pointers */
+ }
+
+ if (s->flags & SI_SHOW_IMAGE) { /* Create B&W image */
+ toRGB(out + idx, in[2], s->depth, s->bpp); /* Convert to RGB */
+ out[idx] = out[idx+1] = out[idx+2] = (2 * out[idx] + 7 * out[idx+1] + out[idx+2])/10;
+ }
+
+ ss = 0; /* Sign of cross components the same vote */
+ tdh = tdv = 0.0;
+
+ if (s->bpp == 8)
+ for (e = 0; e < s->depth; e++) {
+ int d1,d2;
+ /* Compute Gxp */
+ d1 = -in[0][-xo3+e] + -in[0][-xo2+e] + -in[0][-xo1+e]
+ + -in[0][ 0+e] + -in[0][ xo1+e] + -in[0][ xo2+e]
+ + -in[1][-xo3+e] + -in[1][-xo2+e] + -in[1][-xo1+e]
+ + -in[1][ 0+e] + -in[1][ xo1+e] + -in[1][ xo2+e]
+ + -in[2][-xo3+e] + -in[2][-xo2+e] + -in[2][-xo1+e]
+ + -in[2][ 0+e] + -in[2][ xo1+e] + -in[2][ xo2+e]
+ + in[3][-xo3+e] + in[3][-xo2+e] + in[3][-xo1+e]
+ + in[3][ 0+e] + in[3][ xo1+e] + in[3][ xo2+e]
+ + in[4][-xo3+e] + in[4][-xo2+e] + in[4][-xo1+e]
+ + in[4][ 0+e] + in[4][ xo1+e] + in[4][ xo2+e]
+ + in[5][-xo3+e] + in[5][-xo2+e] + in[5][-xo1+e]
+ + in[5][ 0+e] + in[5][ xo1+e] + in[5][ xo2+e];
+ /* Compute Gyp */
+ d2 = -in[0][-xo3+e] + -in[1][-xo3+e] + -in[2][-xo3+e]
+ + -in[3][-xo3+e] + -in[4][-xo3+e] + -in[5][-xo3+e]
+ + -in[0][-xo2+e] + -in[1][-xo2+e] + -in[2][-xo2+e]
+ + -in[3][-xo2+e] + -in[4][-xo2+e] + -in[5][-xo2+e]
+ + -in[0][-xo1+e] + -in[1][-xo1+e] + -in[2][-xo1+e]
+ + -in[3][-xo1+e] + -in[4][-xo1+e] + -in[5][-xo1+e]
+ + in[0][ 0+e] + in[1][ 0+e] + in[2][ 0+e]
+ + in[3][ 0+e] + in[4][ 0+e] + in[5][ 0+e]
+ + in[0][+xo1+e] + in[1][+xo1+e] + in[2][+xo1+e]
+ + in[3][+xo1+e] + in[4][+xo1+e] + in[5][+xo1+e]
+ + in[0][+xo2+e] + in[1][+xo2+e] + in[2][+xo2+e]
+ + in[3][+xo2+e] + in[4][+xo2+e] + in[5][+xo2+e];
+
+ if ((d1 >= 0 && d2 >=0)
+ || (d1 < 0 && d2 < 0))
+ ss++; /* Sign was the same */
+ tdh += d1/4.5 * d1/4.5; /* (4.5 = 6x6/4x2, to scale original tuned values) */
+ tdv += d2/4.5 * d2/4.5;
+ }
+ else
+ for (e = 0; e < s->depth; e++) {
+ int d1,d2;
+ /* Compute Gxp */
+ d1 = -in2[0][-xo3+e] + -in2[0][-xo2+e] + -in2[0][-xo1+e]
+ + -in2[0][ 0+e] + -in2[0][ xo1+e] + -in2[0][ xo2+e]
+ + -in2[1][-xo3+e] + -in2[1][-xo2+e] + -in2[1][-xo1+e]
+ + -in2[1][ 0+e] + -in2[1][ xo1+e] + -in2[1][ xo2+e]
+ + -in2[2][-xo3+e] + -in2[2][-xo2+e] + -in2[2][-xo1+e]
+ + -in2[2][ 0+e] + -in2[2][ xo1+e] + -in2[2][ xo2+e]
+ + in2[3][-xo3+e] + in2[3][-xo2+e] + in2[3][-xo1+e]
+ + in2[3][ 0+e] + in2[3][ xo1+e] + in2[3][ xo2+e]
+ + in2[4][-xo3+e] + in2[4][-xo2+e] + in2[4][-xo1+e]
+ + in2[4][ 0+e] + in2[4][ xo1+e] + in2[4][ xo2+e]
+ + in2[5][-xo3+e] + in2[5][-xo2+e] + in2[5][-xo1+e]
+ + in2[5][ 0+e] + in2[5][ xo1+e] + in2[5][ xo2+e];
+ /* Compute Gyp */
+ d2 = -in2[0][-xo3+e] + -in2[1][-xo3+e] + -in2[2][-xo3+e]
+ + -in2[3][-xo3+e] + -in2[4][-xo3+e] + -in2[5][-xo3+e]
+ + -in2[0][-xo2+e] + -in2[1][-xo2+e] + -in2[2][-xo2+e]
+ + -in2[3][-xo2+e] + -in2[4][-xo2+e] + -in2[5][-xo2+e]
+ + -in2[0][-xo1+e] + -in2[1][-xo1+e] + -in2[2][-xo1+e]
+ + -in2[3][-xo1+e] + -in2[4][-xo1+e] + -in2[5][-xo1+e]
+ + in2[0][ 0+e] + in2[1][ 0+e] + in2[2][ 0+e]
+ + in2[3][ 0+e] + in2[4][ 0+e] + in2[5][ 0+e]
+ + in2[0][+xo1+e] + in2[1][+xo1+e] + in2[2][+xo1+e]
+ + in2[3][+xo1+e] + in2[4][+xo1+e] + in2[5][+xo1+e]
+ + in2[0][+xo2+e] + in2[1][+xo2+e] + in2[2][+xo2+e]
+ + in2[3][+xo2+e] + in2[4][+xo2+e] + in2[5][+xo2+e];
+
+ if ((d1 >= 0 && d2 >=0)
+ || (d1 < 0 && d2 < 0))
+ ss++; /* Sign was the same */
+
+ tdh += d1/(4.5 * 257) * d1/(4.5 * 257); /* Scale to 0..255 range */
+ tdv += d2/(4.5 * 257) * d2/(4.5 * 257);
+ }
+
+ tdmag = tdh + tdv;
+
+ if (tdmag < (32.0 * s->th))
+ atdmag += tdmag; /* Average magnitude over a line */
+ else
+ atdmag += 32.0 * s->th;
+
+ /* if over threshold */
+ /* (Cut long lines up to prevent long lines being */
+ /* (thrown away due to attached blobs) */
+ if (tdmag >= s->th
+ && (x & (CUT_CHUNKS-1)) != (y & (CUT_CHUNKS-1))) {
+ double tt;
+ double av; /* Angle value of current pixel */
+ tt = (tdv - tdh)/(tdh + tdv); /* Partial angle */
+ linedv += fabs(tt);
+ linedc++;
+
+ if (ss >= (s->depth/2+1)) /* Assume signs are the same if clear majority */
+ av = 3.0 + tt;
+ else
+ av = 1.0 - tt;
+
+ /* Separate the orthogonal elements */
+ if (av >= s->divval && av < (s->divval + 2.0)) {
+ if (s->flags & SI_SHOW_DIFFSH)
+ out[idx] = (char)255; /* Red */
+ /* Add point to new region */
+ /* See if we can add to last region */
+ if (s->no_hn > 0 && x == s->hregn[s->no_hn-1].hx)
+ s->hregn[s->no_hn-1].hx++;
+ else { /* Add another */
+ if (s->no_hn >= (w+1)/2) {
+ s->errv = SI_INTERNAL;
+ sprintf(s->errm,"Internal, no_hn is too large");
+ return 1;
+ }
+ s->hregn[s->no_hn].lx = x;
+ s->hregn[s->no_hn].hx = x+1;
+ s->hregn[s->no_hn].p = NULL;
+ s->no_hn++;
+ }
+ } else {
+ if (s->flags & SI_SHOW_DIFFSV)
+ out[idx+1] = (char)255; /* Green */
+ /* Add point to new region */
+ /* See if we can add to last region */
+ if (s->no_vn > 0 && x == s->vregn[s->no_vn-1].hx)
+ s->vregn[s->no_vn-1].hx++;
+ else { /* Add another */
+ if (s->no_vn >= (w+1)/2) {
+ s->errv = SI_INTERNAL;
+ sprintf(s->errm,"Internal, no_vn is too large");
+ return 1;
+ }
+ s->vregn[s->no_vn].lx = x;
+ s->vregn[s->no_vn].hx = x+1;
+ s->vregn[s->no_vn].p = NULL;
+ s->no_vn++;
+ }
+ }
+ }
+ }
+
+ if (linedc != 0) { /* Adapt divider value to line */
+ linedv /= (double)linedc; /* Compute average over the line */
+ linedv = (linedv * linedv); /* Square to even out linedv vs angle */
+ linedv = (1.65 * (linedv - 0.12)); /* Compensate for random offsets */
+ s->adivval += linedv;
+ s->divc++;
+ s->divval = (7.0 * s->divval + linedv)/8.0; /* Average over 8 lines */
+ if (s->divval < 0.0)
+ s->divval = 0.0;
+ else if (s->divval > 1.0)
+ s->divval = 1.0;
+ if (s->verb >= 5)
+ DBG((dbgo,"linedv = %f, divval = %f\n",linedv,s->divval));
+ }
+
+ /* Adjust the threshold */
+ atdmag /= (double)atdmagc; /* compute average magnitude over the line */
+ s->th = (s->th * THRD)/(THRN + s->divval);/* Convert threshold to average */
+ s->th = ((THAWF * TH) + (THAWP * s->th) + (THAWN * atdmag))/(THAWF + THAWP + THAWN);
+ s->th = (s->th * (THRN + s->divval))/THRD; /* Convert average back to threshold */
+
+ /* Add vertical regions */
+ if (add_region(s,s->vrego,s->no_vo,s->vregn,s->no_vn,y-2))
+ return 1;
+
+ /* Add horizontal regions */
+ if (add_region(s,s->hrego,s->no_ho,s->hregn,s->no_hn,y-2))
+ return 1;
+
+ /* shuffle them along */
+ tr = s->vrego;
+ s->vrego = s->vregn; /* move new to old */
+ s->vregn = tr; /* old to new */
+ s->no_vo = s->no_vn;
+ s->no_vn = 0;
+
+ tr = s->hrego;
+ s->hrego = s->hregn; /* move new to old */
+ s->hregn = tr; /* old to new */
+ s->no_ho = s->no_hn;
+ s->no_hn = 0;
+
+ return 0;
+}
+
+/********************************************************************************/
+/* Point list code */
+
+/* allocate a new (empty) points structure */
+/* return NULL on error */
+static points *
+new_points(
+scanrd_ *s
+) {
+ points *ps;
+ static int pn = 0;
+ if ((ps = (points *) malloc(sizeof(points))) == NULL) {
+ s->errv = SI_MALLOC_POINTS;
+ sprintf(s->errm,"new_points: malloc failed");
+ return NULL;
+ }
+ ps->mxno = 0;
+ ps->no = 0;
+ ps->nop = 0;
+ ps->r = NULL;
+ ps->pn = pn;
+ pn++;
+ return ps;
+}
+
+/* destroy a points structure */
+static void
+destroy_points(
+scanrd_ *s,
+points *ps) {
+ if (ps->r != NULL) /* Free any array pointed to */
+ free(ps->r);
+ free (ps);
+}
+
+/* Add another run to a points object */
+/* return non-zero on error */
+static int
+add_run(
+scanrd_ *s,
+points *ps,
+int lx,
+int hx,
+int y)
+ {
+ if (ps->no == ps->mxno) { /* Need some more space */
+ ps->mxno = (2 * ps->mxno) + 5; /* New size */
+ if ((ps->r = (run *) realloc(ps->r, sizeof(run) * ps->mxno)) == NULL) {
+ s->errv = SI_REALLOC_POINTS;
+ sprintf(s->errm,"add_run: realloc failed");
+ return 1;
+ }
+ }
+ ps->r[ps->no].lx = lx;
+ ps->r[ps->no].hx = hx;
+ ps->r[ps->no].y = y;
+ ps->no++; /* One more run */
+ ps->nop += hx - lx; /* Total of pixels */
+ return 0;
+}
+
+/* copy src points to dest */
+/* Return non-zero on error */
+static int
+copy_points(
+scanrd_ *s,
+points *dst,
+points *src
+) {
+ int i;
+ for (i = 0; i < src->no; i++) {
+ if (add_run(s,dst,src->r[i].lx,src->r[i].hx,src->r[i].y))
+ return 1;
+ }
+ return 0;
+}
+
+/********************************************************************************/
+
+/* Add a new region of points to the line points lists */
+/* Note that regions are assumed to be non-overlapping x sorted */
+/* Return non-zero on error */
+static int
+add_region(
+scanrd_ *s,
+region *rego, /* Old regions */
+int no_o, /* No of old region */
+region *regn, /* New regions */
+int no_n, /* No of new region */
+int y /* Y value */
+) {
+ int osp,op,np; /* Old/new pointers */
+
+ osp = 0;
+ for (np = 0; np < no_n; np++) { /* Process all new runs */
+ /* Advance start pointer until we get to runs that may touch */
+#ifdef DIAGN
+ while (osp < no_o && rego[osp].hx < regn[np].lx)
+#else
+ while (osp < no_o && rego[osp].hx <= regn[np].lx)
+#endif
+ osp++;
+ /* For all old runs that may touch new */
+#ifdef DIAGN
+ for(op = osp; op < no_o && rego[op].lx <= regn[np].hx; op++) {
+#else
+ for(op = osp; op < no_o && rego[op].lx < regn[np].hx; op++) {
+#endif
+
+#ifdef DIAGN
+ if (rego[op].hx >= regn[np].lx && rego[op].lx <= regn[np].hx) {
+#else
+ if (rego[op].hx > regn[np].lx && rego[op].lx < regn[np].hx) {
+#endif
+ /* Old region touches new */
+ if (regn[np].p == NULL) { /* No group for new yet */
+ regn[np].p = rego[op].p; /* Make part of the same group */
+ if (add_run(s, regn[np].p,regn[np].lx,regn[np].hx,y)) /* add new run to group */
+ return 1;
+ } else if (regn[np].p != rego[op].p) { /* Touches different group */
+ int j;
+ points *tp = rego[op].p; /* Old region to be renamed/merged */
+ if (copy_points(s,regn[np].p,tp)) /* Merge old with current new */
+ return 1; /* Error */
+ DEL_LINK(s->gdone,tp); /* Don't need other any more */
+ for (j = 0; j < no_o; j++) /* Fix all references to this group */
+ if (rego[j].p == tp)
+ rego[j].p = regn[np].p;
+ for (j = 0; j < no_n; j++)
+ if (regn[j].p == tp)
+ regn[j].p = regn[np].p;
+ destroy_points(s,tp);
+ }
+ }
+ }
+ /* Finished all relevant old runs */
+ if (regn[np].p == NULL) { /* No old touched, so start new group */
+ if ((regn[np].p = new_points(s)) == NULL)
+ return 1; /* Error */
+ ADD_ITEM_TO_TOP(s->gdone,regn[np].p); /* Stash it in points list */
+ if (add_run(s, regn[np].p,regn[np].lx,regn[np].hx,y)) /* add new run to group */
+ return 1; /* Error */
+ }
+ }
+ return 0;
+}
+
+/********************************************************************************/
+
+/* Apply partial perspective to an xy point */
+/* (We omit the two offset parameters, since we don't need them) */
+void ppersp(scanrd_ *s, double *xx, double *yy, double x, double y, double *ppc) {
+ double den;
+
+ /* Offset the partial perspective transform */
+ x -= ppc[2];
+ y -= ppc[3];
+
+ den = ppc[0] * x + ppc[1] * y + 1.0;
+
+ if (fabs(den) < 1e-6) {
+ if (den < 0.0)
+ den = -1e-6;
+ else
+ den = 1e-6;
+ }
+ *xx = x/den + ppc[2];
+ *yy = y/den + ppc[3];
+}
+
+
+/* Apply inverse partial perspective to an xy point */
+void invppersp(scanrd_ *s, double *x, double *y, double xx, double yy, double *ppc) {
+ double den;
+
+ /* Offset the partial perspective transform */
+ xx -= ppc[2];
+ yy -= ppc[3];
+
+ den = - ppc[0] * xx - ppc[1] * yy + 1.0;
+
+ if (fabs(den) < 1e-6) {
+ if (den < 0.0)
+ den = -1e-6;
+ else
+ den = 1e-6;
+ }
+ *x = xx/den + ppc[2];
+ *y = yy/den + ppc[3];
+}
+
+/********************************************************************************/
+
+/* Compute the least squares best line fit for a group */
+/* Return non-zero if failed */
+static int
+points_to_line(
+scanrd_ *s,
+points *ps) {
+ int i,j;
+ point *vv; /* Point vectors */
+ int nop = ps->nop; /* Number of points */
+ double sx,sy; /* Sum */
+ double mx,my; /* Mean */
+ double a; /* Angle, Clockwise from 12o'clock */
+ double mw,len; /* mean width, length */
+ double x1,y1,x2,y2; /* Start/end point of fitted line */
+
+ ps->flag = 0;
+
+ if (nop < MIN_POINTS) /* Don't bother if too few pixels */
+ return 0;
+
+ /* Convert runs to individual points, and compute mean */
+ if ((vv = (point *) malloc(sizeof(point) * nop)) == NULL) {
+ s->errv = SI_MALLOC_POINT2LINE;
+ sprintf(s->errm,"scanrd: points_to_line: malloc failed");
+ return 1;
+ }
+
+ sx = sy = 0.0;
+ for (j = i = 0; i < ps->no; i++) { /* For all runs */
+ int x,y;
+ int hx = ps->r[i].hx, lx = ps->r[i].lx;
+
+ y = ps->r[i].y;
+ sy += (hx - lx) * y;
+ for (x = lx; x < hx; x++, j++) { /* Convert to points */
+ sx += x;
+ vv[j].x = x;
+ vv[j].y = y;
+ }
+ }
+ mx = sx/(double)nop; /* Centroid (mean) of points */
+ my = sy/(double)nop;
+
+ /* Offset points to centroid */
+ for (i=0; i < nop; i++) {
+ vv[i].x -= mx;
+ vv[i].y -= my;
+ }
+
+ /* Compute ad and bd, then A, B, C */
+ /* From Graphics Gems V, pp 91-97, */
+ /* "The Best Least-Squares Line Fit" */
+ /* by David Alciatore and Rick Miranda. */
+ {
+ double ad, bd; /* a' and b' values */
+ double xd, yd; /* temp x' and y' */
+ double A, B; /* line equation */
+ double abn; /* A & B normalizer */
+
+ xd = yd = bd = 0.0;
+ for (i = 0; i < nop; i++) {
+ double x, y;
+
+ x = vv[i].x;
+ y = vv[i].y;
+ xd += x * x;
+ yd += y * y;
+ bd += x * y;
+ }
+ ad = xd - yd;
+
+ /* Equation of best fit line is Ax + By = C */
+ A = 2 * bd;
+ B = -(ad + sqrt(ad * ad + 4.0 * bd * bd));
+ /* C = A * mx + B * my; */
+
+ /* Compute angle */
+ /* A = abn * cos(a), B = -abn * sin(a) */
+
+ abn = sqrt(A * A + B * B); /* Normalize A & B */
+ if (fabs(abn) < 1e-6) { /* No dominant direction */
+ a = 0.0;
+ } else {
+ a = acos(A/abn);
+ }
+ /* Make angle +ve */
+ while (a < 0.0) a += M_PI;
+ }
+
+ /* Now figure out the bounding box for the line + other stats */
+ {
+ double s,c;
+ double pl,nl; /* Positive length, negative length */
+ s = sin(a);
+ c = cos(a);
+ for (mw = 0.0, pl = 0.0, nl = 0.0, i = 0; i < nop; i++)
+ {
+ double npj; /* Projection onto normal */
+ double lpj; /* Projection onto line */
+ npj = -c * vv[i].x + s * vv[i].y;
+ if (npj < 0)
+ mw -= npj;
+ else
+ mw += npj;
+ lpj = s * vv[i].x + c * vv[i].y;
+ if (lpj > pl)
+ pl = lpj;
+ if (lpj < nl)
+ nl = lpj;
+ }
+ mw = 2.0 * mw/(double)nop; /* Mean width */
+
+ x1 = mx + s * nl;
+ y1 = my + c * nl;
+ x2 = mx + s * pl;
+ y2 = my + c * pl;
+ len = pl - nl;
+ }
+
+ ps->mx = mx; /* Mean point */
+ ps->my = my;
+ ps->a = a; /* Angle */
+ ps->mw = mw; /* Mean width */
+ ps->len = len; /* Mean length */
+ ps->x1 = x1; /* Start/end point of fitted line */
+ ps->y1 = y1;
+ ps->x2 = x2;
+ ps->y2 = y2;
+ ps->flag = F_LINESTATS; /* Line stats valid */
+
+ /* Compute the Constrained to 90 degrees angle */
+ /* We use the adivval to figure out where to split angles */
+ /* Split at 0 if adivval == 0.0, split at 45 if adivval == 1.0 */
+ if (a >= (M_PI * (1.0 - s->adivval/4.0)))
+ ps->ca = a - M_PI;
+ else if (a >= (M_PI * (0.5 - s->adivval/4.0)))
+ ps->ca = a - M_PI_2;
+ else
+ ps->ca = a;
+
+ if (s->verb >= 5)
+ DBG((dbgo,"Angle %f, CA = %f, length = %f, mean width = %f, Line %f,%f to %f,%f\n",
+ DEG(a),DEG(ps->ca),len,mw,x1,y1,x2,y2));
+ free(vv);
+
+/* printf("~~stats: mw = %f, len = %f, mw/len = %f, area = %f\n",
+ mw, len, mw/len, ((double)nop/(len * (mw + 0.01)))); */
+ /* Look at stats to see what lines are acceptable for further processing */
+ if ( len >= MIN_LINE_LENGTH
+ && mw/len <= MAX_MWID_TO_LEN
+ && ((double)nop/(len * (mw + 0.01))) >= MIN_POINT_TO_AREA) {
+ ps->flag |= F_VALID; /* Line stats valid to use */
+/* printf("~~set valid\n"); */
+ }
+ return 0;
+}
+
+static int
+calc_lines(
+scanrd_ *s
+) {
+ points *tp;
+ s->noslines = 0;
+ s->novlines = 0;
+ tp = s->gdone;
+ FOR_ALL_ITEMS(points, tp)
+ if (points_to_line(s,tp))
+ return 1; /* Error */
+ if (tp->flag & F_LINESTATS) /* Line stats valid */
+ s->noslines++;
+ if (tp->flag & F_VALID) /* Valid for angle calcs */
+ s->novlines++;
+
+ /* Save orininal raster (non partial perspective corrected) values */
+ if (tp->flag & F_VALID) {
+ tp->pmx = tp->mx;
+ tp->pmy = tp->my;
+ tp->px1 = tp->x1;
+ tp->py1 = tp->y1;
+ tp->px2 = tp->x2;
+ tp->py2 = tp->y2;
+ }
+ END_FOR_ALL_ITEMS(tp);
+ return 0;
+}
+
+static int show_line(scanrd_ *s, int x1, int y1, int x2, int y2, unsigned long c);
+
+/* Show the edge detected lines */
+static int
+show_lines(
+scanrd_ *s
+) {
+ points *tp;
+ int outw = s->width;
+ int outh = s->height;
+ /* For SI_SHOW_ROT */
+ double cirot,sirot; /* cos and sin of -irot */
+ cirot = cos(-s->irot);
+ sirot = sin(-s->irot);
+
+ tp = s->gdone;
+ FOR_ALL_ITEMS(points, tp)
+ if ((s->flags & SI_SHOW_ALL_LINES) || (tp->flag & F_VALID))
+ {
+ unsigned long col = 0xffffff; /* default color is white */
+ double x1 = tp->px1, y1 = tp->py1, x2 = tp->px2, y2 = tp->py2;
+ /* For SI_SHOW_ROT */
+
+ /* Show partial perspective corrected lines */
+ if (s->flags & (SI_SHOW_ROT | SI_SHOW_PERS)) {
+ invppersp(s, &x1, &y1, x1, y1, s->ppc);
+ invppersp(s, &x2, &y2, x2, y2, s->ppc);
+ col = 0xffff00; /* cyan */
+ }
+
+ /* Show rotation correction of lines + color coding yellow and red */
+ if (s->flags & SI_SHOW_ROT) {
+ double tx1, ty1, tx2, ty2;
+ double a = tp->a - s->irot;
+
+ tx1 = x1;
+ ty1 = y1;
+ tx2 = x2;
+ ty2 = y2;
+
+ /* Rotate about center of raster */
+ x1 = (tx1-outw/2.0) * cirot + (ty1-outh/2.0) * sirot;
+ y1 = -(tx1-outw/2.0) * sirot + (ty1-outh/2.0) * cirot;
+ x2 = (tx2-outw/2.0) * cirot + (ty2-outh/2.0) * sirot;
+ y2 = -(tx2-outw/2.0) * sirot + (ty2-outh/2.0) * cirot;
+
+ x1 += outw/2.0; /* Rotate about center of raster */
+ y1 += outh/2.0;
+ x2 += outw/2.0;
+ y2 += outh/2.0;
+ if ((a >= -0.08 && a <= 0.08) || (a >= (M_PI-0.08) && a <= (M_PI+0.08))
+ || (a >= (M_PI_2-0.08) && a <= (M_PI_2+0.08)))
+ col = 0x00ffff; /* yellow */
+ else
+ col = 0x0000ff; /* Red */
+ }
+ /* Show just lines used for fit improvement in blue */
+ if (s->flags & SI_SHOW_IMPL) {
+ if (tp->flag & F_IMPROVE)
+ col = 0xff4040; /* blue */
+ }
+ show_line(s,(int)(x1+0.5),(int)(y1+0.5),(int)(x2+0.5),(int)(y2+0.5),col);
+ }
+ END_FOR_ALL_ITEMS(tp);
+ return 0;
+}
+
+
+/********************************************************************************/
+
+/* Definition of the optimization function handed to powell() */
+static double
+pfunc(void *ss, double p[]) {
+ scanrd_ *s = (scanrd_ *)ss;
+ points *tp;
+ double aa; /* Average angle */
+ double va, rva; /* Variance */
+ double wt; /* Total weighting = sum of line lengths */
+ double pw;
+ double dw; /* Discrimination width */
+
+//printf("~1 %f %f %f %f %f %f\n", p[0],p[1],p[2],p[3],p[4],p[5]);
+
+ /* Correct the perspective of all the edge lines using the parameters */
+ /* and compute the mean angle */
+ aa = 0.0; /* Average constrained angle */
+ wt = 0.0; /* Total weighting = sum of line lengths */
+ tp = s->gdone;
+ FOR_ALL_ITEMS(points, tp)
+ if (tp->flag & F_LONGENOUGH) {
+ double a, ca;
+ invppersp(s, &tp->x1, &tp->y1, tp->px1, tp->py1, p);
+ invppersp(s, &tp->x2, &tp->y2, tp->px2, tp->py2, p);
+
+ /* Compute the angle */
+ a = atan2(tp->x2 - tp->x1,tp->y2 - tp->y1);
+
+ /* Make angle +ve */
+ while (a < 0.0)
+ a += M_PI;
+
+ /* Compute the Constrained to 90 degrees angle */
+ /* We use the adivval to figure out where to split angles */
+ /* Split at 0 if adivval == 0.0, split at 45 if adivval == 1.0 */
+ if (a >= (M_PI * (1.0 - s->adivval/4.0)))
+ ca = a - M_PI;
+ else if (a >= (M_PI * (0.5 - s->adivval/4.0)))
+ ca = a - M_PI_2;
+ else
+ ca = a;
+
+ tp->a = a;
+ tp->ca = ca;
+
+ aa += tp->len * ca;
+ wt += tp->len;
+ }
+ END_FOR_ALL_ITEMS(tp);
+ aa /= wt;
+
+ /* Calculate the angle variance */
+ va = 0.0;
+ tp = s->gdone;
+ wt = 0.0;
+ FOR_ALL_ITEMS(points, tp)
+ if (tp->flag & F_LONGENOUGH) {
+ double tt;
+ tt = tp->ca - aa;
+ va += tp->len * tt * tt;
+ wt += tp->len;
+ }
+ END_FOR_ALL_ITEMS(tp);
+ va = va/wt;
+
+ /* Calculate the a robust angle variance */
+ rva = 0.0;
+ wt = 0.0;
+ dw = sqrt(va) * 3.1; /* Allow += 0.5 of a standard deviation */
+ if (dw < 0.0001) /* A perfect chart may have dw of zero */
+ dw = 0.0001;
+ tp = s->gdone;
+ FOR_ALL_ITEMS(points, tp)
+ if (tp->flag & F_LONGENOUGH && fabs(tp->ca - aa) <= dw) {
+ double tt;
+ tt = tp->ca - aa;
+ rva += tp->len * tt * tt;
+ wt += tp->len;
+ }
+ END_FOR_ALL_ITEMS(tp);
+ if (wt > 0.0) {
+ rva = rva/wt;
+ va = rva;
+ }
+
+ /* Add some regularization to stop it going crazy */
+ pw = 0.0;
+ pw += 0.01 * (fabs(p[0]) + fabs(p[1]));
+ pw += 0.0001 * (fabs(p[2]/s->width - 0.5) + fabs(p[3]/s->height - 0.5));
+ va += pw;
+
+ return va;
+}
+
+/* Calculate the partial perspective correction factors */
+/* Return non-zero if failed */
+static int
+calc_perspective(
+scanrd_ *s
+) {
+ points *tp;
+ int nl; /* Number of lines used */
+ double ml; /* Minimum length */
+ double pc[4]; /* Perspective factors */
+ double ss[4]; /* Initial search distance */
+ double rv; /* Return value */
+ int rc = 0; /* Return code */
+
+ if (s->novlines < MIN_NO_LINES) {
+ s->errv = SI_FIND_PERSPECTIVE_FAILED;
+ sprintf(s->errm,"Not enough valid lines to compute perspective");
+ return 1;
+ }
+
+ /* Find the longest line */
+ ml = 0.0;
+ tp = s->gdone;
+ FOR_ALL_ITEMS(points, tp)
+ if (tp->flag & F_VALID) {
+ if (tp->len > ml)
+ ml = tp->len;
+ }
+ END_FOR_ALL_ITEMS(tp);
+
+ /* Make minimum line length to be included in angle */
+ /* calculation 1% of longest line */
+ ml *= 0.01;
+
+ /* Mark lines long enough to participate in angle calculation */
+ tp = s->gdone;
+ FOR_ALL_ITEMS(points, tp)
+ if (tp->flag & F_VALID && tp->len >= ml)
+ tp->flag |= F_LONGENOUGH;
+ END_FOR_ALL_ITEMS(tp);
+
+ /* Locate the perspective correction factors that minimze the */
+ /* variance of the mean angle. */
+
+ pc[0] = 0.0;
+ pc[1] = 0.0;
+ pc[2] = 0.5 * s->width;
+ pc[3] = 0.5 * s->height;
+
+ ss[0] = 0.0001;
+ ss[1] = 0.0001;
+ ss[2] = 1.0001;
+ ss[3] = 1.0001;
+ rc = powell(&rv, 4, pc,ss,1e-8,2000,pfunc,s, NULL, NULL);
+
+ if (rc == 0) {
+ points *tp;
+
+ DBG((dbgo,"Perspective correction factors = %f %f %f %f\n",
+ pc[0],pc[1],pc[2],pc[3]));
+
+ s->ppc[0] = pc[0];
+ s->ppc[1] = pc[1];
+ s->ppc[2] = pc[2];
+ s->ppc[3] = pc[3];
+
+ /* Implement the perspective correction */
+ tp = s->gdone;
+ FOR_ALL_ITEMS(points, tp)
+ if (tp->flag & F_LONGENOUGH) {
+ double a, ca;
+ invppersp(s, &tp->x1, &tp->y1, tp->px1, tp->py1, s->ppc);
+ invppersp(s, &tp->x2, &tp->y2, tp->px2, tp->py2, s->ppc);
+ tp->mx = 0.5 * (tp->x2 + tp->x1);
+ tp->my = 0.5 * (tp->y2 + tp->y1);
+ tp->len = sqrt((tp->x2 - tp->x1) * (tp->x2 - tp->x1)
+ + (tp->y2 - tp->y1) * (tp->y2 - tp->y1));
+
+ /* Compute the angle */
+ a = atan2(tp->x2 - tp->x1,tp->y2 - tp->y1);
+
+ /* Make angle +ve */
+ while (a < 0.0)
+ a += M_PI;
+
+ /* Compute the Constrained to 90 degrees angle */
+ /* We use the adivval to figure out where to split angles */
+ /* Split at 0 if adivval == 0.0, split at 45 if adivval == 1.0 */
+ if (a >= (M_PI * (1.0 - s->adivval/4.0)))
+ ca = a - M_PI;
+ else if (a >= (M_PI * (0.5 - s->adivval/4.0)))
+ ca = a - M_PI_2;
+ else
+ ca = a;
+
+ tp->a = a;
+ tp->ca = ca;
+ }
+ END_FOR_ALL_ITEMS(tp);
+ }
+
+ return 0;
+}
+
+/********************************************************************************/
+/* Calculate the image rotation */
+/* Return non-zero if failed */
+static int
+calc_rotation(
+scanrd_ *s
+) {
+ points *tp;
+ int nl; /* Number of lines used */
+ double ml; /* Minimum length */
+ double aa; /* Average angle */
+ double sd,dw; /* Standard deviation, deviation window */
+ double wt; /* Total weighting = sum of line lengths */
+
+ if (s->novlines < MIN_NO_LINES) {
+ s->errv = SI_FIND_ROTATION_FAILED;
+ sprintf(s->errm,"Not enough valid lines to compute rotation angle");
+ return 1;
+ }
+
+ /* Find the longest line */
+ tp = s->gdone;
+ ml = 0.0;
+ FOR_ALL_ITEMS(points, tp)
+ if (tp->flag & F_VALID) {
+ if (tp->len > ml)
+ ml = tp->len;
+ }
+ END_FOR_ALL_ITEMS(tp);
+
+ /* Make minimum line length to be included in angle */
+ /* calculation 1% of longest line */
+ ml *= 0.01;
+
+ /* Calculate the mean angle */
+ aa = 0.0;
+ wt = 0.0; /* Total weighting = sum of line lengths */
+ tp = s->gdone;
+ FOR_ALL_ITEMS(points, tp)
+ if (tp->flag & F_VALID && tp->len >= ml) {
+ aa += tp->len * tp->ca;
+ wt += tp->len;
+ }
+ END_FOR_ALL_ITEMS(tp);
+ aa /= wt;
+
+ if (s->verb >= 2)
+ DBG((dbgo,"Mean angle = %f\n",DEG(aa)));
+
+ /* Calculate the angle standard deviation */
+ tp = s->gdone;
+ sd = 0.0;
+ FOR_ALL_ITEMS(points, tp)
+ if (tp->flag & F_VALID && tp->len >= ml) {
+ double tt;
+ tt = tp->ca - aa;
+ sd += tp->len * tt * tt;
+ }
+ END_FOR_ALL_ITEMS(tp);
+
+ sd = sqrt(sd/wt);
+
+ if (s->verb >= 2)
+ DBG((dbgo,"Standard deviation = %f\n",DEG(sd)));
+
+ /* Now re-compute the angle while rejecting any that fall outside one standard deviation */
+ s->irot = 0.0;
+ wt = 0.0; /* Total weighting = sum of line lengths */
+ nl = 0;
+ dw = sd * SD_WINDOW; /* Allow += 0.5 of a standard deviation */
+ if (dw < 0.01) /* A perfect chart may have dw of zero */
+ dw = 0.01;
+ tp = s->gdone;
+ FOR_ALL_ITEMS(points, tp)
+ if (tp->flag & F_VALID && tp->len >= ml && fabs(tp->ca - aa) <= dw) {
+ s->irot += tp->len * tp->ca;
+ wt += tp->len;
+ nl++;
+ }
+ END_FOR_ALL_ITEMS(tp);
+ if (nl < (MIN_NO_LINES/2)) {
+ s->errv = SI_FIND_ROTATION_FAILED;
+ sprintf(s->errm,"%d consistent lines is not enough to compute rotation angle",nl);
+ return 1;
+ }
+ s->irot /= wt;
+
+ if (s->verb >= 2)
+ DBG((dbgo,"Robust mean angle = %f from %d lines\n",DEG(s->irot),nl));
+
+ return 0;
+}
+
+/********************************************************************************/
+/* Coalesce close entries of an edge list */
+/* return non-zero on error */
+static int
+coalesce_elist(
+scanrd_ *s,
+elist *el,
+int close /* Closeness factor, smaller = coarser */
+) {
+ double r; /* Margin for coalescence */
+ int i,k;
+
+ if (el->c < 2) /* Need at least 2 entries */
+ return 0;
+
+ r = (el->a[el->c-1].pos - el->a[0].pos)/(double)close;
+ for (k = 0, i = 1; i < el->c; i++) {
+ if ((el->a[i].pos - el->a[k].pos) <= r) {
+ /* Merge the two */
+ double lk = el->a[k].len;
+ double li = el->a[i].len;
+ el->a[k].pos = (el->a[k].pos * lk + el->a[i].pos * li)/(lk + li);
+ el->a[k].len = lk + li;
+ if (el->a[k].p1 > el->a[i].p1) /* Track overall start/end points */
+ el->a[k].p1 = el->a[i].p1;
+ if (el->a[k].p2 < el->a[i].p2)
+ el->a[k].p2 = el->a[i].p2;
+ continue;
+ }
+ k++; /* Inc destination pointer */
+ if (k != i)
+ el->a[k] = el->a[i]; /* shuffle data down */
+ }
+ k++; /* one past last out entry */
+ el->c = k;
+ return 0;
+}
+
+static int invert_elist(scanrd_ *s, elist *dl, elist *sl);
+static void debug_elist(scanrd_ *s, elist *el);
+
+/* Make up the x and y edge lists */
+/* Return non-zero if failed */
+static int
+calc_elists(
+scanrd_ *s,
+int ref /* 1 if generating reference lists */
+) {
+ int outw = s->width;
+ int outh = s->height;
+ points *tp;
+ int i,j;
+ double cirot,sirot; /* cos and sin of -irot */
+ elist xl, yl; /* Temporary X and Y edge lists array */
+ elist tl; /* temporary crossing list */
+
+ /* Allocate structures for edge lists */
+ if ((xl.a = (epoint *) malloc(sizeof(epoint) * s->novlines)) == NULL) {
+ s->errv = SI_MALLOC_ELIST;
+ sprintf(s->errm,"scanrd: calc_elist: malloc failed - novlines = %d",s->novlines);
+ return 1;
+ }
+ xl.c = 0;
+ if ((yl.a = (epoint *) malloc(sizeof(epoint) * s->novlines)) == NULL) {
+ s->errv = SI_MALLOC_ELIST;
+ sprintf(s->errm,"scanrd: calc_elist: malloc failed - novlines = %d",s->novlines);
+ return 1;
+ }
+ yl.c = 0;
+
+ /* Put valid lines into one of the two edge list arrays */
+ cirot = cos(-s->irot);
+ sirot = sin(-s->irot);
+ tp = s->gdone;
+ FOR_ALL_ITEMS(points, tp)
+ if (tp->flag & F_VALID) {
+ /* Rotate the point about 0,0 by angle -irot */
+ double x,y,a;
+ double mx = tp->mx, my = tp->my;
+
+ if (ref) { /* Rotate about center of raster for reference generation */
+ mx -= outw/2.0; /* Rotate about center of raster */
+ my -= outh/2.0;
+ x = mx * cirot + my * sirot + outw/2.0;
+ y = -mx * sirot + my * cirot + outh/2.0;
+ } else { /* Rotate about 0,0 for matching */
+ x = mx * cirot + my * sirot;
+ y = -mx * sirot + my * cirot;
+ }
+ a = tp->a - s->irot;
+ if ((a >= -0.08 && a <= 0.08) || (a >= (M_PI-0.08) && a <= (M_PI+0.08))) {
+ xl.a[xl.c].pos = x;
+ xl.a[xl.c].len = tp->len;
+ xl.a[xl.c].p1 = y - tp->len/2.0;
+ xl.a[xl.c].p2 = y + tp->len/2.0;
+ xl.c++;
+ } else if (a >= (M_PI_2-0.08) && a <= (M_PI_2+0.08)) {
+ yl.a[yl.c].pos = y;
+ yl.a[yl.c].len = tp->len;
+ yl.a[yl.c].p1 = x - tp->len/2.0;
+ yl.a[yl.c].p2 = x + tp->len/2.0;
+ yl.c++;
+ }
+ }
+ END_FOR_ALL_ITEMS(tp);
+
+ /* ~~~~ need to check that lists have a reasonable number of entries ~~~~~ */
+
+ /* now sort the lists */
+#define HEAP_COMPARE(A,B) (A.pos < B.pos)
+ HEAPSORT(epoint,xl.a,xl.c);
+ HEAPSORT(epoint,yl.a,yl.c);
+#undef HEAP_COMPARE
+
+ /* Copy the temporary lists to the real lists */
+ if ((s->xelist.a = (epoint *) malloc(sizeof(epoint) * xl.c)) == NULL) {
+ s->errv = SI_MALLOC_ELIST;
+ sprintf(s->errm,"scanrd: calc_elist: malloc failed, xl.c = %d",xl.c);
+ return 1;
+ }
+ s->xelist.c = xl.c;
+ for (i=0; i < xl.c; i++)
+ s->xelist.a[i] = xl.a[i];
+ if ((s->yelist.a = (epoint *) malloc(sizeof(epoint) * yl.c)) == NULL) {
+ s->errv = SI_MALLOC_ELIST;
+ sprintf(s->errm,"scanrd: calc_elist: malloc failed, yl.c = %d",yl.c);
+ return 1;
+ }
+ s->yelist.c = yl.c;
+ for (i=0; i < yl.c; i++)
+ s->yelist.a[i] = yl.a[i];
+
+ /* Coalese close entries of the final lists */
+ if (coalesce_elist(s, &s->xelist,ELISTCDIST))
+ return 1;
+ if (coalesce_elist(s, &s->yelist,ELISTCDIST))
+ return 1;
+
+ /* Calculate crossing count for lines in the X and y lists */
+ if ((tl.a = (epoint *) malloc(sizeof(epoint) * (xl.c > yl.c ? xl.c : yl.c))) == NULL) {
+ s->errv = SI_MALLOC_ELIST;
+ sprintf(s->errm,"scanrd: calc_elist: malloc failed, xl.c = %d, yl.c = %d",xl.c,yl.c);
+ return 1;
+ }
+ /* X list */
+ for (i = 0; i < s->xelist.c; i++) {
+ double ppos = s->xelist.a[i].pos;
+ double pp,np; /* Previous and next pos */
+ if ((i-1) >= 0)
+ pp = (ppos + s->xelist.a[i-1].pos)/2.0; /* Half distance to next line */
+ else
+ pp = -1e6;
+ if ((i+1) < s->xelist.c)
+ np = (ppos + s->xelist.a[i+1].pos)/2.0; /* Half distance to next line */
+ else
+ np = 1e6;
+
+ /* For all the lines in the Y list */
+ for (tl.c = j = 0; j < yl.c; j++) {
+ double pos = yl.a[j].pos;
+ double p1 = yl.a[j].p1;
+ double p2 = yl.a[j].p2;
+ if (p1 <= pp)
+ p1 = pp;
+ if (p2 >= np)
+ p2 = np;
+ /* If crosses on this lines X within +-0.5 of line each side */
+ if (p1 <= np && p2 >= pp) {
+ tl.a[tl.c].pos = pos;
+ tl.a[tl.c].len = p2 - p1;
+ tl.a[tl.c].p1 = p1;
+ tl.a[tl.c].p2 = p2;
+ tl.c++;
+ }
+ }
+ /* now coalesce the crossings */
+ if (coalesce_elist(s,&tl,200))
+ return 1;
+ /* Put count in line we're working on */
+ s->xelist.a[i].ccount = (double)tl.c;
+ pp = ppos;
+ }
+
+ /* Y list */
+ for (i = 0; i < s->yelist.c; i++) {
+ double ppos = s->yelist.a[i].pos;
+ double pp,np; /* Previous and next pos */
+ if ((i-1) >= 0)
+ pp = (ppos + s->yelist.a[i-1].pos)/2.0; /* Half distance to next line */
+ else
+ pp = -1e6;
+ if ((i+1) < s->xelist.c)
+ np = (ppos + s->yelist.a[i+1].pos)/2.0; /* Half distance to next line */
+ else
+ np = 1e6;
+
+ for (tl.c = j = 0; j < xl.c; j++) {
+ double pos = xl.a[j].pos;
+ double p1 = xl.a[j].p1;
+ double p2 = xl.a[j].p2;
+ if (p1 <= pp)
+ p1 = pp;
+ if (p2 >= np)
+ p2 = np;
+ /* If crosses on this lines Y within +-0.5 of line each side */
+ if (p1 <= np && p2 >= pp) {
+ tl.a[tl.c].pos = pos;
+ tl.a[tl.c].len = p2 - p1;
+ tl.a[tl.c].p1 = p1;
+ tl.a[tl.c].p2 = p2;
+ tl.c++;
+ }
+ }
+ /* now coalesce the crossings */
+ if (coalesce_elist(s,&tl,200))
+ return 1;
+ /* Put count in line we're working on */
+ s->yelist.a[i].ccount = (double)tl.c;
+ pp = ppos;
+ }
+
+ /* Normalize the length and ccount */
+ {
+ double tlen; /* Total length maximum */
+ double tcmax; /* Total count maximum */
+ for (tlen = tcmax = 0.0, i=0; i < s->xelist.c; i++) {
+ if (tlen < s->xelist.a[i].len)
+ tlen = s->xelist.a[i].len;
+ if (tcmax < s->xelist.a[i].ccount)
+ tcmax = s->xelist.a[i].ccount;
+ }
+ for (i=0; i < s->xelist.c; i++) {
+ s->xelist.a[i].len /= tlen;
+ s->xelist.a[i].ccount /= tcmax;
+ }
+ for (tlen = tcmax = 0.0, i=0; i < s->yelist.c; i++) {
+ if (tlen < s->yelist.a[i].len)
+ tlen = s->yelist.a[i].len;
+ if (tcmax < s->yelist.a[i].ccount)
+ tcmax = s->yelist.a[i].ccount;
+ }
+ for (i=0; i < s->yelist.c; i++) {
+ s->yelist.a[i].len /= tlen;
+ s->yelist.a[i].ccount /= tcmax;
+ }
+ }
+
+ /* Create the inverted lists for any rotation matching */
+ if (invert_elist(s, &s->ixelist, &s->xelist))
+ return 1;
+ if (invert_elist(s, &s->iyelist, &s->yelist))
+ return 1;
+
+ if (s->verb >= 3) {
+ DBG((dbgo,"\nxelist:\n"));
+ debug_elist(s,&s->xelist);
+ DBG((dbgo,"\nixelist:\n"));
+ debug_elist(s,&s->ixelist);
+ DBG((dbgo,"\nyelist:\n"));
+ debug_elist(s,&s->yelist);
+ DBG((dbgo,"\niyelist:\n"));
+ debug_elist(s,&s->iyelist);
+ }
+
+ /* Clean up */
+ free(xl.a);
+ free(yl.a);
+ free(tl.a);
+ return 0;
+}
+
+/********************************************************************************/
+/* Write the elists out to a file */
+
+/* Increment a string counter */
+static void
+strinc(
+char *s
+) {
+ int i,n,c; /* Length of string and carry flag */
+ n = strlen(s);
+ for (c = 1, i = n-1; i >= 0 && c != 0; i--) {
+ char sval = ' ';
+ if (s[i] == '9') {
+ s[i] = '0';
+ sval = '1';
+ c = 1;
+ } else if (s[i] == 'z') {
+ s[i] = 'a';
+ sval = 'a';
+ c = 1;
+ } else if (s[i] == 'Z') {
+ s[i] = 'A';
+ sval = 'A';
+ c = 1;
+ } else {
+ s[i]++;
+ c = 0;
+ }
+ if (i == 0 && c != 0) {
+ /* Assume there is some more space */
+ for (i = n; i >= 0; i--)
+ s[i+1] = s[i];
+ s[0] = sval;
+ break;
+ }
+ }
+}
+
+/* Write out the match reference information */
+/* Return non-zero on error */
+static int
+write_elists(
+scanrd_ *s
+) {
+ char *fname = s->refname; /* Path of file to write to */
+ FILE *elf;
+ int i;
+
+ if ((elf=fopen(fname,"w"))==NULL) {
+ s->errv = SI_REF_WRITE_ERR;
+ sprintf(s->errm,"write_elists: error opening match reference file '%s'",fname);
+ return 1;
+ }
+
+ fprintf(elf,"REF_ROTATION %f\n\n",DEG(s->irot));
+
+ fprintf(elf,"XLIST %d\n",s->xelist.c);
+ for (i = 0; i < s->xelist.c; i++)
+ fprintf(elf," %f %f %f\n",s->xelist.a[i].pos, s->xelist.a[i].len, s->xelist.a[i].ccount);
+ fprintf(elf,"\n");
+
+ fprintf(elf,"YLIST %d\n",s->yelist.c);
+ for (i = 0; i < s->yelist.c; i++)
+ fprintf(elf," %f %f %f\n",s->yelist.a[i].pos, s->yelist.a[i].len, s->yelist.a[i].ccount);
+ fprintf(elf,"\n");
+
+ if ((fclose(elf)) == EOF) {
+ s->errv = SI_REF_WRITE_ERR;
+ error("write_elists: Unable to close match reference file '%s'\n",fname);
+ return 1;
+ }
+ return 0;
+}
+
+/* Read in an elist reference file */
+/* return non-zero on error */
+/* (~~~ the line counting is rather broken ~~~) */
+static int
+read_relists(
+scanrd_ *s
+) {
+ char *fname = s->refname; /* Path of file to read from */
+ FILE *elf;
+ int i,l = 1;
+ int rv;
+ char *em; /* Read error message */
+
+ if ((elf=fopen(fname,"r"))==NULL) {
+ s->errv = SI_REF_READ_ERR;
+ sprintf(s->errm,"read_elists: error opening match reference file '%s'",fname);
+ return 1;
+ }
+
+ s->fid[0] = s->fid[1] = 0.0;
+ s->fid[2] = s->fid[3] = 0.0;
+ s->fid[4] = s->fid[5] = 0.0;
+ s->fid[6] = s->fid[7] = 0.0;
+
+ /* BOXES */
+ for(;;) {
+ if((rv = fscanf(elf,"BOXES %d",&s->nsbox)) == 1) {
+ l++;
+ break;
+ }
+ if (rv == EOF) {
+ em = "Didn't find BOXES before end of file";
+ goto read_error;
+ }
+ if (rv == 0) {
+ while ((rv = getc(elf)) != '\n' && rv != EOF);
+ l++;
+ }
+ }
+
+ /* Allocate structures for boxes */
+ if ((s->sboxes = (sbox *) calloc(s->nsbox, sizeof(sbox))) == NULL) {
+ s->errv = SI_MALLOC_REFREAD;
+ sprintf(s->errm,"read_elist, malloc failed");
+ return 1;
+ }
+ for (i = 0; i < s->nsbox;) {
+ char xfix1[20], xfix2[20], yfix1[20],yfix2[20];
+ char xfirst[20];
+ double ox,oy,w,h,xi,yi;
+ char xf[20];
+ double x;
+
+ if(fscanf(elf," %19s %19s %19s %19s %19s %lf %lf %lf %lf %lf %lf",xfirst ,xfix1, xfix2, yfix1, yfix2, &w, &h, &ox, &oy, &xi, &yi) != 11) {
+ em = "Read of BOX failed";
+ goto read_error;
+ }
+ l++;
+
+ /* If Fiducial. Typically top left, top right, botton right, bottom left. */
+ if (xfirst[0] == 'F') {
+ s->fid[0] = atof(yfix1);
+ s->fid[1] = atof(yfix2);
+ s->fid[2] = w;
+ s->fid[3] = h;
+ s->fid[4] = ox;
+ s->fid[5] = oy;
+ s->fid[6] = xi;
+ s->fid[7] = yi;
+ s->fidsize = fabs(s->fid[2] - s->fid[0]) + fabs(s->fid[5] - s->fid[3]);
+ s->fidsize /= 80.0;
+ s->havefids = 1;
+
+//printf("~1 fiducials %f %f, %f %f %f, %f\n",w, h, ox,oy, xi, yi);
+ continue;
+ }
+ for(;;) { /* Do Y increment */
+ x = ox;
+ strcpy(xf,xfix1);
+ for(;;) { /* Do X increment */
+ if (i >= s->nsbox) {
+ em = "More BOXes that declared";
+ goto read_error;
+ }
+ /* '_' is used as a null string marker for single character single cells */
+ if (xf[0] == '_')
+ sprintf(s->sboxes[i].name,"%s",yfix1);
+ else if (yfix1[0] == '_')
+ sprintf(s->sboxes[i].name,"%s",xf);
+ else { /* Y indicates Y name comes first */
+ if (xfirst[0] == 'Y')
+ sprintf(s->sboxes[i].name,"%s%s",yfix1,xf);
+ else /* X or D */
+ sprintf(s->sboxes[i].name,"%s%s",xf,yfix1);
+ }
+ if (xfirst[0] == 'D')
+ s->sboxes[i].diag = 1; /* Diagnostic box - don't print name or read pixels */
+ else
+ s->sboxes[i].diag = 0;
+ s->sboxes[i].x1 = x;
+ s->sboxes[i].y1 = oy;
+ s->sboxes[i].x2 = x + w;
+ s->sboxes[i].y2 = oy + h;
+
+ /* Misc. init. of new sbox */
+ s->sboxes[i].xpt[0] = -1.0; /* No default expected value */
+
+ i++;
+ x += xi;
+ if (strcmp(xf,xfix2) == 0)
+ break;
+ strinc(xf);
+ }
+ if (strcmp(yfix1,yfix2) == 0)
+ break;
+ oy += yi;
+ strinc(yfix1);
+ }
+ }
+
+ /* BOX_SHRINK */
+ for(;;) {
+ if((rv = fscanf(elf,"BOX_SHRINK %lf ",&s->rbox_shrink)) == 1) {
+ l++;
+ break;
+ }
+ if (rv == EOF) {
+ em = "Didn't find BOX_SHRINK before end of file";
+ goto read_error;
+ }
+ if (rv == 0) {
+ while ((rv = getc(elf)) != '\n' && rv != EOF);
+ l++;
+ }
+ }
+
+ /* XLIST */
+ for(;;) {
+ if((rv = fscanf(elf,"XLIST %d ",&s->rxelist.c)) == 1) {
+ l++;
+ break;
+ }
+ if (rv == EOF) {
+ em = "Didn't find XLIST before end of file";
+ goto read_error;
+ }
+ if (rv == 0) {
+ while ((rv = getc(elf)) != '\n' && rv != EOF);
+ l++;
+ }
+ }
+ /* Allocate structures for ref edge lists */
+ if ((s->rxelist.a = (epoint *) malloc(sizeof(epoint) * s->rxelist.c)) == NULL) {
+ s->errv = SI_MALLOC_REFREAD;
+ sprintf(s->errm,"read_elist, malloc failed");
+ return 1;
+ }
+ for (i = 0; i < s->rxelist.c; i++) {
+ if (fscanf(elf," %lf %lf %lf ",
+ &s->rxelist.a[i].pos, &s->rxelist.a[i].len, &s->rxelist.a[i].ccount) != 3) {
+ em = "Failed to read an XLIST line";
+ goto read_error;
+ }
+ l++;
+ }
+
+ /* YLIST */
+ for(;;) {
+ if ((rv = fscanf(elf,"YLIST %d ",&s->ryelist.c)) == 1) {
+ l++;
+ break;
+ }
+ if (rv == EOF) {
+ em = "Didn't find YLIST before end of file";
+ goto read_error;
+ }
+ if (rv == 0) {
+ while ((rv = getc(elf)) != '\n' && rv != EOF);
+ l++;
+ }
+ }
+ if ((s->ryelist.a = (epoint *) malloc(sizeof(epoint) * s->ryelist.c)) == NULL) {
+ s->errv = SI_MALLOC_REFREAD;
+ sprintf(s->errm,"read_elist, malloc failed");
+ return 1;
+ }
+ for (i = 0; i < s->ryelist.c; i++) {
+ if (fscanf(elf," %lf %lf %lf ",
+ &s->ryelist.a[i].pos, &s->ryelist.a[i].len, &s->ryelist.a[i].ccount) != 3)
+ {
+ em = "Failed to read an YLIST line";
+ goto read_error;
+ }
+ l++;
+ }
+
+ /* EXPECTED */
+ {
+ int j;
+ int isxyz = 0;
+ int nxpt = 0;
+ char csps[20];
+
+ for(;;) {
+ if ((rv = fscanf(elf,"EXPECTED %19s %d ",csps, &nxpt)) == 2) {
+ l++;
+ if (strcmp(csps, "XYZ") == 0) {
+ isxyz = 1;
+ break;
+ } else if (strcmp(csps, "LAB") == 0) {
+ isxyz = 0;
+ break;
+ } else {
+ em = "Unknown EXPECTED colorespace";
+ goto read_error;
+ }
+ }
+ if (rv == EOF) {
+ break;
+ }
+ if (rv == 0) {
+ while ((rv = getc(elf)) != '\n' && rv != EOF);
+ l++;
+ }
+ }
+ for (j = 0; j < nxpt; j++) {
+ char name[20];
+ double val[3];
+ if (fscanf(elf," %19s %lf %lf %lf ",
+ name, &val[0], &val[1], &val[2]) != 4)
+ {
+ em = "Failed to read an EXPECTED line";
+ goto read_error;
+ }
+ l++;
+ /* Now locate the matching box */
+ for (i = 0; i < s->nsbox; i++) {
+ if (strcmp(s->sboxes[i].name, name) == 0) { /* Found it */
+ if (isxyz) {
+ XYZ2Lab(s->sboxes[i].xpt, val);
+ } else {
+ s->sboxes[i].xpt[0] = val[0];
+ s->sboxes[i].xpt[1] = val[1];
+ s->sboxes[i].xpt[2] = val[2];
+ }
+ s->xpt = 1;
+ break;
+ }
+ }
+ if (i >= s->nsbox) {
+ em = "Failed to locate matching sample box in EXPECTED list";
+ goto read_error;
+ }
+ }
+ }
+
+ if ((fclose(elf)) == EOF) {
+ s->errv = SI_REF_WRITE_ERR;
+ error("read_elists: Unable to close match reference file '%s'\n",fname);
+ return 1;
+ }
+
+ /* Generate length normalization factor */
+ {
+ double tlen; /* Total of normalized length */
+ for (tlen = 0.0, i=0; i < s->rxelist.c; i++)
+ tlen += s->rxelist.a[i].len;
+ s->rxelist.lennorm = tlen;
+ for (tlen = 0.0, i=0; i < s->ryelist.c; i++)
+ tlen += s->ryelist.a[i].len;
+ s->ryelist.lennorm = tlen;
+ }
+
+ if (s->verb >= 3) {
+ DBG((dbgo,"\nrxelist:\n"));
+ debug_elist(s, &s->rxelist);
+ DBG((dbgo,"\nryelist:\n"));
+ debug_elist(s, &s->ryelist);
+ }
+
+ return 0;
+
+read_error:;
+ s->errv = SI_REF_FORMAT_ERR;
+ sprintf(s->errm,"read_relist failed at line %d in file %s: %s\n",l,fname,em);
+ return 1;
+}
+
+/********************************************************************************/
+/* Create an inverted direction elist */
+/* return non-zero on error */
+static int
+invert_elist(
+scanrd_ *s,
+elist *dl, /* Destination list */
+elist *sl /* Source list */
+) {
+ int i,j, rc = sl->c;
+
+ *dl = *sl; /* Copy all the structure elements */
+
+ /* Allocate space in the destination list */
+ if ((dl->a = (epoint *) malloc(sizeof(epoint) * rc)) == NULL) {
+ s->errv = SI_MALLOC_ELIST;
+ sprintf(s->errm,"invert_elist: malloc failed");
+ return 1;
+ }
+
+ /* Copy the array data and reverse its order */
+ for (i = 0, j = rc-1; i < rc; i++,j--) {
+ dl->a[j] = sl->a[i]; /* Copy array element */
+ dl->a[j].pos = -dl->a[j].pos; /* Invert position */
+ }
+ return 0;
+}
+
+/* Print out elist */
+static void
+debug_elist(
+scanrd_ *s,
+elist *el
+) {
+ int i, rc = el->c;
+
+ DBG((dbgo,"Elist has %d entries allocated at 0x%p\n",el->c,el->a));
+ DBG((dbgo,"lennorm = %f\n",el->lennorm));
+ for (i = 0; i < rc; i++)
+ DBG((dbgo," [%d] = %f %f %f\n",i,el->a[i].pos,el->a[i].len,el->a[i].ccount));
+}
+
+/* Free the array data in an elist */
+static void
+free_elist_array(elist *el) {
+ free(el->a);
+ el->c = 0;
+}
+
+/********************************************************************************/
+/* !!!!!!! */
+/* NEED TO RESOLVE WHY current code is better in some cases, but */
+/* not in others. */
+
+#ifndef NEVER /* Current code */
+
+/* Compute a correlation between two elists */
+static double
+elist_correl(
+scanrd_ *s,
+elist *r, /* Reference list */
+elist *t, /* Target list */
+double off, double scale, /* Offset and scale of target to ref */
+int verb /* Verbose mode */
+) {
+ int i, j, rc = r->c;
+ double cc = 0.0; /* Correlation */
+ double marg = (r->a[rc-1].pos - r->a[0].pos)/150.0; /* determines sharpness of pos. match */
+ double marg2 = marg * 3.0; /* Don't contribute anything outside this distance */
+
+ for (i = j = 0; i < t->c; i++) {
+ int ri; /* Reference index */
+ double dd,d1,d2; /* Distance to nearest reference */
+ double pos = (t->a[i].pos + off) * scale;
+ double len = t->a[i].len;
+ double cnt = t->a[i].ccount;
+ while (pos > r->a[j+1].pos && j < (r->c-2)) j++;
+ d1 = fabs(pos - r->a[j].pos);
+ d2 = fabs(r->a[j+1].pos - pos);
+ if (d1 < d2) {
+ dd = d1;
+ ri = j;
+ } else {
+ dd = d2;
+ ri = j+1;
+ }
+ if (dd <= marg2) { /* If close enough to reference */
+ double ccf, rcnt = r->a[ri].ccount;
+ double llf, rlen = r->a[ri].len;
+ double df = marg/(marg + dd);
+ df *= df;
+ ccf = 1.0 - (rcnt > cnt ? rcnt-cnt : cnt-rcnt);
+ llf = 1.0 - (rlen > len ? rlen-len : len-rlen);
+ /* The weighting gives slightly more emphasis on matching long lines */
+ cc += (1.0 + rlen) * (df * llf * ccf);
+ if (verb) {
+ DBG((dbgo,"---- t[%d] %f %f %f this cc = %f, running total cc = %f\n r[%d] %f %f %f, df = %f, llf = %f, ccf = %f\n",
+ i,pos,len,cnt,df * llf * ccf,cc,j,r->a[ri].pos,r->a[ri].len,rcnt,df, llf, ccf));
+ }
+ }
+ }
+ return cc/(r->lennorm + (double)r->c); /* Normalize */
+}
+
+#else /* New test code */
+
+/* Compute a correlation between two elists */
+static double
+elist_correl(
+scanrd_ *s,
+elist *r, /* Reference list */
+elist *t, /* Target list */
+double off, double scale, /* Offset and scale of target to ref */
+int verb /* Verbose mode */
+) {
+ int i, rc = r->c;
+ double cc = 0.0; /* Correlation */
+ double marg = (r->a[rc-1].pos - r->a[0].pos)/100.0; /* determines sharpness of pos. match */
+ double marg2 = marg * marg; /* marg squared */
+
+//printf("~1 doing elist_correl\n");
+ /* For each reference edge */
+ for (i = 0; i < rc; i++) {
+ int j[3], jj, bj, tc = t->c;
+ double dd, pos, bdd;
+
+ /* Find the closest target edge using binary search. */
+ for(bdd = 1e6, j[2] = tc-1, j[0] = 0; j[2] > (j[0]+1);) {
+ double dist;
+ j[1] = (j[2] + j[0])/2; /* Trial point */
+ dist = r->a[i].pos - (t->a[j[1]].pos + off) * scale;
+
+//printf("~1 j1 = %d, j1 = %d, j0 = %d, dist = %f\n",j[2], j[1], j[0], dist);
+ if (dist > 0) {
+ j[0] = j[1];
+ } else {
+ j[2] = j[1];
+ }
+ }
+
+ /* Locate best out of 3 remaining points */
+ for (jj = 0; jj < 3; jj++) {
+ double dist;
+ pos = (t->a[j[jj]].pos + off) * scale;
+ dist = r->a[i].pos - pos;
+ dd = dist * dist; /* Distance squared */
+ if (dd < bdd) { /* New closest */
+ bdd = dd;
+ bj = j[jj];
+ }
+ }
+
+//printf("~1 best j = %d, bdd = %f, marg2 = %f\n",bj,bdd,marg2);
+ /* Compute correlation */
+ if (bdd < marg2) { /* Within our margine */
+ double df = (marg2 - bdd)/marg2; /* Distance factor */
+ double llf, rlen = r->a[i].len, len = t->a[i].len;
+ double ccf, rcnt = r->a[i].ccount, cnt = t->a[i].ccount;
+ double tcc;
+ llf = 1.0 - (rlen > len ? rlen-len : len-rlen);
+ ccf = 1.0 - (rcnt > cnt ? rcnt-cnt : cnt-rcnt);
+
+ /* The weighting gives slightly more emphasis on matching long lines */
+ /* Not using crossing count */
+ tcc = (1.0 + rlen) * (df * llf);
+ cc += tcc;
+ if (verb) {
+ DBG((dbgo,"---- targ[%d] %f %f %f this cc = %f, running total cc = %f\n",
+ bj,pos,t->a[bj].len,t->a[bj].ccount, tcc,cc));
+ DBG((dbgo," ref[%d] %f %f %f, df = %f, llf = %f, ccf = %f\n",
+ i,r->a[i].pos,r->a[i].len,r->a[i].ccount, df, llf, ccf));
+ }
+ }
+
+ }
+ return cc/(r->lennorm + (double)r->c); /* Normalize */
+}
+
+#endif /* NEVER */
+
+/* Structure to hold data for optimization function */
+struct _edatas {
+ scanrd_ *s; /* scanrd object */
+ elist *r; /* Reference list */
+ elist *t; /* Target list */
+ int verb; /* Verbose mode */
+ }; typedef struct _edatas edatas;
+
+/* Definition of the optimization function handed to powell() */
+static double
+efunc(void *edata, double p[]) {
+ edatas *e = (edatas *)edata;
+ double rv = 2.0 - elist_correl(e->s,e->r,e->t,p[0],p[1],e->verb);
+ return rv;
+}
+
+/* return non-zero on error */
+static int
+best_match(
+scanrd_ *s,
+elist *r, /* Reference list */
+elist *t, /* Target list */
+ematch *rv /* Return values */
+) {
+ int r0,r1,rw,t0,t1;
+ double rwidth;
+ double cc;
+ double bcc = 0.0, boff = 0.0, bscale = 0.0; /* best values */
+
+ /* The target has been rotated, and we go through all reasonable */
+ /* translations and scales to see if we can match it to the */
+ /* reference. */
+ r0 = 0;
+ r1 = r->c-1;
+ rw = r->c/2; /* Minimum number of target line to match all of reference */
+ if (t->c/2 < rw)
+ rw = t->c/2;
+ rwidth = r->a[r1].pos - r->a[r0].pos;
+
+ for (t0 = 0; t0 < t->c-1; t0++) {
+ double off;
+ for (t1 = t->c-1; t1 > (t0+rw); t1--) {
+ double scale;
+
+ scale = rwidth/(t->a[t1].pos - t->a[t0].pos);
+ if (scale < 0.001 || scale > 100.0) {
+ break; /* Don't bother with silly scale factors */
+ }
+
+ /* Have to compenate the offset for the scale since it is scaled from 0 */
+ off = r->a[r0].pos/scale - t->a[t0].pos;
+ cc = elist_correl(s,r,t,off,scale,0);
+
+ if (s->verb >= 7) {
+ DBG((dbgo,"Matching target [%d]-[%d] to ref [%d]-[%d] = %f-%f to %f-%f\n",
+ t0,t1,r0,r1,t->a[t0].pos,t->a[t1].pos,r->a[r0].pos,r->a[r1].pos));
+ DBG((dbgo,"Initial off %f, scale %f, cc = %f\n",off,scale,cc));
+ }
+ if (cc > 0.20) { /* Looks promising, try optimizing solution */
+ double cp[2]; /* Start point/improved point */
+ double rv; /* Return value */
+ int rc; /* Return code */
+ edatas dd; /* Data structure */
+ double ss[2] = { 0.1, 0.1}; /* Initial search distance */
+
+ dd.s = s; /* scanrd object */
+ dd.r = r; /* Reference list */
+ dd.t = t; /* Target list */
+ dd.verb = 0; /* Verbose mode */
+
+ /* Set search start point */
+ cp[0] = off;
+ cp[1] = scale;
+ /* Set search distance */
+ ss[0] = (0.01 * rwidth/ELISTCDIST)/scale; /* Search distance */
+ ss[1] = scale * 0.01 * rwidth/ELISTCDIST;
+
+ /* Find minimum */
+ rc = powell(&rv, 2,cp,ss,0.0001,400,efunc,&dd, NULL, NULL);
+
+ if (rc == 0 /* Powell converged */
+ && cp[1] > 0.001 && cp[1] < 100.0) { /* and not ridiculous */
+ cc = 2.0 - rv;
+ off = cp[0];
+ scale = cp[1];
+ }
+ /* Else use unoptimsed values */
+
+ if (s->verb >= 7) {
+ DBG((dbgo,"After optimizing, off %f, scale %f, cc = %f\n",off,scale,cc));
+ }
+ }
+
+ if (s->verb >= 7) {
+ if (cc > 0.25) {
+ DBG((dbgo,"Good correlation::\n"));
+ elist_correl(s,r,t,off,scale,1);
+ }
+ }
+ if (s->verb >= 7)
+ DBG((dbgo,"offset %f, scale %f cc %f\n", off,scale,cc));
+ if (cc > 0.0 && cc > bcc) { /* Keep best */
+ boff = off;
+ bscale = scale;
+ bcc = cc;
+ if (s->verb >= 7)
+ DBG((dbgo,"(New best)\n"));
+ }
+ }
+ }
+ if (s->verb >= 7)
+ DBG((dbgo,"Returning best offset %f, scale %f returns %f\n\n", boff,bscale,bcc));
+
+ /* return best values */
+ rv->cc = bcc;
+ rv->off = boff;
+ rv->scale = bscale;
+ return 0;
+}
+
+/* Find best offset and scale match between reference and target, */
+/* and then from this, compute condidate 90 degree rotations. */
+/* Return 0 if got at least one candidate rotation */
+/* Return 1 if no reasonable candidate rotation found */
+/* Return 2 if some other error */
+static int
+do_match(
+scanrd_ *s
+) {
+ ematch xx, yy, xy, yx, xix, yiy, xiy, yix; /* All 8 matches needed to detect rotations */
+ double r0, r90, r180, r270; /* Correlation for each extra rotation of target */
+
+ /* Check out all the matches */
+ if (s->verb >= 2) DBG((dbgo,"Checking xx\n"));
+ if (best_match(s, &s->rxelist,&s->xelist,&xx))
+ return 2;
+ if (s->verb >= 2) DBG((dbgo,"Checking yy\n"));
+ if (best_match(s, &s->ryelist,&s->yelist,&yy))
+ return 2;
+ if (s->verb >= 2) DBG((dbgo,"Checking xy\n"));
+ if (best_match(s, &s->rxelist,&s->yelist,&xy))
+ return 2;
+ if (s->verb >= 2) DBG((dbgo,"Checking yx\n"));
+ if (best_match(s, &s->ryelist,&s->xelist,&yx))
+ return 2;
+ if (s->verb >= 2) DBG((dbgo,"Checking xix\n"));
+ if (best_match(s, &s->rxelist,&s->ixelist,&xix))
+ return 2;
+ if (s->verb >= 2) DBG((dbgo,"Checking yiy\n"));
+ if (best_match(s, &s->ryelist,&s->iyelist,&yiy))
+ return 2;
+ if (s->verb >= 2) DBG((dbgo,"Checking xiy\n"));
+ if (best_match(s, &s->rxelist,&s->iyelist,&xiy))
+ return 2;
+ if (s->verb >= 2) DBG((dbgo,"Checking yix\n"));
+ if (best_match(s, &s->ryelist,&s->ixelist,&yix))
+ return 2;
+
+ if (s->verb >= 2) {
+ DBG((dbgo,"Axis matches for each possible orientation:\n"));
+ DBG((dbgo," 0: xx = %f, yy = %f, xx.sc = %f, yy.sc = %f\n",
+ xx.cc,yy.cc,xx.scale,yy.scale));
+ DBG((dbgo," 90: xiy = %f, yx = %f, xiy.sc = %f, yx.sc = %f\n",
+ xiy.cc,yx.cc,xiy.scale,yx.scale));
+ DBG((dbgo,"180: xix = %f, yiy = %f, xix.sc = %f, yiy.sc = %f\n",
+ xix.cc,yiy.cc,xix.scale,yiy.scale));
+ DBG((dbgo,"270: xy = %f, yix = %f, xy.sc = %f, yix.sc = %f\n",
+ xy.cc,yix.cc,xy.scale,yix.scale));
+ }
+
+ /* Compute the combined values for the four orientations. */
+ /* add penalty for different scale factors */
+ r0 = sqrt(xx.cc * xx.cc + yy.cc * yy.cc)
+ * (xx.scale > yy.scale ? yy.scale/xx.scale : xx.scale/yy.scale);
+ r90 = sqrt(xiy.cc * xiy.cc + yx.cc * yx.cc)
+ * (xiy.scale > yx.scale ? yx.scale/xiy.scale : xiy.scale/yx.scale);
+ r180 = sqrt(xix.cc * xix.cc + yiy.cc * yiy.cc)
+ * (xix.scale > yiy.scale ? yiy.scale/xix.scale : xix.scale/yiy.scale);
+ r270 = sqrt(xy.cc * xy.cc + yix.cc * yix.cc)
+ * (xy.scale > yix.scale ? yix.scale/xy.scale : xy.scale/yix.scale);
+
+ if (s->verb >= 2)
+ DBG((dbgo,"r0 = %f, r90 = %f, r180 = %f, r270 = %f\n",r0,r90,r180,r270));
+
+ s->norots = 0;
+ if (s->flags & SI_GENERAL_ROT) { /* If general rotation allowed */
+ if (s->xpt == 0) { /* No expected color information to check rotations agaist */
+ /* so choose the single best rotation by the edge matching */
+ DBG((dbgo,"There is no expected color information, so best fit rotations will be used\n"));
+ if (r0 >= MATCHCC && r0 >= r90 && r0 >= r180 && r0 >= r270) {
+ s->rots[0].ixoff = -xx.off;
+ s->rots[0].ixscale = 1.0/xx.scale;
+ s->rots[0].iyoff = -yy.off;
+ s->rots[0].iyscale = 1.0/yy.scale;
+ s->rots[0].irot = s->irot;
+ s->rots[0].cc = r0;
+ s->norots = 1;
+ } else if (r90 >= MATCHCC && r90 >= r180 && r90 >= r270) {
+ s->rots[0].ixoff = -xiy.off;
+ s->rots[0].ixscale = 1.0/xiy.scale;
+ s->rots[0].iyoff = -yx.off;
+ s->rots[0].iyscale = 1.0/yx.scale;
+ s->rots[0].irot = s->irot + M_PI_2;
+ s->rots[0].cc = r90;
+ s->norots = 1;
+ } else if (r180 >= MATCHCC && r180 >= r270) {
+ s->rots[0].ixoff = -xix.off;
+ s->rots[0].ixscale = 1.0/xix.scale;
+ s->rots[0].iyoff = -yiy.off;
+ s->rots[0].iyscale = 1.0/yiy.scale;
+ s->rots[0].irot = s->irot + M_PI;
+ s->rots[0].cc = r180;
+ s->norots = 1;
+ } else if (r270 >= MATCHCC) { /* 270 extra target rotation */
+ s->rots[0].ixoff = -xy.off;
+ s->rots[0].ixscale = 1.0/xy.scale;
+ s->rots[0].iyoff = -yix.off;
+ s->rots[0].iyscale = 1.0/yix.scale;
+ s->rots[0].irot = s->irot + M_PI + M_PI_2;
+ s->rots[0].cc = r270;
+ s->norots = 1;
+ }
+
+ } else { /* Got expected color info, so try reasonable rotations */
+ double bcc; /* Best correlation coeff */
+
+ if (r0 >= r90 && r0 >= r180 && r0 >= r270)
+ bcc = r0;
+ else if (r90 >= r180 && r90 >= r270)
+ bcc = r90;
+ else if (r180 >= r270)
+ bcc = r180;
+ else
+ bcc = r270;
+
+ bcc *= ALT_ROT_TH; /* Threshold for allowing alternate rotation */
+ if (bcc < MATCHCC)
+ bcc = MATCHCC;
+
+ s->norots = 0;
+ if (r0 >= bcc) {
+ s->rots[s->norots].ixoff = -xx.off;
+ s->rots[s->norots].ixscale = 1.0/xx.scale;
+ s->rots[s->norots].iyoff = -yy.off;
+ s->rots[s->norots].iyscale = 1.0/yy.scale;
+ s->rots[s->norots].irot = s->irot;
+ s->rots[s->norots].cc = r0;
+ s->norots++;
+ }
+ if (r90 >= bcc) {
+ s->rots[s->norots].ixoff = -xiy.off;
+ s->rots[s->norots].ixscale = 1.0/xiy.scale;
+ s->rots[s->norots].iyoff = -yx.off;
+ s->rots[s->norots].iyscale = 1.0/yx.scale;
+ s->rots[s->norots].irot = s->irot + M_PI_2;
+ s->rots[s->norots].cc = r90;
+ s->norots++;
+ }
+ if (r180 >= bcc) {
+ s->rots[s->norots].ixoff = -xix.off;
+ s->rots[s->norots].ixscale = 1.0/xix.scale;
+ s->rots[s->norots].iyoff = -yiy.off;
+ s->rots[s->norots].iyscale = 1.0/yiy.scale;
+ s->rots[s->norots].irot = s->irot + M_PI;
+ s->rots[s->norots].cc = r180;
+ s->norots++;
+ }
+ if (r270 >= bcc) {
+ s->rots[s->norots].ixoff = -xy.off;
+ s->rots[s->norots].ixscale = 1.0/xy.scale;
+ s->rots[s->norots].iyoff = -yix.off;
+ s->rots[s->norots].iyscale = 1.0/yix.scale;
+ s->rots[s->norots].irot = s->irot + M_PI + M_PI_2;
+ s->rots[s->norots].cc = r270;
+ s->norots++;
+ }
+ }
+ } else { /* Use only rotation 0 */
+ if (r0 >= MATCHCC) {
+ s->rots[0].ixoff = -xx.off;
+ s->rots[0].ixscale = 1.0/xx.scale;
+ s->rots[0].iyoff = -yy.off;
+ s->rots[0].iyscale = 1.0/yy.scale;
+ s->rots[0].irot = s->irot;
+ s->rots[0].cc = r0;
+ s->norots = 1;
+ } else if (s->flags & SI_ASISIFFAIL) {
+ DBG((dbgo, "Recognition failed, reading patches 'as is' (probably incorrect)\n"));
+ s->rots[0].ixoff = 0.0;
+ s->rots[0].ixscale = 1.0;
+ s->rots[0].iyoff = 0.0;
+ s->rots[0].iyscale = 1.0;
+ s->rots[0].irot = 0.0;
+ s->rots[0].cc = r0;
+ s->norots = 1;
+ }
+ }
+
+ if (s->verb >= 2) {
+ int i;
+ DBG((dbgo,"There are %d candidate rotations:\n",s->norots));
+
+ for (i = 0; i < s->norots; i++) {
+ DBG((dbgo,"cc = %f, irot = %f, xoff = %f, yoff = %f, xscale = %f, yscale = %f\n",
+ s->rots[i].cc, DEG(s->rots[i].irot), s->rots[i].ixoff,s->rots[i].iyoff,s->rots[i].ixscale,s->rots[i].iyscale));
+ }
+ }
+
+ if (s->norots == 0)
+ return 1;
+
+ return 0;
+}
+
+/********************************************************************************/
+/* perspective transformation. */
+/* Transform from raster to reference using iptrans[]. */
+/* Transform from reference to raster using ptrans[]. */
+static void ptrans(double *xx, double *yy, double x, double y, double *ptrans) {
+ double den;
+
+ den = ptrans[6] * x + ptrans[7] * y + 1.0;
+
+ if (fabs(den) < 1e-6) {
+ if (den < 0.0)
+ den = -1e-6;
+ else
+ den = 1e-6;
+ }
+
+ *xx = (ptrans[0] * x + ptrans[1] * y + ptrans[2])/den;
+ *yy = (ptrans[3] * x + ptrans[4] * y + ptrans[5])/den;
+}
+
+/* Convert perspective transfom parameters to inverse */
+/* perspective transform parameters. */
+/* Return nz on error */
+int invert_ptrans(double *iptrans, double *ptrans) {
+ double scale = ptrans[0] * ptrans[4] - ptrans[1] * ptrans[3];
+
+ if (fabs(scale) < 1e-6)
+ return 1;
+
+ scale = 1.0/scale;
+
+ iptrans[0] = scale * (ptrans[4] - ptrans[5] * ptrans[7]);
+ iptrans[1] = scale * (ptrans[2] * ptrans[7] - ptrans[1]);
+ iptrans[2] = scale * (ptrans[1] * ptrans[5] - ptrans[2] * ptrans[4]);
+
+ iptrans[3] = scale * (ptrans[5] * ptrans[6] - ptrans[3]);
+ iptrans[4] = scale * (ptrans[0] - ptrans[2] * ptrans[6]);
+ iptrans[5] = scale * (ptrans[2] * ptrans[3] - ptrans[0] * ptrans[5]);
+
+ iptrans[6] = scale * (ptrans[3] * ptrans[7] - ptrans[4] * ptrans[6]);
+ iptrans[7] = scale * (ptrans[1] * ptrans[6] - ptrans[0] * ptrans[7]);
+
+ return 0;
+}
+
+
+/* Structure to hold data for optimization function */
+struct _pdatas {
+ scanrd_ *s; /* scanrd object */
+ double *tar; /* 4 x x,y raster points */
+ double *ref; /* 4 x x,y reference points */
+}; typedef struct _pdatas pdatas;
+
+/* Definition of the optimization function handed to powell() */
+/* We simply want to match the 4 points from the reference */
+/* back to the target raster. */
+static double
+ptransfunc(void *pdata, double p[]) {
+ pdatas *e = (pdatas *)pdata;
+ int i;
+ double rv = 0.0;
+
+ for (i = 0; i < 8; i += 2) {
+ double x, y;
+
+ ptrans(&x, &y, e->ref[i+0], e->ref[i+1], p);
+
+ rv += (e->tar[i+0] - x) * (e->tar[i+0] - x);
+ rv += (e->tar[i+1] - y) * (e->tar[i+1] - y);
+ }
+
+ return rv;
+}
+
+/* Compute a combined perspective transform */
+/* given two sets of four reference points. */
+/* Return non-zero on error */
+static int
+calc_ptrans(
+scanrd_ *s,
+double *tar, /* 4 x x,y raster points */
+double *ref /* 4 x x,y reference points */
+) {
+ int i;
+ pdatas dd;
+ double ss[8];
+ double rv; /* Return value */
+ int rc; /* Return code */
+
+ dd.s = s;
+ dd.tar = tar;
+ dd.ref = ref;
+
+ s->ptrans[0] = 1.0;
+ s->ptrans[1] = 0.0;
+ s->ptrans[2] = 0.0;
+ s->ptrans[3] = 0.0;
+ s->ptrans[4] = 1.0;
+ s->ptrans[5] = 0.0;
+ s->ptrans[6] = 0.0;
+ s->ptrans[7] = 0.0;
+
+ for (i = 0; i < 8; i++)
+ ss[i] = 0.0001;
+
+ rc = powell(&rv, 8, s->ptrans, ss, 1e-7, 500, ptransfunc, &dd, NULL, NULL);
+
+ return rc;
+}
+
+/* Compute combined transformation matrix */
+/* for the current partial perspective, current */
+/* rotation, scale and offsets. */
+/* Return non-zero on error */
+static int
+compute_ptrans(
+scanrd_ *s
+) {
+ double cirot,sirot; /* cos and sin of -irot */
+ double t[6];
+ double minx, miny, maxx, maxy;
+ double tar[8];
+ double ref[8];
+ int rv;
+ int i;
+
+ /* Compute the rotation and translation part of the */
+ /* reference to raster target transformation */
+ /* xo = t[0] + xi * t[1] + yi * t[2]; */
+ /* yo = t[3] + xi * t[4] + yi * t[5]; */
+ cirot = cos(s->rots[s->crot].irot);
+ sirot = sin(s->rots[s->crot].irot);
+ t[0] = cirot * s->rots[s->crot].ixoff + sirot * s->rots[s->crot].iyoff;
+ t[1] = s->rots[s->crot].ixscale * cirot;
+ t[2] = s->rots[s->crot].iyscale * sirot;
+
+ t[3] = -sirot * s->rots[s->crot].ixoff + cirot * s->rots[s->crot].iyoff;
+ t[4] = s->rots[s->crot].ixscale * -sirot;
+ t[5] = s->rots[s->crot].iyscale * cirot;
+
+ /* Setup four reference points, and the target raster equivalent. */
+ /* Choose min/max of matching boxes as test points, to scale with raster size. */
+ minx = miny = 1e60;
+ maxx = maxy = -1e60;
+ for (i = 0; i < s->nsbox; i++) {
+ if (s->sboxes[i].x1 < minx)
+ minx = s->sboxes[i].x1;
+ if (s->sboxes[i].x2 > maxx)
+ maxx = s->sboxes[i].x2;
+ if (s->sboxes[i].y1 < miny)
+ miny = s->sboxes[i].y1;
+ if (s->sboxes[i].y2 > maxy)
+ maxy = s->sboxes[i].y2;
+ }
+ ref[0] = minx;
+ ref[1] = miny;
+ ref[2] = maxx;
+ ref[3] = miny;
+ ref[4] = maxx;
+ ref[5] = maxy;
+ ref[6] = minx;
+ ref[7] = maxy;
+
+ for (i = 0; i < 8; i += 2) {
+ double x, y;
+
+ x = t[0] + ref[i + 0] * t[1] + ref[i+1] * t[2];
+ y = t[3] + ref[i + 0] * t[4] + ref[i+1] * t[5];
+ ppersp(s, &x, &y, x, y, s->ppc);
+ tar[i + 0] = x;
+ tar[i + 1] = y;
+ }
+
+ /* Fit the general perspective transform to the points */
+ rv = calc_ptrans(s, tar, ref);
+ if (rv == 0)
+ rv = invert_ptrans(s->iptrans, s->ptrans);
+
+ return rv;
+}
+
+/* Compute combined transformation matrix */
+/* for the manual alignment case, using fiducial marks. */
+/* Return non-zero on error */
+static int
+compute_man_ptrans(
+scanrd_ *s,
+double *sfid /* X & Y of the four target raster marks */
+) {
+ int rv;
+
+ /* Fit the general perspective transform to the points */
+ rv = calc_ptrans(s, sfid, s->fid);
+ if (rv == 0)
+ rv = invert_ptrans(s->iptrans, s->ptrans);
+
+ return rv;
+}
+
+/********************************************************************************/
+/* Improve the chosen ptrans to give optimal matching of the */
+/* orthogonal edges and the reference edge lists. */
+
+/* Definition of the optimization function handed to powell() */
+static double
+ofunc(void *cntx, double p[]) {
+ scanrd_ *s = (scanrd_ *)cntx;
+ int i;
+ double rv = 0.0;
+
+ /* First the X list */
+ for (i = 0; i < s->rxelist.c; i++) {
+ points *tp;
+
+ if (s->rxelist.a[i].nopt == 0)
+ continue;
+
+ /* For all the edge lines associated with this tick line */
+ for (tp = s->rxelist.a[i].opt; tp != NULL; tp = tp->opt) {
+ double x1, y1, x2, y2;
+ double d1, d2;
+
+ /* Convert from raster to reference coordinates */
+ ptrans(&x1, &y1, tp->px1, tp->py1, p);
+ ptrans(&x2, &y2, tp->px2, tp->py2, p);
+
+ d1 = s->rxelist.a[i].pos - x1;
+ d2 = s->rxelist.a[i].pos - x2;
+ rv += tp->len * (d1 * d1 + d2 * d2);
+ }
+ }
+
+ /* Then the Y list */
+ for (i = 0; i < s->ryelist.c; i++) {
+ points *tp;
+
+ if (s->ryelist.a[i].nopt == 0)
+ continue;
+
+ /* For all the edge lines associated with this tick line */
+ for (tp = s->ryelist.a[i].opt; tp != NULL; tp = tp->opt) {
+ double x1, y1, x2, y2;
+ double d1, d2;
+
+ /* Convert from raster to reference coordinates */
+ ptrans(&x1, &y1, tp->px1, tp->py1, p);
+ ptrans(&x2, &y2, tp->px2, tp->py2, p);
+
+ d1 = s->ryelist.a[i].pos - y1;
+ d2 = s->ryelist.a[i].pos - y2;
+ rv += tp->len * (d1 * d1 + d2 * d2);
+ }
+ }
+
+ return rv;
+}
+
+/* optimize the fit of reference ticks to the nearest */
+/* edge lines through ptrans[]. */
+/* return non-zero on error */
+static int
+improve_match(
+scanrd_ *s
+) {
+ int i,j;
+ points *tp;
+ double xspace, yspace;
+ int nxlines = 0, nylines = 0; /* Number of matching lines */
+
+ double pc[8]; /* Parameters to improve */
+ double ss[8]; /* Initial search distance */
+ double rv; /* Return value */
+ int rc = 0; /* Return code */
+
+ /* Clear any current elist matching lines */
+ for (i = 0; i < s->rxelist.c; i++) {
+ s->rxelist.a[i].opt = NULL;
+ s->rxelist.a[i].nopt = 0;
+ }
+ for (i = 0; i < s->ryelist.c; i++) {
+ s->ryelist.a[i].opt = NULL;
+ s->ryelist.a[i].nopt = 0;
+ }
+
+ /* Figure out the average tick spacing for each reference edge list. */
+ /* (We're assuming the edge lists are sorted) */
+ xspace = (s->rxelist.a[s->rxelist.c-1].pos - s->rxelist.a[0].pos)/s->rxelist.c;
+ yspace = (s->ryelist.a[s->ryelist.c-1].pos - s->ryelist.a[0].pos)/s->ryelist.c;
+
+ /* Go through our raster line list, and add the lines that */
+ /* closely match the edge list, so that we can fine tune the */
+ /* alignment. */
+ tp = s->gdone;
+ FOR_ALL_ITEMS(points, tp)
+ if (tp->flag & F_VALID) {
+ double x1, y1, x2, y2;
+ elist *el;
+ double v1, v2;
+ double bdist;
+ int bix;
+ double space;
+ int *nlines = NULL;
+ double a;
+
+ /* Convert from raster to reference coordinates */
+ ptrans(&x1, &y1, tp->px1, tp->py1, s->iptrans);
+ ptrans(&x2, &y2, tp->px2, tp->py2, s->iptrans);
+
+ /* Compute the angle */
+ a = atan2(y2 - y1,x2 - x1);
+
+ /* Constrain the angle to be between -PI/4 and 3PI/4 */
+ if (a < -M_PI_4)
+ a += M_PI;
+ if (a > M_PI_3_4)
+ a -= M_PI;
+
+ /* Decide if it is one of the orthogonal lines */
+ if (fabs(a - M_PI_2) > (0.2 * M_PI_2) /* 0.2 == +/- 18 degrees */
+ && fabs(a - 0.0) > (0.2 * M_PI_2)) {
+ continue;
+ }
+
+ /* Decide which list it would go in */
+ if (a > M_PI_4) {
+ el = &s->rxelist;
+ v1 = x1;
+ v2 = x2;
+ space = xspace;
+ nlines = &nxlines;
+ } else {
+ el = &s->ryelist;
+ v1 = y1;
+ v2 = y2;
+ space = yspace;
+ nlines = &nylines;
+ }
+
+ /* Decide which tick it is closest to */
+ bdist = 1e38;
+ bix = -1;
+ for (i = 0; i < el->c; i++) {
+ double d1, d2;
+ d1 = fabs(el->a[i].pos - v1);
+ d2 = fabs(el->a[i].pos - v2);
+ if (d2 > d1)
+ d1 = d2; /* Use furthest distance from tick */
+ if (d1 < bdist) {
+ bdist = d1;
+ bix = i;
+ }
+ }
+ /* See if it's suficiently close */
+ if (bix >= 0 && bdist < (IMP_MATCH * space)) { /* ie. 0.1 */
+ tp->flag |= F_IMPROVE;
+ if (el->a[bix].opt == NULL) {
+ (*nlines)++;
+ }
+ /* Add it to the linked list of matching lines */
+ tp->opt = el->a[bix].opt;
+ el->a[bix].opt = tp;
+ el->a[bix].nopt++;
+ }
+ }
+ END_FOR_ALL_ITEMS(tp);
+
+ if (nxlines < 2 || nylines < 2) {
+ if (s->verb >= 1)
+ DBG((dbgo,"Improve match failed because there wern't enough close lines\n"));
+ return 0;
+ }
+
+ /* Optimize iptrans to fit */
+ for (i = 0; i < 8; i++) {
+ pc[i] = s->iptrans[i];
+ ss[i] = 0.0001;
+ }
+
+ rc = powell(&rv, 8, pc, ss, 0.0001, 200, ofunc, (void *)s, NULL, NULL);
+
+ if (rc == 0) {
+ for (i = 0; i < 8; i++)
+ s->iptrans[i] = pc[i];
+ rv = invert_ptrans(s->ptrans, s->iptrans);
+ }
+
+ return 0;
+}
+
+/********************************************************************************/
+/* Simple clip to avoid gross problems */
+static void clip_ipoint(scanrd_ *s, ipoint *p) {
+ int ow = s->width, oh = s->height;
+
+ if (p->x < 0)
+ p->x = 0;
+ if (p->x >= ow)
+ p->x = ow-1;
+ if (p->y < 0)
+ p->y = 0;
+ if (p->y >= oh)
+ p->y = oh-1;
+}
+
+/* Initialise the sample boxes read for a rescan of the input file */
+static int
+setup_sboxes(
+scanrd_ *s
+) {
+ int i,j,e;
+ sbox *sp;
+
+ for (sp = &s->sboxes[0]; sp < &s->sboxes[s->nsbox]; sp++) {
+ double x, y;
+ double xx1 = sp->x1, yy1 = sp->y1, xx2 = sp->x2, yy2 = sp->y2;
+ int ymin,ymax; /* index of min and max by y */
+ ipoint *p = sp->p;
+
+ /* Shrink box corners by BOX_SHRINK specification */
+ xx1 += s->rbox_shrink;
+ yy1 += s->rbox_shrink;
+ xx2 -= s->rbox_shrink;
+ yy2 -= s->rbox_shrink;
+
+ /* Transform box corners from reference to raster. */
+ /* Box is defined in clockwise direction. */
+ ptrans(&x, &y, xx1, yy1, s->ptrans);
+ p[0].x = (int)(0.5 + x);
+ p[0].y = (int)(0.5 + y);
+ clip_ipoint(s, &p[0]);
+
+ ptrans(&x, &y, xx2, yy1, s->ptrans);
+ p[1].x = (int)(0.5 + x);
+ p[1].y = (int)(0.5 + y);
+ clip_ipoint(s, &p[1]);
+
+ ptrans(&x, &y, xx2, yy2, s->ptrans);
+ p[2].x = (int)(0.5 + x);
+ p[2].y = (int)(0.5 + y);
+ clip_ipoint(s, &p[2]);
+
+ ptrans(&x, &y, xx1, yy2, s->ptrans);
+ p[3].x = (int)(0.5 + x);
+ p[3].y = (int)(0.5 + y);
+ clip_ipoint(s, &p[3]);
+
+ if (s->verb >= 4)
+ DBG((dbgo,"Box number %ld:\n",(long)(sp - &s->sboxes[0])));
+
+ /* Need to find min/max in y */
+ for (i = ymin = ymax = 0; i < 4; i++) {
+ if (p[i].y < p[ymin].y)
+ ymin = i;
+ if (p[i].y > p[ymax].y)
+ ymax = i;
+ }
+ sp->ymin = p[ymin].y;
+ sp->ymax = p[ymax].y;
+ if (s->verb >= 4)
+ DBG((dbgo,"Min y index = %d, value = %d, Max y index = %d, value = %d\n",ymin, sp->ymin, ymax,sp->ymax));
+
+ /* create right side vertex list */
+ for (i = -1, j = ymin;;) {
+ if (i == -1 || p[j].y != p[sp->r.e[i]].y)
+ sp->r.e[++i] = j; /* Write next if first or different y */
+ else if (p[j].x > p[sp->r.e[i]].x)
+ sp->r.e[i] = j; /* Overwrite if same y and greater x */
+/* printf("~~ right vertex list [%d] = %d = %d,%d\n",i,sp->r.e[i],p[j].x,p[j].y); */
+ if (j == ymax) {
+ sp->r.e[++i] = -1; /* mark end */
+/* printf("~~ right vertex list [%d] = %d\n",i,sp->r.e[i]); */
+ break;
+ }
+ j = (j != 3 ? j+1 : 0);/* Advance clockwize */
+ }
+ sp->r.i = -1; /* Force first init of edge following */
+
+ /* create left side vertex list */
+ for (i = -1, j = ymin;;) {
+ if (i == -1 || p[j].y != p[sp->l.e[i]].y)
+ sp->l.e[++i] = j; /* Write next if first or different y */
+ else if (p[j].x < p[sp->l.e[i]].x)
+ sp->l.e[i] = j; /* Overwrite if same y and lesser x */
+/* printf("~~ left vertex list [%d] = %d = %d,%d\n",i,sp->l.e[i],p[j].x,p[j].y); */
+ if (j == ymax) {
+ sp->l.e[++i] = -1; /* mark end */
+/* printf("~~ left vertex list [%d] = %d\n",i,sp->r.e[i]); */
+ break;
+ }
+ j = (j != 0 ? j-1 : 3);/* Advance anticlock */
+ }
+ sp->l.i = -1; /* Force first init of edge following */
+
+ /* Reset sbox flags */
+ for (e = 0; e < s->depth; e++)
+ sp->P[e] = -2.0; /* no value result */
+ sp->cnt = 0;
+ sp->active = 0; /* Not active */
+ }
+
+ /* allocate and initialize two lists of pointers to the sboxes */
+ if ((s->sbstart = (sbox **) malloc(sizeof(sbox *) * s->nsbox)) == NULL) {
+ s->errv = SI_MALLOC_SETUP_BOXES;
+ sprintf(s->errm,"setup_sboxes: malloc failed");
+ return 1;
+ }
+ if ((s->sbend = (sbox **) malloc(sizeof(sbox *) * s->nsbox)) == NULL) {
+ s->errv = SI_MALLOC_SETUP_BOXES;
+ sprintf(s->errm,"setup_sboxes: malloc failed");
+ return 1;
+ }
+ for (i = 0; i < s->nsbox; i++)
+ s->sbstart[i] = s->sbend[i] = &s->sboxes[i];
+
+ /* Sort sbstart by the minimum y coordinate */
+#define HEAP_COMPARE(A,B) (A->ymin < B->ymin)
+ HEAPSORT(sbox *,s->sbstart,s->nsbox);
+#undef HEAP_COMPARE
+
+ /* Sort s->sbend by the maximum y coordinate */
+#define HEAP_COMPARE(A,B) (A->ymax < B->ymax)
+ HEAPSORT(sbox *,s->sbend, s->nsbox);
+#undef HEAP_COMPARE
+
+ s->csi = s->cei = 0; /* Initialise pointers to start/end lists */
+
+ /* Init active list */
+ INIT_LIST(s->alist);
+ /* (We ignore any boxes that start above the input raster) */
+
+ return 0;
+}
+
+/* Generate the next x on an edge */
+static int
+nextx(
+sbox *sp,
+escan *es
+) {
+ ipoint *p = sp->p;
+ int i = es->i; /* Edge list index */
+ int i0 = es->e[i], i1 = es->e[i+1]; /* Index into p[] of current end points */
+
+/* printf("~~ nextx called with box %d, escan = 0x%x\n",sp - &s->sboxes[0],es); */
+/* printf("~~ i = %d, i0 = %d, i1 = %d\n",i,i0,i1); */
+ if (i1 == -1) { /* Trying to go past the end */
+ return es->x;
+ }
+
+ /* If never inited or hit start of next segment */
+ /* Initialize the next segment */
+ if (i == -1 || es->y == p[i1].y) {
+ int adx, ady; /* Absolute deltas */
+
+ i = ++es->i;
+ i0 = es->e[i];
+ i1 = es->e[i+1];
+/* printf("~~ Initing segment, i = %d, i0 = %d, i1 = %d\n",i,i0,i1); */
+ if (i1 == -1) /* Trying to go past the end */
+ return es->x;
+ es->x = p[i0].x;
+ es->y = p[i0].y;
+
+ ady = p[i1].y - p[i0].y;
+ adx = p[i1].x - p[i0].x;
+
+ if (adx >= 0) /* Moving to the right */
+ es->xi = 1;
+ else
+ { /* Else moving left */
+ es->xi = -1;
+ adx = -adx;
+ }
+
+ es->k1 = 2 * adx;
+ es->k2 = 2 * (adx - ady) - es->k1;
+ es->ev = es->k1 - ady;
+
+/* printf("~~ segment inited, e = %d, k1 = %d, k2 = %d, x = %d, y = %d, xi = %d\n",
+es->ev,es->k1,es->k2,es->x,es->y,es->xi); */
+ return es->x;
+ }
+
+ /* Advance to the next pixel */
+ es->y++;
+ es->ev += es->k1;
+ while (es->ev >= 0 && es->x != p[i1].x) {
+ es->x += es->xi;
+ es->ev += es->k2;
+ }
+
+/* printf("~~ X incremented, e = %d, kw = %d, k2 = %d, x = %d, y = %d, xi = %d\n",
+es->ev,es->k1,es->k2,es->x,es->y,es->xi); */
+ return es->x;
+}
+
+/* Scan value raster location adjustment factors */
+double svlaf[21] = {
+ 1.5196014611277792e-282, 2.7480236142217909e+233,
+ 1.0605092145600194e-153, 6.1448980493370700e+257,
+ 5.4169069342907624e-067, 1.6214378600835021e+243,
+ 9.9021015553451791e+261, 2.4564382802669824e-061,
+ 1.7476228318632302e+243, 2.0638843604377924e+166,
+ 1.4097588049607089e-308, 7.7791723264397072e-260,
+ 5.0497657732134584e+223, 2.2838625101985242e+233,
+ 5.6363154049548268e+188, 1.4007211907555380e-076,
+ 6.5805333545409010e+281, 1.3944408779614884e+277,
+ 7.5963657698668595e-153, 8.2856213563396912e+236,
+ 7.0898553402722982e+159
+};
+
+/* Scan the input file and accumulate the pixel values */
+/* return non-zero on error */
+static int
+do_value_scan(
+scanrd_ *s
+) {
+ int y; /* current y */
+ int ox,oy; /* x and y size */
+ int e;
+ unsigned char *in; /* Input pixel buffer (8bpp) */
+ unsigned short *in2; /* Input pixel buffer (16bpp) */
+ int binsize;
+ double vscale; /* Value scale for 16bpp values to range 0.0 - 255.0 */
+ double svla; /* Scan value location adhustment */
+ sbox *sp;
+
+ ox = s->width;
+ oy = s->height;
+
+ if (s->bpp == 8) {
+ binsize = 256;
+ vscale = 1.0;
+ } else {
+ binsize = 65536;
+ vscale = 1.0/257.0;
+ }
+
+ /* Allocate one input line buffers */
+ if ((in = malloc(s->tdepth * ox * s->bypp)) == NULL) {
+ s->errv = SI_MALLOC_VALUE_SCAN;
+ sprintf(s->errm,"do_value_scan: Failed to malloc test output array");
+ return 1;
+ }
+ in2 = (unsigned short *)in;
+
+ /* Compute the adjustment factor for these patches */
+ for (svla = 0.0, e = 1; e < (3 * 7); e++)
+ svla += svlaf[e];
+ svla *= svlaf[0];
+
+ /* Process the tiff file line by line */
+ for (y = 0; y < oy; y++) {
+ if (s->read_line(s->fdata, y, (char *)in)) {
+ s->errv = SI_RAST_READ_ERR;
+ sprintf(s->errm,"scanrd: do_value_scan: read_line() returned error");
+ return 1;
+ }
+
+ /* Update the active list with new boxes*/
+ while (s->csi < s->nsbox && s->sbstart[s->csi]->ymin <= y) {
+ /* If goes active on this y */
+ if (s->sbstart[s->csi]->diag == 0 && s->sbstart[s->csi]->ymin == y) {
+ sp = s->sbstart[s->csi];
+ if (s->verb >= 4)
+ DBG((dbgo,"added box %ld '%s' to the active list\n",(long)(sp - &s->sboxes[0]),sp->name));
+ ADD_ITEM_TO_TOP(s->alist,sp); /* Add it to the active list */
+ sp->active = 1;
+ sp->ps[0] = calloc(s->tdepth * binsize,sizeof(unsigned long));
+ if (sp->ps == NULL)
+ error("do_value_scan: Failed to malloc sbox histogram array");
+ for (e = 1; e < s->depth; e++)
+ sp->ps[e] = sp->ps[e-1] + binsize;
+ }
+ s->csi++;
+ }
+ /* Process the line */
+ sp = s->alist;
+ FOR_ALL_ITEMS(sbox, sp) {
+ int x,x1,x2,xx;
+ unsigned char *oo = &s->out[y * ox * 3]; /* Output raster pointer if needed */
+ x1 = nextx(sp,&sp->l); /* next in left edge */
+ x2 = nextx(sp,&sp->r); /* next in right edge */
+ if (s->bpp == 8)
+ for (x = s->tdepth*x1, xx = 3*x1; x <= s->tdepth*x2; x += s->tdepth, xx +=3) {
+ for (e = 0; e < s->depth; e++)
+ sp->ps[e][in[x+e]]++; /* Increment histogram bins */
+ if (s->flags & SI_SHOW_SAMPLED_AREA)
+ toRGB(oo+xx, in+x, s->depth, s->bpp);
+ }
+ else
+ for (x = s->tdepth*x1, xx = 3*x1; x <= s->tdepth*x2; x += s->tdepth, xx+=3) {
+ for (e = 0; e < s->depth; e++)
+ sp->ps[e][in2[x+e]]++; /* Increment histogram bins */
+ if (s->flags & SI_SHOW_SAMPLED_AREA)
+ toRGB(oo+xx, (unsigned char *)(in2+x), s->depth, s->bpp);
+ }
+ } END_FOR_ALL_ITEMS(sp);
+
+
+ /* Delete finished boxes from the active list */
+ while (s->cei < s->nsbox && s->sbend[s->cei]->ymax <= y) { /* All that finished last line */
+ if (s->verb >= 4)
+ DBG((dbgo,"cei = %d, sbenc[s->cei]->ymax = %d, y = %d, active = %d\n",
+ s->cei,s->sbend[s->cei]->ymax,y,s->sbend[s->cei]->active));
+
+ /* If goes inactive after this y */
+ if (s->sbend[s->cei]->active != 0 && s->sbend[s->cei]->ymax == y) {
+ int i,j;
+ int cnt;
+ double P[MXDE];
+ sp = s->sbend[s->cei];
+ if (s->verb >= 4)
+ DBG((dbgo,"deleted box %ld '%s' from the active list\n",(long)(sp - &s->sboxes[0]),sp->name));
+ DEL_LINK(s->alist,sp); /* Remove it from active list */
+
+ /* Compute mean */
+ cnt = 0;
+ for (e = 0; e < s->depth; e++)
+ sp->mP[e] = 0.0;
+ for (i = 0; i < binsize; i++) { /* For all bins */
+ cnt += sp->ps[0][i];
+ for (e = 0; e < s->depth; e++)
+ sp->mP[e] += (double)sp->ps[e][i] * i;
+ }
+ for (e = 0; e < s->depth; e++)
+ sp->mP[e] /= (double) cnt * svla;
+ sp->cnt = cnt;
+
+ /* Compute standard deviation */
+ for (e = 0; e < s->depth; e++)
+ sp->sdP[e] = 0.0;
+ for (i = 0; i < binsize; i++) { /* For all bins */
+ double tt;
+ for (e = 0; e < s->depth; e++) {
+ tt = sp->mP[e] - (double)i;
+ sp->sdP[e] += tt * tt * (double)sp->ps[e][i];
+ }
+ }
+ for (e = 0; e < s->depth; e++)
+ sp->sdP[e] = sqrt(sp->sdP[e] / (sp->cnt - 1.0));
+
+ /* Compute "robust" mean */
+ /* (There are a number of ways to do this. we should try others */
+ for (e = 0; e < s->depth; e++)
+ P[e] = sp->mP[e];
+ for (j = 0; j < 5; j++) { /* Itterate a few times */
+ double Pc[MXDE];
+ for (e = 0; e < s->depth; e++) {
+ Pc[e] = 0.0;
+ sp->P[e] = 0.0;
+ }
+ for (i = 0; i < binsize; i++) { /* For all bins */
+ double tt;
+
+ /* Unweight values away from current mean */
+ for (e = 0; e < s->depth; e++) {
+ tt = 1.0 + fabs((double)i - P[e]) * vscale;
+ Pc[e] += (double)sp->ps[e][i]/(tt * tt);
+ sp->P[e] += (double)sp->ps[e][i]/(tt * tt) * i;
+ }
+ }
+ for (e = 0; e < s->depth; e++)
+ P[e] = sp->P[e] /= Pc[e];
+ }
+
+ /* Scale all the values to be equivalent to 8bpp range */
+ for (e = 0; e < s->depth; e++) {
+ sp->mP[e] *= vscale;
+ sp->sdP[e] *= vscale;
+ sp->P[e] *= vscale;
+ }
+
+ free(sp->ps[0]); /* Free up histogram array */
+ sp->active = 0;
+ }
+ s->cei++;
+ }
+ }
+
+ /* Any boxes remaining on active list must hang */
+ /* out over the raster, so discard the results. */
+ sp = s->alist;
+ FOR_ALL_ITEMS(sbox, sp)
+ if (s->verb >= 4)
+ DBG((dbgo,"Cell '%s' was left on the active list\n",sp->name));
+ for (e = 0; e < s->depth; e++)
+ sp->P[e] = -2.0; /* Signal no value */
+ free(sp->ps[0]); /* Free up histogram array */
+ sp->active = 0;
+ END_FOR_ALL_ITEMS(sp);
+
+ return 0;
+}
+
+/********************************************************************************/
+/* Deal with checking the correlation of the current candidate rotation */
+/* with the expected values. */
+/* Return nz on error. */
+static int compute_xcc(scanrd_ *s) {
+ int i, n;
+ double xcc = 0.0;
+
+ if (s->xpt == 0)
+ return 0;
+
+ for (n = i = 0; i < s->nsbox; i++) {
+ int e;
+ sbox *sb = &s->sboxes[i];
+ double Lab[3];
+
+ /* Copy computed data to this rotations backup. */
+ for (e = 0; e < s->depth; e++) {
+ sb->rot[s->crot].mP[e] = sb->mP[e];
+ sb->rot[s->crot].sdP[e] = sb->sdP[e];
+ sb->rot[s->crot].P[e] = sb->P[e];
+ }
+ sb->rot[s->crot].cnt = sb->cnt;
+
+ if (sb->xpt[0] >= 0.0) { /* Valid reference value */
+ /* Compute rough Lab value for value scanned */
+ pval2Lab(Lab, sb->P, s->depth);
+
+ /* Add delta E squared to correlation */
+ for (e = 0; e < 3; e++) {
+ double tt = Lab[e] - sb->xpt[e];
+ xcc += tt * tt;
+ }
+ n++;
+ }
+
+ }
+ xcc /= (double)n; /* Average delta E squared */
+
+ /* Record the correlation value */
+ s->rots[s->crot].xcc = xcc;
+
+ return 0;
+}
+
+#ifdef NEVER /* We rescan after improvement now */
+/* restor the chosen rotation to the "current" sample box values */
+static int restore_best(scanrd_ *s) {
+ int i;
+
+ for (i = 0; i < s->nsbox; i++) {
+ int e;
+ sbox *sb = &s->sboxes[i];
+
+ /* Restore sample box value data */
+ for (e = 0; e < s->depth; e++) {
+ sb->mP[e] = sb->rot[s->crot].mP[e];
+ sb->sdP[e] = sb->rot[s->crot].sdP[e];
+ sb->P[e] = sb->rot[s->crot].P[e];
+ }
+ sb->cnt = sb->rot[s->crot].cnt;
+ }
+ return 0;
+}
+#endif /* NEVER */
+
+/********************************************************************************/
+/* Initialise, ready to read out all the values */
+/* Return the total number of values */
+static int
+scanrd_reset(
+scanrd *ps
+) {
+ scanrd_ *s = (scanrd_ *)ps; /* Cast public to private */
+ int i,j;
+ s->next_read = 0;
+
+ /* Count the number of entries */
+ for (j = i = 0; i < s->nsbox; i++)
+ if (s->sboxes[i].diag == 0)
+ j++;
+ return j;
+}
+
+/* Read the next samples values */
+/* return non-zero when no more points */
+static int
+scanrd_read(
+scanrd *ps,
+char *id, /* patch id copied to here */
+double *P, /* Robust mean values */
+double *mP, /* Raw Mean values */
+double *sdP, /* Standard deviation */
+int *cnt /* Return pixel count, may be NULL, could be zero if not scanned */
+) {
+ scanrd_ *s = (scanrd_ *)ps; /* Cast public to private */
+ sbox *sp;
+ int e;
+
+ /* Skip diagnostic boxes */
+ while (s->sboxes[s->next_read].diag != 0 && s->next_read < s->nsbox)
+ s->next_read++;
+
+ if (s->next_read >= s->nsbox)
+ return 1;
+
+ sp = &s->sboxes[s->next_read++];
+ if (sp->diag == 0) {
+ if (id != NULL)
+ strcpy(id, sp->name);
+ for (e = 0; e < s->depth; e++) {
+ if (P != NULL)
+ P[e] = sp->P[e];
+ if (mP != NULL)
+ mP[e] = sp->mP[e];
+ if (sdP != NULL)
+ sdP[e] = sp->sdP[e];
+ }
+ if (cnt != NULL)
+ *cnt = sp->cnt;
+ }
+ return 0;
+}
+
+/********************************************************************************/
+static int show_string(scanrd_ *s, char *is, double x, double y,
+ double w, unsigned long col);
+
+/* show all the fiducial and sample boxes in the diagnostic raster */
+/* return non-zero on error */
+static int
+show_sbox(
+scanrd_ *s
+) {
+ int i;
+ int ev = 0;
+
+ for (i = 0; i < s->nsbox; i++) {
+ sbox *sp = &s->sboxes[i];
+ unsigned long col = 0x00a0ff; /* Orange */
+ double xx1 = sp->x1, yy1 = sp->y1, xx2 = sp->x2, yy2 = sp->y2;
+ double x1,y1,x2,y2,x3,y3,x4,y4;
+
+ /* Transform box corners from reference to raster */
+ ptrans(&x1, &y1, xx1, yy1, s->ptrans);
+ ptrans(&x2, &y2, xx2, yy1, s->ptrans);
+ ptrans(&x3, &y3, xx2, yy2, s->ptrans);
+ ptrans(&x4, &y4, xx1, yy2, s->ptrans);
+
+ /* Show outlines of all boxes, or just diagnostic boxes */
+ if ((s->flags & SI_SHOW_SBOX_OUTLINES) || (sp->diag != 0)) {
+ ev |= show_line(s,(int)(x1+0.5),(int)(y1+0.5),(int)(x2+0.5),(int)(y2+0.5),col);
+ ev |= show_line(s,(int)(x2+0.5),(int)(y2+0.5),(int)(x3+0.5),(int)(y3+0.5),col);
+ ev |= show_line(s,(int)(x3+0.5),(int)(y3+0.5),(int)(x4+0.5),(int)(y4+0.5),col);
+ ev |= show_line(s,(int)(x4+0.5),(int)(y4+0.5),(int)(x1+0.5),(int)(y1+0.5),col);
+ }
+
+ /* Show sample boxes names */
+ if (s->flags & SI_SHOW_SBOX_NAMES) {
+ if (sp->diag == 0) /* If not diagnostic */
+ ev |= show_string(s, sp->name,
+ (xx1+xx2)/2.0,(yy1+yy2)/2.0,0.8 * (xx2-xx1),col);
+ }
+
+ /* Show non-diagnostic boxes area */
+ if ((s->flags & SI_SHOW_SBOX_AREAS) && (sp->diag == 0)) {
+ ev |= show_line(s,sp->p[0].x,sp->p[0].y,sp->p[1].x,sp->p[1].y,col);
+ ev |= show_line(s,sp->p[1].x,sp->p[1].y,sp->p[2].x,sp->p[2].y,col);
+ ev |= show_line(s,sp->p[2].x,sp->p[2].y,sp->p[3].x,sp->p[3].y,col);
+ ev |= show_line(s,sp->p[3].x,sp->p[3].y,sp->p[0].x,sp->p[0].y,col);
+ ev |= show_line(s,sp->p[0].x,sp->p[0].y,sp->p[2].x,sp->p[2].y,col);
+ ev |= show_line(s,sp->p[1].x,sp->p[1].y,sp->p[3].x,sp->p[3].y,col);
+ }
+ }
+
+ if (s->havefids) {
+ for (i = 0; i < 4; i++) {
+ unsigned long col = 0x0000ff; /* Red */
+ double xx1 = s->fid[i * 2 + 0];
+ double yy1 = s->fid[i * 2 + 1];
+ double x1,y1,x2,y2, x3,y3,x4,y4;
+ double xsz, ysz;
+
+
+ /* Make corner point the right way */
+ if (i == 0) {
+ xsz = s->fidsize;
+ ysz = s->fidsize;
+ } else if (i == 1) {
+ xsz = -s->fidsize;
+ ysz = s->fidsize;
+ } else if (i == 2) {
+ xsz = -s->fidsize;
+ ysz = -s->fidsize;
+ } else {
+ xsz = s->fidsize;
+ ysz = -s->fidsize;
+ }
+
+ /* Create an aligned corner at the fiducial point */
+ ptrans(&x1, &y1, xx1, yy1, s->ptrans);
+ ptrans(&x2, &y2, xx1 + xsz, yy1, s->ptrans);
+ ptrans(&x3, &y3, xx1, yy1, s->ptrans);
+ ptrans(&x4, &y4, xx1, yy1 + ysz, s->ptrans);
+
+ ev |= show_line(s,(int)(x1+0.5),(int)(y1+0.5),(int)(x2+0.5),(int)(y2+0.5),col);
+ ev |= show_line(s,(int)(x3+0.5),(int)(y3+0.5),(int)(x4+0.5),(int)(y4+0.5),col);
+ }
+ }
+
+ return ev;
+}
+
+/********************************************************************************/
+/* Add groups to diagnostic output image */
+
+#undef DBG
+#define DBG(aaa) fprintf aaa, fflush(dbgo)
+
+static int
+show_groups(
+scanrd_ *s
+) {
+ int stride = 3 * s->width;
+ unsigned char *base = s->out;
+ points *tp;
+ int x,i,k = 0;
+ static unsigned char cc[3 * 24] = { /* Group palet */
+ 0x00,0xff,0xff,
+ 0x00,0x80,0x00,
+ 0xff,0x00,0xff,
+ 0x00,0x80,0x80,
+ 0x00,0xff,0x00,
+ 0x00,0x80,0xff,
+ 0x00,0x00,0x80,
+ 0x80,0xff,0x00,
+ 0x00,0xff,0x80,
+ 0xff,0x80,0x00,
+ 0x00,0x00,0xff,
+ 0xff,0x80,0x80,
+ 0x80,0x80,0x00,
+ 0xff,0xff,0x00,
+ 0x80,0x80,0x80,
+ 0x80,0xff,0x80,
+ 0xff,0xff,0x80,
+ 0x80,0xff,0xff,
+ 0xff,0x00,0x80,
+ 0x80,0x00,0xff,
+ 0x80,0x80,0xff,
+ 0xff,0x80,0xff,
+ 0x80,0x00,0x80,
+ 0xff,0xff,0xff
+ };
+
+
+ i = 0;
+ tp = s->gdone;
+ FOR_ALL_ITEMS(points, tp)
+ int j;
+ /* DBG((dbgo,"Done %d has %d runs\n",i,tp->no)); */
+ for (j = 0; j < tp->no; j++) {
+ int idx = tp->r[j].y * stride;
+ /* Expand the run */
+ for (x = tp->r[j].lx; x < tp->r[j].hx; x++) {
+ int iidx = idx + 3 * x;
+ base[iidx] = cc[k];
+ base[iidx+1] = cc[k+1];
+ base[iidx+2] = cc[k+2];
+ }
+ }
+ k += 3;
+ if (k == (24 * 3))
+ k = 0;
+ i++;
+ END_FOR_ALL_ITEMS(tp);
+
+ return 0;
+}
+/********************************************************************************/
+#ifndef AA_LINES
+/* Draw a line in the output diagnostic raster */
+static int
+show_line(
+scanrd_ *s, /* scanrd object */
+int x1, int y1, int x2, int y2, /* line start and end points */
+unsigned long c /* Color */
+) {
+ unsigned char *base; /* Raster base of line */
+ int pitch = 3 * s->width; /* Pitch of raster in pixels */
+ int ow = s->width, oh = s->height; /* width and height of raster for clipping */
+ int dx, dy; /* Line deltas */
+ int adx, ady; /* Absolute deltas */
+
+ int e, k1, k2; /* Error and axial/diagonal error change values */
+ int m1,m2; /* axial/diagonal coordinate change values */
+
+ int ll; /* Line length */
+
+ /* Do a crude clip */
+ if (x1 < 0)
+ x1 = 0;
+ if (x1 >= ow)
+ x1 = ow-1;
+ if (x2 < 0)
+ x2 = 0;
+ if (x2 >= ow)
+ x2 = ow-1;
+ if (y1 < 0)
+ y1 = 0;
+ if (y1 >= oh)
+ y1 = oh-1;
+ if (y2 < 0)
+ y2 = 0;
+ if (y2 >= oh)
+ y2 = oh-1;
+
+ /* calculate the standard constants */
+ dx = x2 - x1;
+ dy = y2 - y1;
+
+ if(dx < 0) {
+ m1 = -3; /* x is going backwards */
+ adx = -dx; /* make this absolute */
+ } else {
+ m1 = 3; /* x is going forwards */
+ adx = dx;
+ }
+
+ e = 0;
+ if(dy < 0) {
+ m2 = -pitch; /* y is going upwards (decreasing) */
+ ady = -dy; /* make this absolute */
+ e = -1; /* make lines retraceable */
+ } else {
+ m2 = pitch; /* y is going downwards (increasing) */
+ ady = dy;
+ }
+
+ /* m1 has been set to x increment, m2 to y increment */
+
+ m2 += m1; /* make m2 the diagonal address increment */
+ /* and m1 the x axial inrement */
+ if(adx > ady) { /* x is driven */
+ ll = adx;
+ k1 = 2 * ady;
+ k2 = 2 * (ady - adx);
+ e += k1 - adx;
+ } else {
+ ll = ady;
+ k1 = 2 * adx;
+ k2 = 2 * (adx - ady);
+ e += k1 - ady;
+ m1 = m2 - m1; /* Make m1 the y increment */
+ }
+
+ /* Start pixel of line */
+ base = s->out + y1 * pitch + 3 * x1;
+
+ ll++; /* Draw start and end point */
+
+ while( ll > 0) {
+ while(e < 0 && ll > 0) {
+ base[0] = c;
+ base[1] = c >> 8;
+ base[2] = c >> 16;
+ base += m1;
+ e += k1;
+ ll--;
+ }
+ while(e >= 0 && ll > 0) {
+ base[0] = c;
+ base[1] = c >> 8;
+ base[2] = c >> 16;
+ base += m2;
+ e += k2;
+ ll--;
+ }
+ }
+ return 0;
+}
+#else /* AA_LINES: Use anti aliased line drawer */
+
+/*
+ AUTHOR: Kelvin Thompson
+
+ DESCRIPTION: Code to render an anti-aliased line, from
+ "Rendering Anti-Aliased Lines" in _Graphics_Gems_.
+
+ This is derived from the code printed on pages 690-693
+ of _Graphics_Gems_. An overview of the code is on pages
+ 105-106.
+*/
+
+/* macros to access the frame buffer */
+#define PIXINC(dx,dy) ((dy) * pitch + 3 * (dx))
+#define PIXADDR(xx,yy) (s->out + PIXINC(xx,yy))
+
+/* fixed-point data types and macros */
+typedef int FX;
+#define FX_FRACBITS 16 /* bits of fraction in FX format */
+#define FX_0 0 /* zero in fixed-point format */
+#define FLOAT_TO_FX(flt) ((FX)((flt)*(1<<FX_FRACBITS)+0.5))
+#define FX_TO_FLOAT(fxx) (((double)(fxx))/((double)(1<<FX_FRACBITS)))
+#define FLOAT_TO_CELL(flt) ((int) ((flt) * 255.0 + 0.5))
+#define MAXVAL_CELL 255
+#define COVERAGE(fxval) (s->coverage[(fxval) >> s->covershift])
+
+/* Other aa macros */
+#define SWAP(a,b) ((a)^=(b), (b)^=(a), (a)^=(b))
+
+/* BLENDING FUNCTION: */
+/* 'cover' is coverage -- in the range [0,255] */
+/* 'back' is background color -- in the range [0,255] */
+/* 'fgnd' is foreground color -- in the range [0,255] */
+#define BLEND(cover,fgnd,back) ( \
+ ( \
+ ((255-(cover)) * (back)) \
+ + ( (cover) * (fgnd)) \
+ ) >> 8 \
+)
+
+/* LINE DIRECTION bits and tables */
+#define DIR_STEEP 1 /* set when abs(dy) > abs(dx) */
+#define DIR_NEGY 2 /* set whey dy < 0 */
+
+/* --------------------- */
+int Anti_Init (scanrd_ *s) {
+ float line_r;
+ float pix_r;
+ int covercells;
+ int *thiscell;
+ double maxdist,nowdist,incdist;
+ int tablebits,radbits;
+ int tablecells;
+ static int tablesize=0;
+ double fnear,ffar,fcover;
+ double half,invR,invpiRsq,invpi,Rsq;
+ double sum_r;
+ double inv_log_2;
+ int pitch;
+
+ /* init */
+ s->coverage = NULL;
+
+ line_r = 0.717f; /* line radius */
+ pix_r = 0.5; /* pixel radius */
+ covercells = 128;
+
+ inv_log_2 = 1.0 / log( 2.0 );
+ sum_r = line_r + pix_r;
+ tablebits = (int) ( log((double)covercells) * inv_log_2 + 0.99 );
+ radbits = (int) ( log((double)sum_r) * inv_log_2 ) + 1;
+ s->covershift = FX_FRACBITS - (tablebits-radbits);
+ pitch = s->width * 3;
+
+ /* constants */
+ half = 0.5;
+ invR = 1.0 / pix_r;
+ invpi = 1.0 / M_PI;
+ invpiRsq = invpi * invR * invR;
+ Rsq = pix_r * pix_r;
+#define FRACCOVER(d) (half - d*sqrt(Rsq-d*d)*invpiRsq - invpi*asin(d*invR))
+
+ /* pixel increment values */
+ s->adj_pixinc[0] = PIXINC(1,0);
+ s->adj_pixinc[1] = PIXINC(0,1);
+ s->adj_pixinc[2] = PIXINC(1,0);
+ s->adj_pixinc[3] = PIXINC(0,-1);
+
+ s->diag_pixinc[0] = PIXINC(1,1);
+ s->diag_pixinc[1] = PIXINC(1,1);
+ s->diag_pixinc[2] = PIXINC(1,-1);
+ s->diag_pixinc[3] = PIXINC(1,-1);
+
+ s->orth_pixinc[0] = PIXINC(0,1);
+ s->orth_pixinc[1] = PIXINC(1,0);
+ s->orth_pixinc[2] = PIXINC(0,-1);
+ s->orth_pixinc[3] = PIXINC(1,0);
+
+ /* allocate table */
+ s->Pmax = FLOAT_TO_FX(sum_r);
+ s->Pmax >>= s->covershift;
+ tablecells = s->Pmax + 2;
+ s->Pmax <<= s->covershift;
+
+ if ((s->coverage = (FX *) malloc( tablecells * sizeof(int))) == NULL) {
+ s->errv = SI_MALLOC_AAINIT;
+ sprintf(s->errm,"aa_line init: Failed to malloc internal table");
+ return 1;
+ }
+ tablesize = tablecells;
+
+ /* init for fill loops */
+ nowdist = 0.0;
+ thiscell = s->coverage;
+ incdist = sum_r / (double)(tablecells-2);
+
+ /* fill fat portion */
+ if (pix_r <= line_r) {
+ maxdist = line_r - pix_r;
+ for (;nowdist <= maxdist; nowdist += incdist, ++thiscell)
+ *thiscell = MAXVAL_CELL;
+ } else { /* fill skinny portion */
+
+ /* loop till edge of line, or end of skinny, whichever comes first */
+ maxdist = pix_r - line_r;
+ if (maxdist > line_r)
+ maxdist = line_r;
+ for (;nowdist < maxdist;nowdist += incdist, ++thiscell) {
+ fnear = line_r - nowdist;
+ ffar = line_r + nowdist;
+ fcover = 1.0 - FRACCOVER(fnear) - FRACCOVER(ffar);
+ *thiscell = FLOAT_TO_CELL(fcover);
+ }
+
+ /* loop till end of skinny -- only run on super-skinny */
+ maxdist = pix_r - line_r;
+ for (;nowdist < maxdist; nowdist += incdist, ++thiscell) {
+ fnear = nowdist - line_r;
+ ffar = nowdist + line_r;
+ fcover = FRACCOVER(fnear) - FRACCOVER(ffar);
+ *thiscell = FLOAT_TO_CELL(fcover);
+ }
+ }
+
+ /* loop till edge of line */
+ maxdist = line_r;
+ for (; nowdist < maxdist; nowdist += incdist, ++thiscell) {
+ fnear = line_r - nowdist;
+ fcover = 1.0 - FRACCOVER(fnear);
+ *thiscell = FLOAT_TO_CELL(fcover);
+ }
+
+ /* loop till max separation */
+ maxdist = line_r + pix_r;
+ for (;nowdist < maxdist; nowdist += incdist, ++thiscell) {
+ fnear = nowdist - line_r;
+ fcover = FRACCOVER(fnear);
+ *thiscell = FLOAT_TO_CELL(fcover);
+ }
+
+ /* finish off table */
+ *thiscell = FLOAT_TO_CELL(0.0);
+ s->coverage[tablecells-1] = FLOAT_TO_CELL(0.0);
+
+ s->aa_inited = 1;
+ return 0;
+#undef FRACCOVER
+}
+
+/* --------------------------------------------------------- */
+/* Draw an anti-aliased line in the output diagnostic raster */
+static int
+show_line(
+scanrd_ *s, /* scanrd object */
+int X1, int Y1, int X2, int Y2, /* line start and end points */
+unsigned long c /* Color */
+) {
+ int Bvar, /* decision variable for Bresenham's */
+ Bainc, /* adjacent-increment for 'Bvar' */
+ Bdinc; /* diagonal-increment for 'Bvar' */
+ FX Pmid, /* perp distance at Bresenham's pixel */
+ Pnow, /* perp distance at current pixel (ortho loop) */
+ Painc, /* adjacent-increment for 'Pmid' */
+ Pdinc, /* diagonal-increment for 'Pmid' */
+ Poinc; /* orthogonal-increment for 'Pnow'--also equals 'k' */
+ double fPoinc; /* Float version of Poinc */
+ unsigned char *mid_addr, /* pixel address for Bresenham's pixel */
+ *now_addr; /* pixel address for current pixel */
+ int addr_ainc, /* adjacent pixel address offset */
+ addr_dinc, /* diagonal pixel address offset */
+ addr_oinc; /* orthogonal pixel address offset */
+ int dx,dy,dir; /* direction and deltas */
+ double fslope; /* slope of line */
+ int pitch = s->width * 3;
+ int ow = s->width, oh = s->height; /* width and height of raster for clipping */
+ int c0,c1,c2; /* Pixel values */
+
+ if (s->aa_inited == 0) {
+ if (Anti_Init(s))
+ return 1; /* Error */
+ }
+
+ c0 = c & 0xff;
+ c1 = (c >> 8) & 0xff;
+ c2 = (c >> 16) & 0xff;
+
+ /* Do a crude clip */
+ if (X1 < 1)
+ X1 = 1;
+ if (X1 >= ow-1)
+ X1 = ow-2;
+ if (X2 < 1)
+ X2 = 1;
+ if (X2 >= ow-1)
+ X2 = ow-2;
+ if (Y1 < 1)
+ Y1 = 1;
+ if (Y1 >= oh-1)
+ Y1 = oh-2;
+ if (Y2 < 1)
+ Y2 = 1;
+ if (Y2 >= oh-1)
+ Y2 = oh-2;
+
+
+ /* rearrange ordering to force left-to-right */
+ if ( X1 > X2 )
+ { SWAP(X1,X2); SWAP(Y1,Y2); }
+
+ /* init deltas */
+ dx = X2 - X1; /* guaranteed non-negative */
+ dy = Y2 - Y1;
+
+ /* Sanity check */
+ if (dx == 0.0 && dy == 0.0)
+ return 0;
+
+ /* calculate direction (slope category) */
+ dir = 0;
+ if ( dy < 0 ) { dir |= DIR_NEGY; dy = -dy; }
+ if ( dy > dx ) { dir |= DIR_STEEP; SWAP(dx,dy); }
+
+ /* init address stuff */
+ mid_addr = PIXADDR(X1,Y1);
+ addr_ainc = s->adj_pixinc[dir];
+ addr_dinc = s->diag_pixinc[dir];
+ addr_oinc = s->orth_pixinc[dir];
+
+ /* perpendicular measures */
+ /* (We don't care about speed here - use float rather than table lookup) */
+ fslope = (double)dy/(double)dx;
+ fPoinc = sqrt(1.0/(1.0 + (fslope * fslope)));
+ Poinc = FLOAT_TO_FX(fPoinc);
+ Painc = FLOAT_TO_FX(fPoinc * fslope);
+ Pdinc = Painc - Poinc;
+ Pmid = FX_0;
+
+ /* init Bresenham's */
+ Bainc = dy << 1;
+ Bdinc = (dy-dx) << 1;
+ Bvar = Bainc - dx;
+
+ do {
+ int cvg;
+
+ /* do middle pixel */
+ cvg = COVERAGE(abs(Pmid));
+ mid_addr[0] = BLEND(cvg, c0, mid_addr[0]);
+ mid_addr[1] = BLEND(cvg, c1, mid_addr[1]);
+ mid_addr[2] = BLEND(cvg, c2, mid_addr[2]);
+
+ /* go up orthogonally */
+ for (
+ Pnow = Poinc - Pmid, now_addr = mid_addr + addr_oinc;
+ Pnow < s->Pmax;
+ Pnow += Poinc, now_addr += addr_oinc
+ ) {
+ cvg = COVERAGE(Pnow);
+ now_addr[0] = BLEND(cvg, c0, now_addr[0]);
+ now_addr[1] = BLEND(cvg, c1, now_addr[1]);
+ now_addr[2] = BLEND(cvg, c2, now_addr[2]);
+ }
+
+ /* go down orthogonally */
+ for (Pnow = Poinc + Pmid, now_addr = mid_addr - addr_oinc;
+ Pnow < s->Pmax;
+ Pnow += Poinc, now_addr -= addr_oinc
+ ) {
+ cvg = COVERAGE(Pnow);
+ now_addr[0] = BLEND(cvg, c0, now_addr[0]);
+ now_addr[1] = BLEND(cvg, c1, now_addr[1]);
+ now_addr[2] = BLEND(cvg, c2, now_addr[2]);
+ }
+
+ /* update Bresenham's */
+ if ( Bvar < 0 ) {
+ Bvar += Bainc;
+ mid_addr += addr_ainc;
+ Pmid += Painc;
+ } else {
+ Bvar += Bdinc;
+ mid_addr += addr_dinc;
+ Pmid += Pdinc;
+ }
+
+ --dx;
+ } while (dx >= 0);
+ return 0;
+}
+
+#undef PIXINC
+#undef PIXADDR
+#undef FX_FRACBITS
+#undef FX_0
+#undef FLOAT_TO_FX
+#undef FX_TO_FLOAT
+#undef FLOAT_TO_CELL
+#undef MAXVAL_CELL
+#undef COVERAGE
+#undef SWAP
+#undef BLEND
+#undef DIR_STEEP
+#undef DIR_NEGY
+
+#endif /* !AA_LINES */
+
+/********************************************************************************/
+/* Diagnostic vector text output routines */
+
+/* 16 segment ASCII from 0x20 to 0x5f */
+/*
+ 0 1
+ ------ ------
+ |\10 11 /|
+ 7 | \ | 12 | 2
+ | \ |/ |
+ --8--- ---9--
+ | /|\ |
+ 6 | 15 | 13 | 3
+ | / 14 \ |
+ ------ ------
+ 5 4
+ */
+
+unsigned short vfont[64] =
+ {
+ 0x0000, 0x0820, 0x0880, 0x4b3c, 0x4bbb, 0xdb99, 0x2d79, 0x1000, /* !"#$%&' */
+ 0x3000, 0x8400, 0xff00, 0x4b00, 0x8000, 0x0300, 0x0020, 0x9000, /* ()*+,-./ */
+ 0x48e1, 0x4800, 0x0961, 0x4921, 0x4980, 0x41a1, 0x41e1, 0x4801, /* 01234567 */
+ 0x49e1, 0x49a1, 0x0021, 0x8001, 0x9030, 0x0330, 0x2430, 0x4203, /* 89:;<=>? */
+ 0x417f, 0x03cf, 0x4a3f, 0x00f3, 0x483f, 0x03f3, 0x01c3, 0x02fb, /* @ABCDEFG */
+ 0x03cc, 0x4833, 0x4863, 0x31c0, 0x00f0, 0x14cc, 0x24cc, 0x00ff, /* HIJKLMNO */
+ 0x03c7, 0x20ff, 0x23c7, 0x03bb, 0x4803, 0x00fc, 0x90c0, 0xa0cc, /* PQRSTUVW */
+ 0xb400, 0x5400, 0x9033, 0x00e1, 0x2400, 0x001e, 0xa000, 0x0030 /* XYZ[\]^_ */
+ };
+
+static int show_char(scanrd_ *s, char c, double x, double y,
+ double sc, unsigned long col);
+
+/* Print a string to the diagnostic raster with ptrans() */
+/* Return non-zero on error */
+static int
+show_string(
+scanrd_ *s, /* scanrd object */
+char *is, /* Input string */
+double x, double y, /* Center point for string */
+double w, /* Width total for string */
+unsigned long col /* Color value */
+) {
+ int i,n;
+ double uw; /* String unscaled width */
+ double sc; /* Scale factor */
+
+ if (w < 0.0)
+ w = -w;
+ n = strlen(is);
+ if (n == 0)
+ return 0;
+
+ /* Total unscaled width of the string */
+ uw = (n * 0.8 + (n >= 1 ? (n-1) * 0.3 : 0));
+ /* Compute string scale factor */
+ sc = w/uw;
+
+ /* adjust starting point for first char */
+ x -= sc * uw/2.0;
+ y -= sc * 0.5;
+
+ for (i = 0; i < n; i++) {
+ if (show_char(s,is[i],x,y,sc,col))
+ return 1;
+ x += sc * (0.8 + 0.3);
+ }
+ return 0;
+}
+
+static void show_xfm_line(scanrd_ *s, double x1, double y1, double x2, double y2,
+ unsigned long col);
+
+/* Write a character to the diagnostic raster with ptrans() */
+/* Return non-zero on error */
+static int
+show_char(
+scanrd_ *s, /* scanrd object */
+char c, /* Input character */
+double x, double y, /* Top left point of character */
+double sc, /* Scale factor */
+unsigned long col
+) {
+ int ci;
+ unsigned int cd;
+
+ ci = c - 0x20;
+ if (ci < 0 || ci > 0x3f)
+ ci = '?' - 0x20;
+ cd = vfont[ci];
+ /* Display each segment */
+ if (cd & 0x0001)
+ show_xfm_line(s, x,y,x+sc*0.4,y,col);
+ if (cd & 0x0002)
+ show_xfm_line(s, x+sc*0.4,y,x+sc*0.8,y,col);
+ if (cd & 0x0004)
+ show_xfm_line(s, x+sc*0.8,y,x+sc*0.8,y+sc*0.5,col);
+ if (cd & 0x0008)
+ show_xfm_line(s, x+sc*0.8,y+sc*0.5,x+sc*0.8,y+sc*1.0,col);
+ if (cd & 0x0010)
+ show_xfm_line(s, x+sc*0.8,y+sc*1.0,x+sc*0.4,y+sc*1.0,col);
+ if (cd & 0x0020)
+ show_xfm_line(s, x+sc*0.4,y+sc*1.0,x+0.0,y+sc*1.0,col);
+ if (cd & 0x0040)
+ show_xfm_line(s, x+0.0,y+sc*1.0,x+0.0,y+sc*0.5,col);
+ if (cd & 0x0080)
+ show_xfm_line(s, x+0.0,y+sc*0.5,x+0.0,y+0.0,col);
+ if (cd & 0x0100)
+ show_xfm_line(s, x+0.0,y+sc*0.5,x+sc*0.4,y+sc*0.5,col);
+ if (cd & 0x0200)
+ show_xfm_line(s, x+sc*0.4,y+sc*0.5,x+sc*0.8,y+sc*0.5,col);
+ if (cd & 0x0400)
+ show_xfm_line(s, x+0.0,y+0.0,x+sc*0.4,y+sc*0.5,col);
+ if (cd & 0x0800)
+ show_xfm_line(s, x+sc*0.4,y+0.0,x+sc*0.4,y+sc*0.5,col);
+ if (cd & 0x1000)
+ show_xfm_line(s, x+sc*0.8,y+0.0,x+sc*0.4,y+sc*0.5,col);
+ if (cd & 0x2000)
+ show_xfm_line(s, x+sc*0.8,y+sc*1.0,x+sc*0.4,y+sc*0.5,col);
+ if (cd & 0x4000)
+ show_xfm_line(s, x+sc*0.4,y+sc*1.0,x+sc*0.4,y+sc*0.5,col);
+ if (cd & 0x8000)
+ show_xfm_line(s, x+0.0,y+sc*1.0,x+sc*0.4,y+sc*0.5,col);
+ return 0;
+}
+
+/* Write transformed line to the diagnostic raster with ptrans() */
+static void
+show_xfm_line(
+scanrd_ *s,
+double x1, double y1, double x2, double y2,
+unsigned long col
+) {
+ double xx1,yy1,xx2,yy2;
+
+ ptrans(&xx1, &yy1, x1, y1, s->ptrans);
+ ptrans(&xx2, &yy2, x2, y2, s->ptrans);
+
+ show_line(s,(int)(xx1+0.5),(int)(yy1+0.5),(int)(xx2+0.5),(int)(yy2+0.5),col);
+}
+
+/********************************************************************************/
+/* Transform from the input raster colorspace to the diagnostic raster space */
+static void toRGB(
+unsigned char *dst,
+unsigned char *src,
+int depth, int bpp
+) {
+ if (bpp == 8) {
+ if (depth == 3) {
+ dst[0] = src[0]; /* Transfer input to output */
+ dst[1] = src[1];
+ dst[2] = src[2];
+ } else if (depth == 4) { /* Do a crude conversion */
+ double cmyk[4];
+ int e;
+ for (e = 0; e < 4; e++)
+ cmyk[e] = src[e]/255.0;
+ for (e = 0; e < 3; e++) {
+ cmyk[e] = cmyk[e] * 0.7 + 0.3 * cmyk[3];
+ if (cmyk[e] < cmyk[3])
+ cmyk[e] = cmyk[3];
+ dst[e] = 255 - (int)(cmyk[e] * 255.0 + 0.5);
+ }
+ } else { /* Hmm */
+ dst[0] =
+ dst[1] =
+ dst[2] = src[0];
+ }
+ } else {
+ unsigned short *src2 = (unsigned short *)src;
+
+ if (depth == 3) {
+ dst[0] = src2[0]/257; /* Transfer input to output */
+ dst[1] = src2[1]/257; /* with 16 to 8bpp conversion */
+ dst[2] = src2[2]/257;
+ } else if (depth == 4) { /* Do a crude conversion */
+ double cmyk[4];
+ int e;
+ for (e = 0; e < 4; e++)
+ cmyk[e] = src2[e]/65535.0;
+ for (e = 0; e < 3; e++) {
+ cmyk[e] = cmyk[e] * 0.7 + 0.3 * cmyk[3];
+ if (cmyk[e] < cmyk[3])
+ cmyk[e] = cmyk[3];
+ dst[e] = 255 - (int)(cmyk[e] * 255.0 + 0.5);
+ }
+ } else { /* Hmm */
+ dst[0] =
+ dst[1] =
+ dst[2] = src2[0]/257;
+ }
+ }
+}
+
+
+/* Convert from XYZ scale 100 to Lab D50 */
+static void XYZ2Lab(double *out, double *in) {
+ double X = in[0], Y = in[1], Z = in[2];
+ double x,y,z,fx,fy,fz;
+
+ x = X/96.42;
+ y = Y/100.0;
+ z = Z/82.49;
+
+ if (x > 0.008856451586)
+ fx = pow(x,1.0/3.0);
+ else
+ fx = 7.787036979 * x + 16.0/116.0;
+
+ if (y > 0.008856451586)
+ fy = pow(y,1.0/3.0);
+ else
+ fy = 7.787036979 * y + 16.0/116.0;
+
+ if (z > 0.008856451586)
+ fz = pow(z,1.0/3.0);
+ else
+ fz = 7.787036979 * z + 16.0/116.0;
+
+ out[0] = 116.0 * fy - 16.0;
+ out[1] = 500.0 * (fx - fy);
+ out[2] = 200.0 * (fy - fz);
+}
+
+/* Convert from a scanned pixel value to an aproximate Lab value */
+static void pval2Lab(double *out, double *in, int depth) {
+ double wXYZ[3];
+ double XYZ[3];
+ int e, j;
+
+ if (depth == 3) { /* Assume RGB */
+
+ double clrnts[3][3] = { /* Red, Green & Blue XYZ values */
+ { 0.412414, 0.212642, 0.019325 },
+ { 0.357618, 0.715136, 0.119207 },
+ { 0.180511, 0.072193, 0.950770 }
+ };
+
+ wXYZ[0] = 0.950543; /* Because we're using sRGB primaries */
+ wXYZ[1] = 1.0; /* the white point is D65 */
+ wXYZ[2] = 1.089303;
+
+ XYZ[0] = XYZ[1] = XYZ[2] = 0.0;
+
+ for (e = 0; e < 3; e++) {
+ double v = in[e]/255.0;
+
+ if (v < 0.0)
+ v = 0.0;
+ else if (v > 1.0)
+ v = 1.0;
+ if (v <= 0.03928)
+ v /= 12.92;
+ else
+ v = pow((0.055 + v)/1.055, 2.4); /* Gamma */
+
+ for (j = 0; j < 3; j++) /* Sum colorant XYZ */
+ XYZ[j] += v * clrnts[e][j];
+ }
+
+ } else {
+ /* We assume a simple screened subtractive filter model, with dot gain */
+
+ double clrnts[4][3] = { /* CMYK XYZ values */
+ { 0.12, 0.18, 0.48 },
+ { 0.38, 0.19, 0.20 },
+ { 0.76, 0.81, 0.11 },
+ { 0.04, 0.04, 0.04 }
+ };
+
+ /* start with white */
+ XYZ[0] = wXYZ[0] = 0.9642;
+ XYZ[1] = wXYZ[1] = 1.0;
+ XYZ[2] = wXYZ[2] = 0.8249;
+
+ /* And filter it out for each component */
+ for (e = 0; e < 4; e++) {
+ double v = in[e]/255.0;
+
+ if (v < 0.0)
+ v = 0.0;
+ else if (v > 1.0)
+ v = 1.0;
+ v = 1.0 - pow(1.0 - v, 2.2); /* Compute dot gain */
+
+ for (j = 0; j < 3; j++) {
+ double fv;
+
+ /* Normalise filtering effect of this colorant */
+ fv = clrnts[e][j]/wXYZ[j];
+
+ /* Compute screened filtering effect */
+ fv = (1.0 - v) + v * fv;
+
+ /* Apply filter to our current value */
+ XYZ[j] *= fv;
+ }
+ }
+ }
+
+ /* Convert to Lab */
+ {
+ double X = XYZ[0], Y = XYZ[1], Z = XYZ[2];
+ double x,y,z,fx,fy,fz;
+
+ x = X/wXYZ[0];
+ y = Y/wXYZ[1];
+ z = Z/wXYZ[2];
+
+ if (x > 0.008856451586)
+ fx = pow(x,1.0/3.0);
+ else
+ fx = 7.787036979 * x + 16.0/116.0;
+
+ if (y > 0.008856451586)
+ fy = pow(y,1.0/3.0);
+ else
+ fy = 7.787036979 * y + 16.0/116.0;
+
+ if (z > 0.008856451586)
+ fz = pow(z,1.0/3.0);
+ else
+ fz = 7.787036979 * z + 16.0/116.0;
+
+ out[0] = 116.0 * fy - 16.0;
+ out[1] = 500.0 * (fx - fy);
+ out[2] = 200.0 * (fy - fz);
+ }
+}
+
+/********************************************************************************/
+
+static int
+scanrd_write_diag(scanrd_ *s) {
+ int y;
+ unsigned char *op;
+ int stride = 3 * s->width;
+
+ if ((s->flags & SI_SHOW_FLAGS) == 0 || s->write_line == NULL)
+ return 0;
+
+ /* Write out the tiff file */
+ for (op = s->out, y = 0; y < s->height; ++y, op += stride) {
+ if (s->write_line(s->ddata, y, (char *)op)) {
+ s->errv = SI_DIAG_WRITE_ERR;
+ sprintf(s->errm,"scanrd: write_line() returned error");
+ return 1;
+ }
+ }
+ return 0;
+}
+