4e497185814d58d9e20c6987b0da71b789bd46ca
[jalview.git] / src2 / javajs / img / JpgEncoder.java
1 // Version 1.0a
2 // Copyright (C) 1998, James R. Weeks and BioElectroMech.
3 // Visit BioElectroMech at www.obrador.com.  Email James@obrador.com.
4
5 // See license.txt for details about the allowed used of this software.
6 // This software is based in part on the work of the Independent JPEG Group.
7 // See IJGreadme.txt for details about the Independent JPEG Group's license.
8
9 // This encoder is inspired by the Java Jpeg encoder by Florian Raemy,
10 // studwww.eurecom.fr/~raemy.
11 // It borrows a great deal of code and structure from the Independent
12 // Jpeg Group's Jpeg 6a library, Copyright Thomas G. Lane.
13 // See license.txt for details 
14
15 /*
16  * JpegEncoder and its associated classes are Copyright (c) 1998, James R. Weeks and BioElectroMech
17  * see(Jmol/src/com/obrador/license.txt)
18  * 
19  * javjs.img.JpegEncoder.java was adapted by Bob Hanson
20  * 
21  * for Jmol in the following ways:
22  * 
23  * 1) minor coding efficiencies were made in some for() loops.
24  * 2) methods not used by Jmol were commented out
25  * 3) method and variable signatures were modified to provide 
26  *    more appropriate method privacy.
27  * 4) additions for Java2Script compatibility 
28  * 
29  * Original files are maintained in the Jmol.src.com.obrador package, but
30  * these original files are not distributed with Jmol.
31  *   
32 */
33
34 package javajs.img;
35
36 import java.io.IOException;
37 import java.util.Map;
38
39 import javajs.img.ImageEncoder;
40 import javajs.util.AU;
41 import javajs.util.OC;
42
43 /**
44  * JpegEncoder - The JPEG main program which performs a jpeg compression of an
45  * image.
46  * 
47  *  A system to allow the full Jmol state -- regardless of length -- 
48  *  to be encoded in a set of APP1 (FFE1) tags.
49  *  But we have to be careful about line ends for backward compatibility. 
50  *  This solution is not 100% effective, because some data lines may in principle be 
51  *  Very large and may not contain new lines for more than 65500 characters, 
52  *  But that would be very unusual. Perhaps a huge data set loaded from a 
53  *  string. Introduced in Jmol 12.1.36. Bob Hanson
54  *  
55  * See org.com.obrador.license.txt
56  * 
57  */
58
59 public class JpgEncoder extends ImageEncoder {
60
61   // this string will GENERALLY appear at the end of lines and be escaped 
62   private static final int CONTINUE_MAX = 65500; // some room to spare here. 
63   private static final int CONTINUE_MAX_BUFFER = CONTINUE_MAX + 10; // never break up last 10 bytes
64
65   private JpegObj jpegObj;
66   private Huffman huf;
67   private DCT dct;
68   protected int defaultQuality = 100;
69   private String applicationTag;
70
71   public JpgEncoder() {
72
73   }
74
75   @Override
76   protected void setParams(Map<String, Object> params) {
77     if (quality <= 0)
78       quality = (params.containsKey("qualityJPG") ? ((Integer) params.get("qualityJPG")).intValue() : defaultQuality);
79     jpegObj = new JpegObj();
80     jpegObj.comment = (String) params.get("comment");
81     applicationTag = (String) params.get("jpgAppTag");
82   }
83
84   @Override
85   protected void generate() throws IOException {
86     jpegObj.imageWidth = width;
87     jpegObj.imageHeight = height;
88     dct = new DCT(quality);
89     huf = new Huffman(width, height);
90     if (jpegObj == null)
91       return;
92     jpegObj.getYCCArray(pixels);
93     String longState = writeHeaders(jpegObj, dct);
94     writeCompressedData(jpegObj, dct, huf);
95     writeMarker(eoi);
96     if (longState != null) {
97       byte[] b = longState.getBytes();
98       out.write(b, 0, b.length);
99     }
100   }
101
102   private void writeCompressedData(JpegObj jpegObj, DCT dct, Huffman huf) {
103     int i, j, r, c, a, b;
104     int comp, xpos, ypos, xblockoffset, yblockoffset;
105     float inputArray[][];
106     float dctArray1[][] = new float[8][8];
107     double dctArray2[][] = new double[8][8];
108     int dctArray3[] = new int[8 * 8];
109
110     /*
111      * This method controls the compression of the image.
112      * Starting at the upper left of the image, it compresses 8x8 blocks
113      * of data until the entire image has been compressed.
114      */
115
116     int lastDCvalue[] = new int[jpegObj.numberOfComponents];
117     //int zeroArray[] = new int[64]; // initialized to hold all zeros
118     //int Width = 0, Height = 0;
119     //int nothing = 0, not;
120     int minBlockWidth, minBlockHeight;
121     // This initial setting of MinBlockWidth and MinBlockHeight is done to
122     // ensure they start with values larger than will actually be the case.
123     minBlockWidth = ((huf.imageWidth % 8 != 0) ? (int) (Math
124         .floor(huf.imageWidth / 8.0) + 1) * 8 : huf.imageWidth);
125     minBlockHeight = ((huf.imageHeight % 8 != 0) ? (int) (Math
126         .floor(huf.imageHeight / 8.0) + 1) * 8 : huf.imageHeight);
127     for (comp = 0; comp < jpegObj.numberOfComponents; comp++) {
128       minBlockWidth = Math.min(minBlockWidth, jpegObj.blockWidth[comp]);
129       minBlockHeight = Math.min(minBlockHeight, jpegObj.blockHeight[comp]);
130     }
131     xpos = 0;
132     for (r = 0; r < minBlockHeight; r++) {
133       for (c = 0; c < minBlockWidth; c++) {
134         xpos = c * 8;
135         ypos = r * 8;
136         for (comp = 0; comp < jpegObj.numberOfComponents; comp++) {
137           //Width = JpegObj.BlockWidth[comp];
138           //Height = JpegObj.BlockHeight[comp];
139           inputArray = jpegObj.components[comp];
140           int vsampF = jpegObj.vsampFactor[comp];
141           int hsampF = jpegObj.hsampFactor[comp];
142           int qNumber = jpegObj.qtableNumber[comp];
143           int dcNumber = jpegObj.dctableNumber[comp];
144           int acNumber = jpegObj.actableNumber[comp];
145
146           for (i = 0; i < vsampF; i++) {
147             for (j = 0; j < hsampF; j++) {
148               xblockoffset = j * 8;
149               yblockoffset = i * 8;
150               for (a = 0; a < 8; a++) {
151                 for (b = 0; b < 8; b++) {
152
153                   // I believe this is where the dirty line at the bottom of
154                   // the image is coming from.
155                   // I need to do a check here to make sure I'm not reading past
156                   // image data.
157                   // This seems to not be a big issue right now. (04/04/98)
158
159                   dctArray1[a][b] = inputArray[ypos + yblockoffset + a][xpos
160                       + xblockoffset + b];
161                 }
162               }
163               // The following code commented out because on some images this technique
164               // results in poor right and bottom borders.
165               // if ((!JpegObj.lastColumnIsDummy[comp] || c < Width - 1) &&
166               //       (!JpegObj.lastRowIsDummy[comp] || r < Height - 1)) {
167               dctArray2 = DCT.forwardDCT(dctArray1);
168               dctArray3 = DCT.quantizeBlock(dctArray2, dct.divisors[qNumber]);
169               // }
170               // else {
171               //   zeroArray[0] = dctArray3[0];
172               //   zeroArray[0] = lastDCvalue[comp];
173               //   dctArray3 = zeroArray;
174               // }
175               huf.HuffmanBlockEncoder(out, dctArray3, lastDCvalue[comp],
176                   dcNumber, acNumber);
177               lastDCvalue[comp] = dctArray3[0];
178             }
179           }
180         }
181       }
182     }
183     huf.flushBuffer(out);
184   }
185
186   private static byte[] eoi = { (byte) 0xFF, (byte) 0xD9 };
187
188   private static byte[] jfif = new byte[] {
189   /* JFIF[0] =*/(byte) 0xff,
190   /* JFIF[1] =*/(byte) 0xe0,
191   /* JFIF[2] =*/0,
192   /* JFIF[3] =*/16,
193   /* JFIF[4] =*/(byte) 0x4a, //'J'
194       /* JFIF[5] =*/(byte) 0x46, //'F'
195       /* JFIF[6] =*/(byte) 0x49, //'I'
196       /* JFIF[7] =*/(byte) 0x46, //'F'
197       /* JFIF[8] =*/0,
198       /* JFIF[9] =*/1,
199       /* JFIF[10] =*/0,
200       /* JFIF[11] =*/0,
201       /* JFIF[12] =*/0,
202       /* JFIF[13] =*/1,
203       /* JFIF[14] =*/0,
204       /* JFIF[15] =*/1,
205       /* JFIF[16] =*/0,
206       /* JFIF[17] =*/0 };
207
208   private static byte[] soi = { (byte) 0xFF, (byte) 0xD8 };
209
210   private String writeHeaders(JpegObj jpegObj, DCT dct) {
211     int i, j, index, offset;
212     int tempArray[];
213
214     // the SOI marker
215     writeMarker(soi);
216
217     // The order of the following headers is quite inconsequential.
218     // the JFIF header
219     writeArray(jfif);
220
221     // Comment Header
222     String comment = null;
223     if (jpegObj.comment != null && jpegObj.comment.length() > 0)
224       writeString(jpegObj.comment, (byte) 0xE1); // App data 1
225     writeString(
226         "JPEG Encoder Copyright 1998, James R. Weeks and BioElectroMech.\n\n",
227         (byte) 0xFE);
228
229     // The DQT header
230     // 0 is the luminance index and 1 is the chrominance index
231     byte dqt[] = new byte[134];
232     dqt[0] = (byte) 0xFF;
233     dqt[1] = (byte) 0xDB;
234     dqt[2] = 0;
235     dqt[3] = (byte) 132;
236     offset = 4;
237     for (i = 0; i < 2; i++) {
238       dqt[offset++] = (byte) ((0 << 4) + i);
239       tempArray = dct.quantum[i];
240       for (j = 0; j < 64; j++) {
241         dqt[offset++] = (byte) tempArray[Huffman.jpegNaturalOrder[j]];
242       }
243     }
244     writeArray(dqt);
245
246     // Start of Frame Header
247     byte sof[] = new byte[19];
248     sof[0] = (byte) 0xFF;
249     sof[1] = (byte) 0xC0;
250     sof[2] = 0;
251     sof[3] = 17;
252     sof[4] = (byte) jpegObj.precision;
253     sof[5] = (byte) ((jpegObj.imageHeight >> 8) & 0xFF);
254     sof[6] = (byte) ((jpegObj.imageHeight) & 0xFF);
255     sof[7] = (byte) ((jpegObj.imageWidth >> 8) & 0xFF);
256     sof[8] = (byte) ((jpegObj.imageWidth) & 0xFF);
257     sof[9] = (byte) jpegObj.numberOfComponents;
258     index = 10;
259     for (i = 0; i < sof[9]; i++) {
260       sof[index++] = (byte) jpegObj.compID[i];
261       sof[index++] = (byte) ((jpegObj.hsampFactor[i] << 4) + jpegObj.vsampFactor[i]);
262       sof[index++] = (byte) jpegObj.qtableNumber[i];
263     }
264     writeArray(sof);
265
266     WriteDHTHeader(Huffman.bitsDCluminance, Huffman.valDCluminance);
267     WriteDHTHeader(Huffman.bitsACluminance, Huffman.valACluminance);
268     WriteDHTHeader(Huffman.bitsDCchrominance, Huffman.valDCchrominance);
269     WriteDHTHeader(Huffman.bitsACchrominance, Huffman.valACchrominance);
270
271     // Start of Scan Header
272     byte sos[] = new byte[14];
273     sos[0] = (byte) 0xFF;
274     sos[1] = (byte) 0xDA;
275     sos[2] = 0;
276     sos[3] = 12;
277     sos[4] = (byte) jpegObj.numberOfComponents;
278     index = 5;
279     for (i = 0; i < sos[4]; i++) {
280       sos[index++] = (byte) jpegObj.compID[i];
281       sos[index++] = (byte) ((jpegObj.dctableNumber[i] << 4) + jpegObj.actableNumber[i]);
282     }
283     sos[index++] = (byte) jpegObj.ss;
284     sos[index++] = (byte) jpegObj.se;
285     sos[index++] = (byte) ((jpegObj.ah << 4) + jpegObj.al);
286     writeArray(sos);
287     return comment;
288   }
289
290   private void writeString(String s, byte id) {
291     int len = s.length();
292     int i0 = 0;
293     String suffix = applicationTag;
294     while (i0 < len) {
295       int nBytes = len - i0;
296       if (nBytes > CONTINUE_MAX_BUFFER) {
297         nBytes = CONTINUE_MAX;
298         // but break only at line breaks
299         int pt = s.lastIndexOf('\n', i0 + nBytes);
300         if (pt > i0 + 1)
301           nBytes = pt - i0;
302       }
303       if (i0 + nBytes == len)
304         suffix = "";
305       writeTag(nBytes + suffix.length(), id);
306       writeArray(s.substring(i0, i0 + nBytes).getBytes());
307       if (suffix.length() > 0)
308         writeArray(suffix.getBytes());
309       i0 += nBytes;
310     }
311   }
312
313   private void writeTag(int length, byte id) {
314     length += 2;
315     byte com[] = new byte[4];
316     com[0] = (byte) 0xFF;
317     com[1] = id;
318     com[2] = (byte) ((length >> 8) & 0xFF);
319     com[3] = (byte) (length & 0xFF);
320     writeArray(com);
321   }
322
323   void WriteDHTHeader(int[] bits, int[] val) {
324     // hansonr@stolaf.edu: simplified code.
325     byte[] dht;
326     int bytes = 0;
327     for (int j = 1; j < 17; j++)
328       bytes += bits[j];
329     dht = new byte[21 + bytes];
330     dht[0] = (byte) 0xFF;
331     dht[1] = (byte) 0xC4;
332     int index = 4;
333     for (int j = 0; j < 17; j++)
334       dht[index++] = (byte) bits[j];
335     for (int j = 0; j < bytes; j++)
336       dht[index++] = (byte) val[j];
337     dht[2] = (byte) (((index - 2) >> 8) & 0xFF);
338     dht[3] = (byte) ((index - 2) & 0xFF);
339     writeArray(dht);
340   }
341
342   void writeMarker(byte[] data) {
343     out.write(data, 0, 2);
344   }
345
346   void writeArray(byte[] data) {
347     out.write(data, 0, data.length);
348   }
349
350 }
351
352 // This class incorporates quality scaling as implemented in the JPEG-6a
353 // library.
354
355 /*
356  * DCT - A Java implementation of the Discreet Cosine Transform
357  */
358
359 class DCT {
360
361   /**
362    * DCT Block Size - default 8
363    */
364   private final static int N = 8;
365   private final static int NN = N * N;
366
367   /**
368    * Image Quality (0-100) - default 80 (good image / good compression)
369    */
370   //public int QUALITY = 80;
371
372   int[][] quantum = AU.newInt2(2);
373   double[][] divisors = AU.newDouble2(2);
374
375   /**
376    * Quantitization Matrix for luminace.
377    */
378   private int quantum_luminance[] = new int[NN];
379   private double DivisorsLuminance[] = new double[NN];
380
381   /**
382    * Quantitization Matrix for chrominance.
383    */
384   private int quantum_chrominance[] = new int[NN];
385   private double DivisorsChrominance[] = new double[NN];
386
387   /**
388    * Constructs a new DCT object. Initializes the cosine transform matrix these
389    * are used when computing the DCT and it's inverse. This also initializes the
390    * run length counters and the ZigZag sequence. Note that the image quality
391    * can be worse than 25 however the image will be extemely pixelated, usually
392    * to a block size of N.
393    * 
394    * @param quality
395    *        The quality of the image (0 worst - 100 best)
396    * 
397    */
398   DCT(int quality) {
399     initMatrix(quality);
400   }
401
402   /*
403    * This method sets up the quantization matrix for luminance and
404    * chrominance using the Quality parameter.
405    */
406   private void initMatrix(int quality) {
407     // converting quality setting to that specified in the jpeg_quality_scaling
408     // method in the IJG Jpeg-6a C libraries
409
410     quality = (quality < 1 ? 1 : quality > 100 ? 100 : quality);
411     quality = (quality < 50 ? 5000 / quality : 200 - quality * 2);
412
413     // Creating the luminance matrix
414
415     quantum_luminance[0] = 16;
416     quantum_luminance[1] = 11;
417     quantum_luminance[2] = 10;
418     quantum_luminance[3] = 16;
419     quantum_luminance[4] = 24;
420     quantum_luminance[5] = 40;
421     quantum_luminance[6] = 51;
422     quantum_luminance[7] = 61;
423     quantum_luminance[8] = 12;
424     quantum_luminance[9] = 12;
425     quantum_luminance[10] = 14;
426     quantum_luminance[11] = 19;
427     quantum_luminance[12] = 26;
428     quantum_luminance[13] = 58;
429     quantum_luminance[14] = 60;
430     quantum_luminance[15] = 55;
431     quantum_luminance[16] = 14;
432     quantum_luminance[17] = 13;
433     quantum_luminance[18] = 16;
434     quantum_luminance[19] = 24;
435     quantum_luminance[20] = 40;
436     quantum_luminance[21] = 57;
437     quantum_luminance[22] = 69;
438     quantum_luminance[23] = 56;
439     quantum_luminance[24] = 14;
440     quantum_luminance[25] = 17;
441     quantum_luminance[26] = 22;
442     quantum_luminance[27] = 29;
443     quantum_luminance[28] = 51;
444     quantum_luminance[29] = 87;
445     quantum_luminance[30] = 80;
446     quantum_luminance[31] = 62;
447     quantum_luminance[32] = 18;
448     quantum_luminance[33] = 22;
449     quantum_luminance[34] = 37;
450     quantum_luminance[35] = 56;
451     quantum_luminance[36] = 68;
452     quantum_luminance[37] = 109;
453     quantum_luminance[38] = 103;
454     quantum_luminance[39] = 77;
455     quantum_luminance[40] = 24;
456     quantum_luminance[41] = 35;
457     quantum_luminance[42] = 55;
458     quantum_luminance[43] = 64;
459     quantum_luminance[44] = 81;
460     quantum_luminance[45] = 104;
461     quantum_luminance[46] = 113;
462     quantum_luminance[47] = 92;
463     quantum_luminance[48] = 49;
464     quantum_luminance[49] = 64;
465     quantum_luminance[50] = 78;
466     quantum_luminance[51] = 87;
467     quantum_luminance[52] = 103;
468     quantum_luminance[53] = 121;
469     quantum_luminance[54] = 120;
470     quantum_luminance[55] = 101;
471     quantum_luminance[56] = 72;
472     quantum_luminance[57] = 92;
473     quantum_luminance[58] = 95;
474     quantum_luminance[59] = 98;
475     quantum_luminance[60] = 112;
476     quantum_luminance[61] = 100;
477     quantum_luminance[62] = 103;
478     quantum_luminance[63] = 99;
479
480     AANscale(DivisorsLuminance, quantum_luminance, quality);
481
482     // Creating the chrominance matrix
483
484     for (int i = 4; i < 64; i++)
485       quantum_chrominance[i] = 99;
486
487     quantum_chrominance[0] = 17;
488     quantum_chrominance[1] = 18;
489     quantum_chrominance[2] = 24;
490     quantum_chrominance[3] = 47;
491
492     quantum_chrominance[8] = 18;
493     quantum_chrominance[9] = 21;
494     quantum_chrominance[10] = 26;
495     quantum_chrominance[11] = 66;
496
497     quantum_chrominance[16] = 24;
498     quantum_chrominance[17] = 26;
499     quantum_chrominance[18] = 56;
500
501     quantum_chrominance[24] = 47;
502     quantum_chrominance[25] = 66;
503
504     AANscale(DivisorsChrominance, quantum_chrominance, quality);
505
506     // quantum and Divisors are objects used to hold the appropriate matices
507
508     quantum[0] = quantum_luminance;
509     quantum[1] = quantum_chrominance;
510
511     divisors[0] = DivisorsLuminance;
512     divisors[1] = DivisorsChrominance;
513
514   }
515
516   private final static double[] AANscaleFactor = { 1.0, 1.387039845,
517       1.306562965, 1.175875602, 1.0, 0.785694958, 0.541196100, 0.275899379 };
518
519   static private void AANscale(double[] divisors, int[] values, int quality) {
520
521     for (int j = 0; j < 64; j++) {
522       int temp = (values[j] * quality + 50) / 100;
523       values[j] = (temp < 1 ? 1 : temp > 255 ? 255 : temp);
524     }
525
526     for (int i = 0, index = 0; i < 8; i++)
527       for (int j = 0; j < 8; j++, index++)
528         // The divisors for the LL&M method (the slow integer method used in
529         // jpeg 6a library).  This method is currently (04/04/98) incompletely
530         // implemented.
531         // DivisorsLuminance[index] = ((double) quantum_luminance[index]) << 3;
532         // The divisors for the AAN method (the float method used in jpeg 6a library.
533         divisors[index] = (0.125 / (values[index] * AANscaleFactor[i] * AANscaleFactor[j]));
534   }
535
536   /*
537    * This method preforms forward DCT on a block of image data using
538    * the literal method specified for a 2-D Discrete Cosine Transform.
539    * It is included as a curiosity and can give you an idea of the
540    * difference in the compression result (the resulting image quality)
541    * by comparing its output to the output of the AAN method below.
542    * It is ridiculously inefficient.
543    */
544
545   // For now the final output is unusable.  The associated quantization step
546   // needs some tweaking.  If you get this part working, please let me know.
547   /*
548     public double[][] forwardDCTExtreme(float input[][])
549     {
550       double output[][] = new double[N][N];
551       int v, u, x, y;
552       for (v = 0; v < 8; v++) {
553         for (u = 0; u < 8; u++) {
554           for (x = 0; x < 8; x++) {
555             for (y = 0; y < 8; y++) {
556               output[v][u] += input[x][y] * 
557                   Math.cos(((double)(2*x + 1)*(double)u*Math.PI)/16)*
558                   Math.cos(((double)(2*y + 1)*(double)v*Math.PI)/16);
559             }
560           }
561           output[v][u] *= (0.25)*((u == 0) ? (1.0/Math.sqrt(2)) : (double) 1.0)*((v == 0) ? (1.0/Math.sqrt(2)) : (double) 1.0);
562         }
563       }
564       return output;
565     }
566     
567   */
568   /*
569    * This method preforms a DCT on a block of image data using the AAN
570    * method as implemented in the IJG Jpeg-6a library.
571    */
572   static double[][] forwardDCT(float input[][]) {
573     double output[][] = new double[N][N];
574     double tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
575     double tmp10, tmp11, tmp12, tmp13;
576     double z1, z2, z3, z4, z5, z11, z13;
577     // Subtracts 128 from the input values
578     for (int i = 0; i < 8; i++)
579       for (int j = 0; j < 8; j++)
580         output[i][j] = (input[i][j] - 128.0);
581     // input[i][j] -= 128;
582
583     for (int i = 0; i < 8; i++) {
584       tmp0 = output[i][0] + output[i][7];
585       tmp7 = output[i][0] - output[i][7];
586       tmp1 = output[i][1] + output[i][6];
587       tmp6 = output[i][1] - output[i][6];
588       tmp2 = output[i][2] + output[i][5];
589       tmp5 = output[i][2] - output[i][5];
590       tmp3 = output[i][3] + output[i][4];
591       tmp4 = output[i][3] - output[i][4];
592
593       tmp10 = tmp0 + tmp3;
594       tmp13 = tmp0 - tmp3;
595       tmp11 = tmp1 + tmp2;
596       tmp12 = tmp1 - tmp2;
597
598       output[i][0] = tmp10 + tmp11;
599       output[i][4] = tmp10 - tmp11;
600
601       z1 = (tmp12 + tmp13) * 0.707106781;
602       output[i][2] = tmp13 + z1;
603       output[i][6] = tmp13 - z1;
604
605       tmp10 = tmp4 + tmp5;
606       tmp11 = tmp5 + tmp6;
607       tmp12 = tmp6 + tmp7;
608
609       z5 = (tmp10 - tmp12) * 0.382683433;
610       z2 = 0.541196100 * tmp10 + z5;
611       z4 = 1.306562965 * tmp12 + z5;
612       z3 = tmp11 * 0.707106781;
613
614       z11 = tmp7 + z3;
615       z13 = tmp7 - z3;
616
617       output[i][5] = z13 + z2;
618       output[i][3] = z13 - z2;
619       output[i][1] = z11 + z4;
620       output[i][7] = z11 - z4;
621     }
622
623     for (int i = 0; i < 8; i++) {
624       tmp0 = output[0][i] + output[7][i];
625       tmp7 = output[0][i] - output[7][i];
626       tmp1 = output[1][i] + output[6][i];
627       tmp6 = output[1][i] - output[6][i];
628       tmp2 = output[2][i] + output[5][i];
629       tmp5 = output[2][i] - output[5][i];
630       tmp3 = output[3][i] + output[4][i];
631       tmp4 = output[3][i] - output[4][i];
632
633       tmp10 = tmp0 + tmp3;
634       tmp13 = tmp0 - tmp3;
635       tmp11 = tmp1 + tmp2;
636       tmp12 = tmp1 - tmp2;
637
638       output[0][i] = tmp10 + tmp11;
639       output[4][i] = tmp10 - tmp11;
640
641       z1 = (tmp12 + tmp13) * 0.707106781;
642       output[2][i] = tmp13 + z1;
643       output[6][i] = tmp13 - z1;
644
645       tmp10 = tmp4 + tmp5;
646       tmp11 = tmp5 + tmp6;
647       tmp12 = tmp6 + tmp7;
648
649       z5 = (tmp10 - tmp12) * 0.382683433;
650       z2 = 0.541196100 * tmp10 + z5;
651       z4 = 1.306562965 * tmp12 + z5;
652       z3 = tmp11 * 0.707106781;
653
654       z11 = tmp7 + z3;
655       z13 = tmp7 - z3;
656
657       output[5][i] = z13 + z2;
658       output[3][i] = z13 - z2;
659       output[1][i] = z11 + z4;
660       output[7][i] = z11 - z4;
661     }
662
663     return output;
664   }
665
666   /*
667    * This method quantitizes data and rounds it to the nearest integer.
668    */
669   static int[] quantizeBlock(double inputData[][], double[] divisorsCode) {
670     int outputData[] = new int[NN];
671     for (int i = 0, index = 0; i < 8; i++)
672       for (int j = 0; j < 8; j++, index++)
673         // The second line results in significantly better compression.
674         outputData[index] = (int) (Math.round(inputData[i][j]
675             * divisorsCode[index]));
676     //                        outputData[index] = (int)(((inputData[i][j] * (((double[]) (Divisors[code]))[index])) + 16384.5) -16384);
677     return outputData;
678   }
679
680   /*
681    * This is the method for quantizing a block DCT'ed with forwardDCTExtreme
682    * This method quantitizes data and rounds it to the nearest integer.
683    */
684
685   /*
686
687     public double[][] forwardDCTExtreme(float input[][])
688     {
689       double output[][] = new double[N][N];
690       int v, u, x, y;
691       for (v = 0; v < 8; v++) {
692         for (u = 0; u < 8; u++) {
693           for (x = 0; x < 8; x++) {
694             for (y = 0; y < 8; y++) {
695               output[v][u] += input[x][y] * 
696                   Math.cos(((double)(2*x + 1)*(double)u*Math.PI)/16)*
697                   Math.cos(((double)(2*y + 1)*(double)v*Math.PI)/16);
698             }
699           }
700           output[v][u] *= (0.25)*((u == 0) ? (1.0/Math.sqrt(2)) : (double) 1.0)*((v == 0) ? (1.0/Math.sqrt(2)) : (double) 1.0);
701         }
702       }
703       return output;
704     }
705
706    */
707   /*
708     public int[] quantizeBlockExtreme(double inputData[][], int code)
709     {
710       int outputData[] = new int[NN];
711       int i, j;
712       int index;
713       index = 0;
714       for (i = 0; i < 8; i++) {
715         for (j = 0; j < 8; j++) {
716           outputData[index] = (int)(Math.round(inputData[i][j] / (((int[]) (quantum[code]))[index])));
717           index++;
718         }
719       }
720       
721       return outputData;
722     }
723   */
724 }
725
726 // This class was modified by James R. Weeks on 3/27/98.
727 // It now incorporates Huffman table derivation as in the C jpeg library
728 // from the IJG, Jpeg-6a.
729
730 class Huffman {
731   private int bufferPutBits, bufferPutBuffer;
732   int imageHeight;
733   int imageWidth;
734   private int dc_matrix0[][];
735   private int ac_matrix0[][];
736   private int dc_matrix1[][];
737   private int ac_matrix1[][];
738   private int[][][] dc_matrix;
739   private int[][][] ac_matrix;
740   //private int code;
741   int numOfDCTables;
742   int numOfACTables;
743   final static int[] bitsDCluminance = { 0x00, 0, 1, 5, 1, 1, 1, 1, 1, 1, 0, 0,
744       0, 0, 0, 0, 0 };
745   final static int[] valDCluminance = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 };
746   final static int[] bitsDCchrominance = { 0x01, 0, 3, 1, 1, 1, 1, 1, 1, 1, 1,
747       1, 0, 0, 0, 0, 0 };
748   final static int[] valDCchrominance = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 };
749   final static int[] bitsACluminance = { 0x10, 0, 2, 1, 3, 3, 2, 4, 3, 5, 5, 4,
750       4, 0, 0, 1, 0x7d };
751   final static int[] valACluminance = { 0x01, 0x02, 0x03, 0x00, 0x04, 0x11,
752       0x05, 0x12, 0x21, 0x31, 0x41, 0x06, 0x13, 0x51, 0x61, 0x07, 0x22, 0x71,
753       0x14, 0x32, 0x81, 0x91, 0xa1, 0x08, 0x23, 0x42, 0xb1, 0xc1, 0x15, 0x52,
754       0xd1, 0xf0, 0x24, 0x33, 0x62, 0x72, 0x82, 0x09, 0x0a, 0x16, 0x17, 0x18,
755       0x19, 0x1a, 0x25, 0x26, 0x27, 0x28, 0x29, 0x2a, 0x34, 0x35, 0x36, 0x37,
756       0x38, 0x39, 0x3a, 0x43, 0x44, 0x45, 0x46, 0x47, 0x48, 0x49, 0x4a, 0x53,
757       0x54, 0x55, 0x56, 0x57, 0x58, 0x59, 0x5a, 0x63, 0x64, 0x65, 0x66, 0x67,
758       0x68, 0x69, 0x6a, 0x73, 0x74, 0x75, 0x76, 0x77, 0x78, 0x79, 0x7a, 0x83,
759       0x84, 0x85, 0x86, 0x87, 0x88, 0x89, 0x8a, 0x92, 0x93, 0x94, 0x95, 0x96,
760       0x97, 0x98, 0x99, 0x9a, 0xa2, 0xa3, 0xa4, 0xa5, 0xa6, 0xa7, 0xa8, 0xa9,
761       0xaa, 0xb2, 0xb3, 0xb4, 0xb5, 0xb6, 0xb7, 0xb8, 0xb9, 0xba, 0xc2, 0xc3,
762       0xc4, 0xc5, 0xc6, 0xc7, 0xc8, 0xc9, 0xca, 0xd2, 0xd3, 0xd4, 0xd5, 0xd6,
763       0xd7, 0xd8, 0xd9, 0xda, 0xe1, 0xe2, 0xe3, 0xe4, 0xe5, 0xe6, 0xe7, 0xe8,
764       0xe9, 0xea, 0xf1, 0xf2, 0xf3, 0xf4, 0xf5, 0xf6, 0xf7, 0xf8, 0xf9, 0xfa };
765   final static int[] bitsACchrominance = { 0x11, 0, 2, 1, 2, 4, 4, 3, 4, 7, 5,
766       4, 4, 0, 1, 2, 0x77 };
767   final static int[] valACchrominance = { 0x00, 0x01, 0x02, 0x03, 0x11, 0x04,
768       0x05, 0x21, 0x31, 0x06, 0x12, 0x41, 0x51, 0x07, 0x61, 0x71, 0x13, 0x22,
769       0x32, 0x81, 0x08, 0x14, 0x42, 0x91, 0xa1, 0xb1, 0xc1, 0x09, 0x23, 0x33,
770       0x52, 0xf0, 0x15, 0x62, 0x72, 0xd1, 0x0a, 0x16, 0x24, 0x34, 0xe1, 0x25,
771       0xf1, 0x17, 0x18, 0x19, 0x1a, 0x26, 0x27, 0x28, 0x29, 0x2a, 0x35, 0x36,
772       0x37, 0x38, 0x39, 0x3a, 0x43, 0x44, 0x45, 0x46, 0x47, 0x48, 0x49, 0x4a,
773       0x53, 0x54, 0x55, 0x56, 0x57, 0x58, 0x59, 0x5a, 0x63, 0x64, 0x65, 0x66,
774       0x67, 0x68, 0x69, 0x6a, 0x73, 0x74, 0x75, 0x76, 0x77, 0x78, 0x79, 0x7a,
775       0x82, 0x83, 0x84, 0x85, 0x86, 0x87, 0x88, 0x89, 0x8a, 0x92, 0x93, 0x94,
776       0x95, 0x96, 0x97, 0x98, 0x99, 0x9a, 0xa2, 0xa3, 0xa4, 0xa5, 0xa6, 0xa7,
777       0xa8, 0xa9, 0xaa, 0xb2, 0xb3, 0xb4, 0xb5, 0xb6, 0xb7, 0xb8, 0xb9, 0xba,
778       0xc2, 0xc3, 0xc4, 0xc5, 0xc6, 0xc7, 0xc8, 0xc9, 0xca, 0xd2, 0xd3, 0xd4,
779       0xd5, 0xd6, 0xd7, 0xd8, 0xd9, 0xda, 0xe2, 0xe3, 0xe4, 0xe5, 0xe6, 0xe7,
780       0xe8, 0xe9, 0xea, 0xf2, 0xf3, 0xf4, 0xf5, 0xf6, 0xf7, 0xf8, 0xf9, 0xfa };
781
782   /*
783    * jpegNaturalOrder[i] is the natural-order position of the i'th element
784    * of zigzag order.
785    */
786   final static int[] jpegNaturalOrder = { 0, 1, 8, 16, 9, 2, 3, 10, 17, 24, 32,
787       25, 18, 11, 4, 5, 12, 19, 26, 33, 40, 48, 41, 34, 27, 20, 13, 6, 7, 14,
788       21, 28, 35, 42, 49, 56, 57, 50, 43, 36, 29, 22, 15, 23, 30, 37, 44, 51,
789       58, 59, 52, 45, 38, 31, 39, 46, 53, 60, 61, 54, 47, 55, 62, 63, };
790
791   Huffman(int width, int height) {
792     initHuf();
793     imageWidth = width;
794     imageHeight = height;
795
796   }
797
798   /**
799    * HuffmanBlockEncoder run length encodes and Huffman encodes the quantized
800    * data.
801    * 
802    * @param out
803    * @param zigzag
804    * @param prec
805    * @param dcCode
806    * @param acCode
807    **/
808
809   void HuffmanBlockEncoder(OC out, int zigzag[], int prec,
810                            int dcCode, int acCode) {
811     int temp, temp2, nbits, k, r, i;
812
813     numOfDCTables = 2;
814     numOfACTables = 2;
815
816     int[][] matrixDC = dc_matrix[dcCode];
817     int[][] matrixAC = ac_matrix[acCode];
818
819     // The DC portion
820
821     temp = temp2 = zigzag[0] - prec;
822     if (temp < 0) {
823       temp = -temp;
824       temp2--;
825     }
826     nbits = 0;
827     while (temp != 0) {
828       nbits++;
829       temp >>= 1;
830     }
831     //        if (nbits > 11) nbits = 11;
832     bufferIt(out, matrixDC[nbits][0], matrixDC[nbits][1]);
833     // The arguments in bufferIt are code and size.
834     if (nbits != 0) {
835       bufferIt(out, temp2, nbits);
836     }
837
838     // The AC portion
839
840     r = 0;
841
842     for (k = 1; k < 64; k++) {
843       if ((temp = zigzag[jpegNaturalOrder[k]]) == 0) {
844         r++;
845       } else {
846         while (r > 15) {
847           bufferIt(out, matrixAC[0xF0][0], matrixAC[0xF0][1]);
848           r -= 16;
849         }
850         temp2 = temp;
851         if (temp < 0) {
852           temp = -temp;
853           temp2--;
854         }
855         nbits = 1;
856         while ((temp >>= 1) != 0) {
857           nbits++;
858         }
859         i = (r << 4) + nbits;
860         bufferIt(out, matrixAC[i][0], matrixAC[i][1]);
861         bufferIt(out, temp2, nbits);
862
863         r = 0;
864       }
865     }
866
867     if (r > 0) {
868       bufferIt(out, matrixAC[0][0], matrixAC[0][1]);
869     }
870
871   }
872
873   // Uses an integer long (32 bits) buffer to store the Huffman encoded bits
874   // and sends them to out by the byte.
875
876   void bufferIt(OC out, int code, int size) {
877     int putBuffer = code;
878     int putBits = bufferPutBits;
879
880     putBuffer &= (1 << size) - 1;
881     putBits += size;
882     putBuffer <<= 24 - putBits;
883     putBuffer |= bufferPutBuffer;
884
885     while (putBits >= 8) {
886       int c = ((putBuffer >> 16) & 0xFF);
887       out.writeByteAsInt(c);
888       if (c == 0xFF) {
889         out.writeByteAsInt(0);
890       }
891       putBuffer <<= 8;
892       putBits -= 8;
893     }
894     bufferPutBuffer = putBuffer;
895     bufferPutBits = putBits;
896
897   }
898
899   void flushBuffer(OC out) {
900     int putBuffer = bufferPutBuffer;
901     int putBits = bufferPutBits;
902     while (putBits >= 8) {
903       int c = ((putBuffer >> 16) & 0xFF);
904       out.writeByteAsInt(c);
905       if (c == 0xFF) {
906         out.writeByteAsInt(0);
907       }
908       putBuffer <<= 8;
909       putBits -= 8;
910     }
911     if (putBits > 0) {
912       int c = ((putBuffer >> 16) & 0xFF);
913       out.writeByteAsInt(c);
914     }
915   }
916
917   /*
918    * Initialisation of the Huffman codes for Luminance and Chrominance.
919    * This code results in the same tables created in the IJG Jpeg-6a
920    * library.
921    */
922
923   private void initHuf() {
924     dc_matrix0 = new int[12][2];
925     dc_matrix1 = new int[12][2];
926     ac_matrix0 = new int[255][2];
927     ac_matrix1 = new int[255][2];
928     dc_matrix = AU.newInt3(2, -1);
929     ac_matrix = AU.newInt3(2, -1);
930     int p, l, i, lastp, si, code;
931     int[] huffsize = new int[257];
932     int[] huffcode = new int[257];
933
934     /*
935      * init of the DC values for the chrominance
936      * [][0] is the code   [][1] is the number of bit
937      */
938
939     p = 0;
940     for (l = 1; l <= 16; l++) {
941       //      for (i = 1; i <= bitsDCchrominance[l]; i++)
942       for (i = bitsDCchrominance[l]; --i >= 0;) {
943         huffsize[p++] = l; //that's an "el", not a "one"
944       }
945     }
946     huffsize[p] = 0;
947     lastp = p;
948
949     code = 0;
950     si = huffsize[0];
951     p = 0;
952     while (huffsize[p] != 0) {
953       while (huffsize[p] == si) {
954         huffcode[p++] = code;
955         code++;
956       }
957       code <<= 1;
958       si++;
959     }
960
961     for (p = 0; p < lastp; p++) {
962       dc_matrix1[valDCchrominance[p]][0] = huffcode[p];
963       dc_matrix1[valDCchrominance[p]][1] = huffsize[p];
964     }
965
966     /*
967      * Init of the AC huffman code for the chrominance
968      * matrix [][][0] is the code & matrix[][][1] is the number of bit needed
969      */
970
971     p = 0;
972     for (l = 1; l <= 16; l++) {
973       for (i = bitsACchrominance[l]; --i >= 0;)
974       //      for (i = 1; i <= bitsACchrominance[l]; i++)
975       {
976         huffsize[p++] = l;
977       }
978     }
979     huffsize[p] = 0;
980     lastp = p;
981
982     code = 0;
983     si = huffsize[0];
984     p = 0;
985     while (huffsize[p] != 0) {
986       while (huffsize[p] == si) {
987         huffcode[p++] = code;
988         code++;
989       }
990       code <<= 1;
991       si++;
992     }
993
994     for (p = 0; p < lastp; p++) {
995       ac_matrix1[valACchrominance[p]][0] = huffcode[p];
996       ac_matrix1[valACchrominance[p]][1] = huffsize[p];
997     }
998
999     /*
1000      * init of the DC values for the luminance
1001      * [][0] is the code   [][1] is the number of bit
1002      */
1003     p = 0;
1004     for (l = 1; l <= 16; l++) {
1005       //      for (i = 1; i <= bitsDCluminance[l]; i++)
1006       for (i = bitsDCluminance[l]; --i >= 0;) {
1007         huffsize[p++] = l;
1008       }
1009     }
1010     huffsize[p] = 0;
1011     lastp = p;
1012
1013     code = 0;
1014     si = huffsize[0];
1015     p = 0;
1016     while (huffsize[p] != 0) {
1017       while (huffsize[p] == si) {
1018         huffcode[p++] = code;
1019         code++;
1020       }
1021       code <<= 1;
1022       si++;
1023     }
1024
1025     for (p = 0; p < lastp; p++) {
1026       dc_matrix0[valDCluminance[p]][0] = huffcode[p];
1027       dc_matrix0[valDCluminance[p]][1] = huffsize[p];
1028     }
1029
1030     /*
1031      * Init of the AC huffman code for luminance
1032      * matrix [][][0] is the code & matrix[][][1] is the number of bit
1033      */
1034
1035     p = 0;
1036     for (l = 1; l <= 16; l++) {
1037       //      for (i = 1; i <= bitsACluminance[l]; i++)
1038       for (i = bitsACluminance[l]; --i >= 0;) {
1039         huffsize[p++] = l;
1040       }
1041     }
1042     huffsize[p] = 0;
1043     lastp = p;
1044
1045     code = 0;
1046     si = huffsize[0];
1047     p = 0;
1048     while (huffsize[p] != 0) {
1049       while (huffsize[p] == si) {
1050         huffcode[p++] = code;
1051         code++;
1052       }
1053       code <<= 1;
1054       si++;
1055     }
1056     for (int q = 0; q < lastp; q++) {
1057       ac_matrix0[valACluminance[q]][0] = huffcode[q];
1058       ac_matrix0[valACluminance[q]][1] = huffsize[q];
1059     }
1060
1061     dc_matrix[0] = dc_matrix0;
1062     dc_matrix[1] = dc_matrix1;
1063     ac_matrix[0] = ac_matrix0;
1064     ac_matrix[1] = ac_matrix1;
1065   }
1066
1067 }
1068
1069 /*
1070  * JpegInfo - Given an image, sets default information about it and divides
1071  * it into its constituant components, downsizing those that need to be.
1072  */
1073
1074 class JpegObj {
1075   String comment;
1076   int imageHeight;
1077   int imageWidth;
1078   int blockWidth[];
1079   int blockHeight[];
1080
1081   int precision = 8;
1082   int numberOfComponents = 3;
1083   float[][][] components;
1084   int[] compID = { 1, 2, 3 };
1085   int[] hsampFactor = { 1, 1, 1 };
1086   int[] vsampFactor = { 1, 1, 1 };
1087   int[] qtableNumber = { 0, 1, 1 };
1088   int[] dctableNumber = { 0, 1, 1 };
1089   int[] actableNumber = { 0, 1, 1 };
1090   private boolean[] lastColumnIsDummy = { false, false, false };
1091   private boolean[] lastRowIsDummy = { false, false, false };
1092   int ss = 0;
1093   int se = 63;
1094   int ah = 0;
1095   int al = 0;
1096   private int compWidth[];
1097   private int compHeight[];
1098   private int maxHsampFactor;
1099   private int maxVsampFactor;
1100
1101   public JpegObj() {
1102     components = AU.newFloat3(numberOfComponents, -1);
1103     compWidth = new int[numberOfComponents];
1104     compHeight = new int[numberOfComponents];
1105     blockWidth = new int[numberOfComponents];
1106     blockHeight = new int[numberOfComponents];
1107   }
1108
1109   /*
1110    * This method creates and fills three arrays, Y, Cb, and Cr using the
1111    * input image.
1112    */
1113
1114   void getYCCArray(int[] pixels) {
1115     // In order to minimize the chance that grabPixels will throw an exception
1116     // it may be necessary to grab some pixels every few scanlines and process
1117     // those before going for more.  The time expense may be prohibitive.
1118     // However, for a situation where memory overhead is a concern, this may be
1119     // the only choice.
1120     maxHsampFactor = 1;
1121     maxVsampFactor = 1;
1122     for (int y = 0; y < numberOfComponents; y++) {
1123       maxHsampFactor = Math.max(maxHsampFactor, hsampFactor[y]);
1124       maxVsampFactor = Math.max(maxVsampFactor, vsampFactor[y]);
1125     }
1126     for (int y = 0; y < numberOfComponents; y++) {
1127       compWidth[y] = (((imageWidth % 8 != 0) ? ((int) Math
1128           .ceil(imageWidth / 8.0)) * 8 : imageWidth) / maxHsampFactor)
1129           * hsampFactor[y];
1130       if (compWidth[y] != ((imageWidth / maxHsampFactor) * hsampFactor[y])) {
1131         lastColumnIsDummy[y] = true;
1132       }
1133       // results in a multiple of 8 for compWidth
1134       // this will make the rest of the program fail for the unlikely
1135       // event that someone tries to compress an 16 x 16 pixel image
1136       // which would of course be worse than pointless
1137       blockWidth[y] = (int) Math.ceil(compWidth[y] / 8.0);
1138       compHeight[y] = (((imageHeight % 8 != 0) ? ((int) Math
1139           .ceil(imageHeight / 8.0)) * 8 : imageHeight) / maxVsampFactor)
1140           * vsampFactor[y];
1141       if (compHeight[y] != ((imageHeight / maxVsampFactor) * vsampFactor[y])) {
1142         lastRowIsDummy[y] = true;
1143       }
1144       blockHeight[y] = (int) Math.ceil(compHeight[y] / 8.0);
1145     }
1146     float Y[][] = new float[compHeight[0]][compWidth[0]];
1147     float Cr1[][] = new float[compHeight[0]][compWidth[0]];
1148     float Cb1[][] = new float[compHeight[0]][compWidth[0]];
1149     //float Cb2[][] = new float[compHeight[1]][compWidth[1]];
1150     //float Cr2[][] = new float[compHeight[2]][compWidth[2]];
1151     for (int pt = 0, y = 0; y < imageHeight; ++y) {
1152       for (int x = 0; x < imageWidth; ++x, pt++) {
1153         int p = pixels[pt];
1154         int r = ((p >> 16) & 0xff);
1155         int g = ((p >> 8) & 0xff);
1156         int b = (p & 0xff);
1157         // The following three lines are a more correct color conversion but
1158         // the current conversion technique is sufficient and results in a higher
1159         // compression rate.
1160         // Y[y][x] = 16 + (float)(0.8588*(0.299 * (float)r + 0.587 * (float)g + 0.114 * (float)b ));
1161         // Cb1[y][x] = 128 + (float)(0.8784*(-0.16874 * (float)r - 0.33126 * (float)g + 0.5 * (float)b));
1162         // Cr1[y][x] = 128 + (float)(0.8784*(0.5 * (float)r - 0.41869 * (float)g - 0.08131 * (float)b));
1163         Y[y][x] = (float) ((0.299 * r + 0.587 * g + 0.114 * b));
1164         Cb1[y][x] = 128 + (float) ((-0.16874 * r - 0.33126 * g + 0.5 * b));
1165         Cr1[y][x] = 128 + (float) ((0.5 * r - 0.41869 * g - 0.08131 * b));
1166       }
1167     }
1168
1169     // Need a way to set the H and V sample factors before allowing downsampling.
1170     // For now (04/04/98) downsampling must be hard coded.
1171     // Until a better downsampler is implemented, this will not be done.
1172     // Downsampling is currently supported.  The downsampling method here
1173     // is a simple box filter.
1174
1175     components[0] = Y;
1176     //        Cb2 = DownSample(Cb1, 1);
1177     components[1] = Cb1;
1178     //        Cr2 = DownSample(Cr1, 2);
1179     components[2] = Cr1;
1180   }
1181   /*  
1182     float[][] DownSample(float[][] C, int comp)
1183     {
1184       int inrow, incol;
1185       int outrow, outcol;
1186       float output[][];
1187       int bias;
1188       inrow = 0;
1189       incol = 0;
1190       int cHeight = compHeight[comp];
1191       int cWidth = compWidth[comp];
1192       output = new float[cHeight][cWidth];
1193       
1194       for (outrow = 0; outrow < cHeight; outrow++) {
1195         bias = 1;
1196         for (outcol = 0; outcol < cWidth; outcol++) {
1197           output[outrow][outcol] = (C[inrow][incol++] + C[inrow++][incol--] 
1198                    + C[inrow][incol++] + C[inrow--][incol++] + bias)/(float)4.0;
1199           bias ^= 3;
1200         }
1201         inrow += 2;
1202         incol = 0;
1203       }
1204       return output;
1205     }
1206   */
1207
1208 }