Tweaks to pattern scaling.

Spun out of bug 706373, though not believed to be a fix for it.

When we calculate pattern matrixes and bboxes etc, floating point
accuracy conspires to work against us, in particular where we
end up mapping boxes back and forth via inverse matrices.

The rounding is such that we can end up expanding the number of
repeats that we have to render.

We try to be smarter about that (or at least to give saner bboxes
and matrices that are easier to debug) by employing some smart
("sane") rounding.

2 files changed, 141 insertions, 11 deletions

diff --git a/base/gsptype1.c b/base/gsptype1.c
index 5353da27a..cbc85d80f 100644
--- a/base/gsptype1.c
+++ b/base/gsptype1.c
@@ -48,6 +48,9 @@
 #define ADJUST_SCALE_BY_GS_TRADITION 0  /* Old code = 1 */
 #define ADJUST_AS_ADOBE 1               /* Old code = 0 *//* This one is closer to Adobe. */
 
+#define fastfloor(x) (((int)(x)) - (((x)<0) && ((x) != (float)(int)(x))))
+#define fastfloord(x) (((int)(x)) - (((x)<0) && ((x) != (double)(int)(x))))
+
 /* GC descriptors */
 private_st_pattern1_template();
 public_st_pattern1_instance();
@@ -130,6 +133,72 @@ gs_makepattern(gs_client_color * pcc, const gs_pattern1_template_t * pcp,
     return gs_pattern1_make_pattern(pcc, (const gs_pattern_template_t *)pcp,
                                     pmat, pgs, mem);
 }
+
+/* Limited accuracy due to Floating point implementation limits can
+ * cause us headaches due to values not being representable.
+ * This is particular bad when we start using matrix maths (and
+ * inverse matrixes in particular) to map bboxes into device space.
+ * We therefore have a set of functions to 'sanely' round stuff.
+ * Essentially, any device space location that is sufficiently
+ * close to an exact pixel boundary will be clamped to that boundary.
+ */
+#define SANE_THRESHOLD 0.001f
+
+static int
+sane_ceil(float f)
+{
+    int i = (int)f;
+
+    if (f - i < SANE_THRESHOLD)
+        return i;
+    return i+1;
+}
+static float
+sane_clamp_float(float f)
+{
+    int i = (int)fastfloor(f);
+
+    if (f - i < SANE_THRESHOLD)
+        return (float)i;
+    else if (f - i > (1-SANE_THRESHOLD))
+        return (float)(i+1);
+    return f;
+}
+static double
+sane_clamp_double(double d)
+{
+    int i = (int)fastfloord(d);
+
+    if (d - i < SANE_THRESHOLD)
+        return (double)i;
+    else if (d - i > (1-SANE_THRESHOLD))
+        return (double)(i+1);
+    return d;
+}
+static void
+sane_clamp_rect(gs_rect *r)
+{
+    double x0 = sane_clamp_double(r->p.x);
+    double y0 = sane_clamp_double(r->p.y);
+    double x1 = sane_clamp_double(r->q.x);
+    double y1 = sane_clamp_double(r->q.y);
+
+    /* Be careful not to round stuff to zero, because this breaks fts_15_1529.pdf */
+    if (x0 != x1)
+        r->p.x = x0, r->q.x = x1;
+    if (y0 != y1)
+        r->p.y = y0, r->q.y = y1;
+}
+static void
+sane_clamp_matrix(gs_matrix *m)
+{
+    m->xx = sane_clamp_float(m->xx);
+    m->xy = sane_clamp_float(m->xy);
+    m->yx = sane_clamp_float(m->yx);
+    m->yy = sane_clamp_float(m->yy);
+    m->tx = sane_clamp_float(m->tx);
+    m->ty = sane_clamp_float(m->ty);
+}
 static int
 gs_pattern1_make_pattern(gs_client_color * pcc,
                          const gs_pattern_template_t * ptemp,
@@ -219,8 +288,8 @@ gs_pattern1_make_pattern(gs_client_color * pcc,
     } else if (inst.templat.TilingType == 2) {
         /* Always round up for TilingType 2, as we don't want any
          * content to be lost. */
-        inst.size.x = (int)ceil(bbw);
-        inst.size.y = (int)ceil(bbh);
+        inst.size.x = sane_ceil(bbw);
+        inst.size.y = sane_ceil(bbh);
     } else {
         /* For TilingType's other than 2 allow us to round up or down
          * to whatever is nearer. The scale we do later prevents us
@@ -471,6 +540,9 @@ clamp_pattern_bbox(gs_pattern1_instance_t * pinst, gs_rect * pbbox,
     code = gs_bbox_transform_inverse(&dev_page, pmat, &pat_page);
     if (code < 0)
         return code;
+    /* So pat_page is the bbox for the region in pattern space that will
+     * be mapped forwards to cover the page. We want to find which tiles
+     * are required to cover this area. */
     /*
      * Determine the location of the pattern origin in device coordinates.
      */
@@ -479,6 +551,15 @@ clamp_pattern_bbox(gs_pattern1_instance_t * pinst, gs_rect * pbbox,
      * Determine our starting point.  We start with a postion that puts the
      * pattern below and to the left of the page (in pattern space) and scan
      * until the pattern is above and right of the page.
+     *
+     * So the right hand edge of each tile is:
+     *
+     *  xstep * n + pinst->templat.BBox.q.x
+     *
+     * and we want the largest n s.t. that is <= pat_page.p.x. i.e.
+     *
+     *  xstep * n <= pat_page.p.x - pinst->templat.BBox.q.x < xstep *n+1
+     *  n <= (pat_page.p.x - pinst->templat.BBox.q.x) / xstep < n+1
      */
     ixpat = (int) floor((pat_page.p.x - pinst->templat.BBox.q.x) / xstep);
     iystart = (int) floor((pat_page.p.y - pinst->templat.BBox.q.y) / ystep);
@@ -605,6 +686,8 @@ compute_inst_matrix(gs_pattern1_instance_t * pinst,
     if (code < 0)
         return code;
 
+    sane_clamp_rect(pbbox);
+
     *pbbw = pbbox->q.x - pbbox->p.x;
     *pbbh = pbbox->q.y - pbbox->p.y;
 
@@ -662,6 +745,8 @@ compute_inst_matrix(gs_pattern1_instance_t * pinst,
     pinst->step_matrix.yx = yx;
     pinst->step_matrix.yy = yy;
 
+    sane_clamp_matrix(&pinst->step_matrix);
+
     /*
      * Some applications produce patterns that are larger than the page.
      * If the bounding box for the pattern is larger than the page. clamp
diff --git a/base/gxp1fill.c b/base/gxp1fill.c
index 519e853cf..9fadb3057 100644
--- a/base/gxp1fill.c
+++ b/base/gxp1fill.c
@@ -1,4 +1,4 @@
-/* Copyright (C) 2001-2022 Artifex Software, Inc.
+/* Copyright (C) 2001-2023 Artifex Software, Inc.
    All Rights Reserved.
 
    This software is provided AS-IS with no warranty, either express or
@@ -140,6 +140,45 @@ tile_fill_init(tile_fill_state_t * ptfs, const gx_device_color * pdevc,
     return tile_clip_initialize(ptfs->cdev, ptfs->tmask, dev, px, py);
 }
 
+static int
+threshold_ceil(float f, float e)
+{
+    int i = (int)fastfloor(f);
+
+    if (f - i < e)
+        return i;
+    return i+1;
+}
+static int
+threshold_floor(float f, float e)
+{
+    int i = (int)fastfloor(f);
+
+    if (f - i > 1-e)
+        return i+1;
+    return i;
+}
+
+/* Return a conservative approximation to the maximum expansion caused by
+ * a given matrix. */
+static float
+matrix_expansion(gs_matrix *m)
+{
+    float e, f;
+
+    e = fabs(m->xx);
+    f = fabs(m->xy);
+    if (f > e)
+        e = f;
+    f = fabs(m->yx);
+    if (f > e)
+        e = f;
+    f = fabs(m->yy);
+    if (f > e)
+        e = f;
+
+    return e*e;
+}
 /*
  * Fill with non-standard X and Y stepping.
  * ptile is pdevc->colors.pattern.{m,p}_tile.
@@ -173,6 +212,12 @@ tile_by_steps(tile_fill_state_t * ptfs, int x0, int y0, int w0, int h0,
         double bbw = ptile->bbox.q.x - ptile->bbox.p.x;
         double bbh = ptile->bbox.q.y - ptile->bbox.p.y;
         double u0, v0, u1, v1;
+        /* Any difference smaller than error is guaranteed to result in
+         * less than a pixels difference in the output. Use this as a
+         * threshold when rounding to allow for inaccuracies in
+         * floating point maths. This enables us to avoid doing more
+         * repeats than we need to. */
+        float error = 1/matrix_expansion(&step_matrix);
 
         bbox.p.x = x0, bbox.p.y = y0;
         bbox.q.x = x1, bbox.q.y = y1;
@@ -190,16 +235,16 @@ tile_by_steps(tile_fill_state_t * ptfs, int x0, int y0, int w0, int h0,
          * each pixel of the rectangle being filled with *every* pattern
          * that overlaps it, not just *some* pattern copy.
          */
-        u0 = ibbox.p.x - max(ptile->bbox.p.x, 0) - 0.000001;
-        v0 = ibbox.p.y - max(ptile->bbox.p.y, 0) - 0.000001;
-        u1 = ibbox.q.x - min(ptile->bbox.q.x, 0) + 0.000001;
-        v1 = ibbox.q.y - min(ptile->bbox.q.y, 0) + 0.000001;
+        u0 = ibbox.p.x - max(ptile->bbox.p.x, 0);
+        v0 = ibbox.p.y - max(ptile->bbox.p.y, 0);
+        u1 = ibbox.q.x - min(ptile->bbox.q.x, 0);
+        v1 = ibbox.q.y - min(ptile->bbox.q.y, 0);
         if (!ptile->is_simple)
             u0 -= bbw, v0 -= bbh, u1 += bbw, v1 += bbh;
-        i0 = (int)fastfloor(u0);
-        j0 = (int)fastfloor(v0);
-        i1 = (int)ceil(u1);
-        j1 = (int)ceil(v1);
+        i0 = threshold_floor(u0, error);
+        j0 = threshold_floor(v0, error);
+        i1 = threshold_ceil(u1, error);
+        j1 = threshold_ceil(v1, error);
     }
     if_debug4m('T', mem, "[T]i=(%d,%d) j=(%d,%d)\n", i0, i1, j0, j1);
     for (i = i0; i < i1; i++)