From 59283bdf6d439fdcd1a0ab07b318b48031091b34 Mon Sep 17 00:00:00 2001
From: Andreas Gal <andreas.gal@gmail.com>
Date: Mon, 30 Jan 2012 01:38:28 -0800
Subject: [PATCH 1/2] implement sampled functions based on the PDF spec

---
 src/function.js | 126 ++++++++++++++++--------------------------------
 1 file changed, 42 insertions(+), 84 deletions(-)

diff --git a/src/function.js b/src/function.js
index 26b8fe679..3093b0a99 100644
--- a/src/function.js
+++ b/src/function.js
@@ -125,109 +125,67 @@ var PDFFunction = (function PDFFunctionClosure() {
       else
         decode = toMultiArray(decode);
 
-      // Precalc the multipliers
-      var inputMul = new Float64Array(inputSize);
-      for (var i = 0; i < inputSize; ++i) {
-        inputMul[i] = (encode[i][1] - encode[i][0]) /
-                  (domain[i][1] - domain[i][0]);
-      }
-
-      var idxMul = new Int32Array(inputSize);
-      idxMul[0] = outputSize;
-      for (i = 1; i < inputSize; ++i) {
-        idxMul[i] = idxMul[i - 1] * size[i - 1];
-      }
-
-      var nSamples = outputSize;
-      for (i = 0; i < inputSize; ++i)
-          nSamples *= size[i];
-
       var samples = this.getSampleArray(size, outputSize, bps, str);
 
       return [
         CONSTRUCT_SAMPLED, inputSize, domain, encode, decode, samples, size,
-        outputSize, bps, range, inputMul, idxMul, nSamples
+        outputSize, Math.pow(2, bps) - 1, range
       ];
     },
 
     constructSampledFromIR: function pdfFunctionConstructSampledFromIR(IR) {
-      var inputSize = IR[1];
-      var domain = IR[2];
-      var encode = IR[3];
-      var decode = IR[4];
-      var samples = IR[5];
-      var size = IR[6];
-      var outputSize = IR[7];
-      var bps = IR[8];
-      var range = IR[9];
-      var inputMul = IR[10];
-      var idxMul = IR[11];
-      var nSamples = IR[12];
+      // See chapter 3, page 109 of the PDF reference
+      function interpolate(x, xmin, xmax, ymin, ymax) {
+        return ymin + ((x - xmin) * ((ymax - ymin) / (xmax - xmin)));
+      }
 
       return function constructSampledFromIRResult(args) {
-        if (inputSize != args.length)
+        // See chapter 3, page 110 of the PDF reference.
+        var m = IR[1];
+        var domain = IR[2];
+        var encode = IR[3];
+        var decode = IR[4];
+        var samples = IR[5];
+        var size = IR[6];
+        var n = IR[7];
+        var mask = IR[8];
+        var range = IR[9];
+
+        if (m != args.length)
           error('Incorrect number of arguments: ' + inputSize + ' != ' +
                 args.length);
-        // Most of the below is a port of Poppler's implementation.
-        // TODO: There's a few other ways to do multilinear interpolation such
-        // as piecewise, which is much faster but an approximation.
-        var out = new Float64Array(outputSize);
-        var x;
-        var e = new Array(inputSize);
-        var efrac0 = new Float64Array(inputSize);
-        var efrac1 = new Float64Array(inputSize);
-        var sBuf = new Float64Array(1 << inputSize);
-        var i, j, k, idx, t;
 
-        // map input values into sample array
-        for (i = 0; i < inputSize; ++i) {
-          x = (args[i] - domain[i][0]) * inputMul[i] + encode[i][0];
-          if (x < 0) {
-            x = 0;
-          } else if (x > size[i] - 1) {
-            x = size[i] - 1;
-          }
-          e[i] = [Math.floor(x), 0];
-          if ((e[i][1] = e[i][0] + 1) >= size[i]) {
-            // this happens if in[i] = domain[i][1]
-            e[i][1] = e[i][0];
-          }
-          efrac1[i] = x - e[i][0];
-          efrac0[i] = 1 - efrac1[i];
-        }
+        var x = args;
+        var y = new Float64Array(n * m);
 
-        // for each output, do m-linear interpolation
-        for (i = 0; i < outputSize; ++i) {
+        // Map x_i to y_j for 0 <= i < m using the sampled function.
+        for (var i = 0; i < m; ++i) {
+          // x_i' = min(max(x_i, Domain_2i), Domain_2i+1)
+          var domain_2i = domain[2 * i];
+          var domain_2i_1 = domain[2 * i + 1];
+          var xi = Math.min(Math.max(x[i], domain_2i), domain_2i_1);
 
-          // pull 2^m values out of the sample array
-          for (j = 0; j < (1 << inputSize); ++j) {
-            idx = i;
-            for (k = 0, t = j; k < inputSize; ++k, t >>= 1) {
-              idx += idxMul[k] * (e[k][t & 1]);
-            }
-            if (idx >= 0 && idx < nSamples) {
-              sBuf[j] = samples[idx];
-            } else {
-              sBuf[j] = 0; // TODO Investigate if this is what Adobe does
-            }
-          }
+          // e_i = Interpolate(x_i', Domain_2i, Domain_2i+1, Encode_2i, Encode_2i+1)
+          var e = interpolate(xi, domain_2i, domain_2i_1, encode[2 * i], encode[2 * i + 1]);
 
-          // do m sets of interpolations
-          for (j = 0, t = (1 << inputSize); j < inputSize; ++j, t >>= 1) {
-            for (k = 0; k < t; k += 2) {
-              sBuf[k >> 1] = efrac0[j] * sBuf[k] + efrac1[j] * sBuf[k + 1];
-            }
-          }
+          // e_i' = min(max(e_i, 0), Size_i - 1)
+          e = Math.min(Math.max(e, 0), size[i] - 1);
 
-          // map output value to range
-          out[i] = (sBuf[0] * (decode[i][1] - decode[i][0]) + decode[i][0]);
-          if (out[i] < range[i][0]) {
-            out[i] = range[i][0];
-          } else if (out[i] > range[i][1]) {
-            out[i] = range[i][1];
+          var in = i * n;
+
+          for (var j = 0; j < n; ++j) {
+            // average the two nearest neighbors in the sampling table
+            var rj = (samples[Math.floor(e) * n + j] + samples[Math.ceil(e) * n + j]) / 2;
+
+            // r_j' = Interpolate(r_j, 0, 2^BitsPerSample - 1, Decode_2j, Decode_2j+1)
+            rj = interpolate(rj, 0, mask, 1, decode[2 * j], decode[2 * j + 1]);
+
+            // y_j = min(max(r_j, range_2j, range_2j+1)
+            y[in + j] = Math.min(Math.max(rj, range[2 * j], range[2 * j + 1]));
           }
         }
-        return out;
+
+        return y;
       }
     },
 

From 5fba376a336e57e5961a239f387d99b1cfb11c71 Mon Sep 17 00:00:00 2001
From: notmasteryet <async.processingjs@yahoo.com>
Date: Tue, 31 Jan 2012 20:19:44 -0600
Subject: [PATCH 2/2] Fixing interpolation (continuation of #1143)

---
 src/function.js | 66 ++++++++++++++++++++++++++++++++++++-------------
 1 file changed, 49 insertions(+), 17 deletions(-)

diff --git a/src/function.js b/src/function.js
index 3093b0a99..4f81158f0 100644
--- a/src/function.js
+++ b/src/function.js
@@ -156,33 +156,65 @@ var PDFFunction = (function PDFFunctionClosure() {
                 args.length);
 
         var x = args;
-        var y = new Float64Array(n * m);
 
+        // Building the cube vertices: its part and sample index
+        // http://rjwagner49.com/Mathematics/Interpolation.pdf
+        var cubeVertices = 1 << m;
+        var cubeN = new Float64Array(cubeVertices);
+        var cubeVertex = new Uint32Array(cubeVertices);
+        for (var j = 0; j < cubeVertices; j++)
+          cubeN[j] = 1;
+
+        var k = n, pos = 1;
         // Map x_i to y_j for 0 <= i < m using the sampled function.
         for (var i = 0; i < m; ++i) {
           // x_i' = min(max(x_i, Domain_2i), Domain_2i+1)
-          var domain_2i = domain[2 * i];
-          var domain_2i_1 = domain[2 * i + 1];
+          var domain_2i = domain[i][0];
+          var domain_2i_1 = domain[i][1];
           var xi = Math.min(Math.max(x[i], domain_2i), domain_2i_1);
 
-          // e_i = Interpolate(x_i', Domain_2i, Domain_2i+1, Encode_2i, Encode_2i+1)
-          var e = interpolate(xi, domain_2i, domain_2i_1, encode[2 * i], encode[2 * i + 1]);
+          // e_i = Interpolate(x_i', Domain_2i, Domain_2i+1,
+          //                   Encode_2i, Encode_2i+1)
+          var e = interpolate(xi, domain_2i, domain_2i_1,
+                              encode[i][0], encode[i][1]);
 
           // e_i' = min(max(e_i, 0), Size_i - 1)
-          e = Math.min(Math.max(e, 0), size[i] - 1);
+          var size_i = size[i];
+          e = Math.min(Math.max(e, 0), size_i - 1);
 
-          var in = i * n;
-
-          for (var j = 0; j < n; ++j) {
-            // average the two nearest neighbors in the sampling table
-            var rj = (samples[Math.floor(e) * n + j] + samples[Math.ceil(e) * n + j]) / 2;
-
-            // r_j' = Interpolate(r_j, 0, 2^BitsPerSample - 1, Decode_2j, Decode_2j+1)
-            rj = interpolate(rj, 0, mask, 1, decode[2 * j], decode[2 * j + 1]);
-
-            // y_j = min(max(r_j, range_2j, range_2j+1)
-            y[in + j] = Math.min(Math.max(rj, range[2 * j], range[2 * j + 1]));
+          // Adjusting the cube: N and vertex sample index
+          var e0 = e < size_i - 1 ? Math.floor(e) : e - 1; // e1 = e0 + 1;
+          var n0 = e0 + 1 - e; // (e1 - e) / (e1 - e0);
+          var n1 = e - e0; // (e - e0) / (e1 - e0);
+          var offset0 = e0 * k;
+          var offset1 = offset0 + k; // e1 * k
+          for (var j = 0; j < cubeVertices; j++) {
+            if (j & pos) {
+              cubeN[j] *= n1;
+              cubeVertex[j] += offset1;
+            } else {
+              cubeN[j] *= n0;
+              cubeVertex[j] += offset0;
+            }
           }
+
+          k *= size_i;
+          pos <<= 1;
+        }
+
+        var y = new Float64Array(n);
+        for (var j = 0; j < n; ++j) {
+          // Sum all cube vertices' samples portions
+          var rj = 0;
+          for (var i = 0; i < cubeVertices; i++)
+            rj += samples[cubeVertex[i] + j] * cubeN[i];
+
+          // r_j' = Interpolate(r_j, 0, 2^BitsPerSample - 1,
+          //                    Decode_2j, Decode_2j+1)
+          rj = interpolate(rj, 0, 1, decode[j][0], decode[j][1]);
+
+          // y_j = min(max(r_j, range_2j), range_2j+1)
+          y[j] = Math.min(Math.max(rj, range[j][0]), range[j][1]);
         }
 
         return y;