• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /*
2  *   Colorspace conversions for CUPS.
3  *
4  *   Copyright 2007-2011 by Apple Inc.
5  *   Copyright 1993-2006 by Easy Software Products.
6  *
7  *   The color saturation/hue matrix stuff is provided thanks to Mr. Paul
8  *   Haeberli at "http://www.sgi.com/grafica/matrix/index.html".
9  *
10  *   These coded instructions, statements, and computer programs are the
11  *   property of Apple Inc. and are protected by Federal copyright
12  *   law.  Distribution and use rights are outlined in the file "COPYING"
13  *   which should have been included with this file.
14  *
15  * Contents:
16  *
17  *   cupsImageCMYKToBlack()         - Convert CMYK data to black.
18  *   cupsImageCMYKToCMY()           - Convert CMYK colors to CMY.
19  *   cupsImageCMYKToCMYK()          - Convert CMYK colors to CMYK.
20  *   cupsImageCMYKToRGB()           - Convert CMYK colors to device-dependent
21  *                                    RGB.
22  *   cupsImageCMYKToWhite()         - Convert CMYK colors to luminance.
23  *   cupsImageLut()                 - Adjust all pixel values with the given
24  *                                    LUT.
25  *   cupsImageRGBAdjust()           - Adjust the hue and saturation of the
26  *                                    given RGB colors.
27  *   cupsImageRGBToBlack()          - Convert RGB data to black.
28  *   cupsImageRGBToCMY()            - Convert RGB colors to CMY.
29  *   cupsImageRGBToCMYK()           - Convert RGB colors to CMYK.
30  *   cupsImageRGBToRGB()            - Convert RGB colors to device-dependent
31  *                                    RGB.
32  *   cupsImageRGBToWhite()          - Convert RGB colors to luminance.
33  *   cupsImageSetProfile()          - Set the device color profile.
34  *   cupsImageSetRasterColorSpace() - Set the destination colorspace.
35  *   cupsImageWhiteToBlack()        - Convert luminance colors to black.
36  *   cupsImageWhiteToCMY()          - Convert luminance colors to CMY.
37  *   cupsImageWhiteToCMYK()         - Convert luminance colors to CMYK.
38  *   cupsImageWhiteToRGB()          - Convert luminance data to RGB.
39  *   cupsImageWhiteToWhite()        - Convert luminance colors to device-
40  *                                    dependent luminance.
41  *   cielab()                       - Map CIE Lab transformation...
42  *   huerotate()                    - Rotate the hue, maintaining luminance.
43  *   ident()                        - Make an identity matrix.
44  *   mult()                         - Multiply two matrices.
45  *   rgb_to_lab()                   - Convert an RGB color to CIE Lab.
46  *   rgb_to_xyz()                   - Convert an RGB color to CIE XYZ.
47  *   saturate()                     - Make a saturation matrix.
48  *   xform()                        - Transform a 3D point using a matrix...
49  *   xrotate()                      - Rotate about the x (red) axis...
50  *   yrotate()                      - Rotate about the y (green) axis...
51  *   zrotate()                      - Rotate about the z (blue) axis...
52  *   zshear()                       - Shear z using x and y...
53  */
54 
55 /*
56  * Include necessary headers...
57  */
58 
59 #include "image-private.h"
60 
61 
62 /*
63  * Define some math constants that are required...
64  */
65 
66 #ifndef M_PI
67 #  define M_PI		3.14159265358979323846
68 #endif /* !M_PI */
69 
70 #ifndef M_SQRT2
71 #  define M_SQRT2	1.41421356237309504880
72 #endif /* !M_SQRT2 */
73 
74 #ifndef M_SQRT1_2
75 #  define M_SQRT1_2	0.70710678118654752440
76 #endif /* !M_SQRT1_2 */
77 
78 /*
79  * CIE XYZ whitepoint...
80  */
81 
82 #define D65_X	(0.412453 + 0.357580 + 0.180423)
83 #define D65_Y	(0.212671 + 0.715160 + 0.072169)
84 #define D65_Z	(0.019334 + 0.119193 + 0.950227)
85 
86 
87 /*
88  * Lookup table structure...
89  */
90 
91 typedef int cups_clut_t[3][256];
92 
93 
94 /*
95  * Local globals...
96  */
97 
98 static int		cupsImageHaveProfile = 0;
99 					/* Do we have a color profile? */
100 static int		*cupsImageDensity;
101 					/* Ink/marker density LUT */
102 static cups_clut_t	*cupsImageMatrix;
103 					/* Color transform matrix LUT */
104 static cups_cspace_t	cupsImageColorSpace = CUPS_CSPACE_RGB;
105 					/* Destination colorspace */
106 
107 
108 /*
109  * Local functions...
110  */
111 
112 static float	cielab(float x, float xn);
113 static void	huerotate(float [3][3], float);
114 static void	ident(float [3][3]);
115 static void	mult(float [3][3], float [3][3], float [3][3]);
116 static void	rgb_to_lab(cups_ib_t *val);
117 static void	rgb_to_xyz(cups_ib_t *val);
118 static void	saturate(float [3][3], float);
119 static void	xform(float [3][3], float, float, float, float *, float *, float *);
120 static void	xrotate(float [3][3], float, float);
121 static void	yrotate(float [3][3], float, float);
122 static void	zrotate(float [3][3], float, float);
123 static void	zshear(float [3][3], float, float);
124 
125 
126 /*
127  * 'cupsImageCMYKToBlack()' - Convert CMYK data to black.
128  */
129 
130 void
cupsImageCMYKToBlack(const cups_ib_t * in,cups_ib_t * out,int count)131 cupsImageCMYKToBlack(
132     const cups_ib_t *in,		/* I - Input pixels */
133     cups_ib_t       *out,		/* I - Output pixels */
134     int             count)		/* I - Number of pixels */
135 {
136   int	k;				/* Black value */
137 
138 
139   if (cupsImageHaveProfile)
140     while (count > 0)
141     {
142       k = (31 * in[0] + 61 * in[1] + 8 * in[2]) / 100 + in[3];
143 
144       if (k < 255)
145         *out++ = cupsImageDensity[k];
146       else
147         *out++ = cupsImageDensity[255];
148 
149       in += 4;
150       count --;
151     }
152   else
153     while (count > 0)
154     {
155       k = (31 * in[0] + 61 * in[1] + 8 * in[2]) / 100 + in[3];
156 
157       if (k < 255)
158         *out++ = k;
159       else
160         *out++ = 255;
161 
162       in += 4;
163       count --;
164     }
165 }
166 
167 
168 /*
169  * 'cupsImageCMYKToCMY()' - Convert CMYK colors to CMY.
170  */
171 
172 void
cupsImageCMYKToCMY(const cups_ib_t * in,cups_ib_t * out,int count)173 cupsImageCMYKToCMY(
174     const cups_ib_t *in,		/* I - Input pixels */
175     cups_ib_t       *out,		/* I - Output pixels */
176     int             count)		/* I - Number of pixels */
177 {
178   int	c, m, y, k;			/* CMYK values */
179   int	cc, cm, cy;			/* Calibrated CMY values */
180 
181 
182   if (cupsImageHaveProfile)
183     while (count > 0)
184     {
185       c = *in++;
186       m = *in++;
187       y = *in++;
188       k = *in++;
189 
190       cc = cupsImageMatrix[0][0][c] +
191            cupsImageMatrix[0][1][m] +
192 	   cupsImageMatrix[0][2][y] + k;
193       cm = cupsImageMatrix[1][0][c] +
194            cupsImageMatrix[1][1][m] +
195 	   cupsImageMatrix[1][2][y] + k;
196       cy = cupsImageMatrix[2][0][c] +
197            cupsImageMatrix[2][1][m] +
198 	   cupsImageMatrix[2][2][y] + k;
199 
200       if (cc < 0)
201         *out++ = 0;
202       else if (cc > 255)
203         *out++ = cupsImageDensity[255];
204       else
205         *out++ = cupsImageDensity[cc];
206 
207       if (cm < 0)
208         *out++ = 0;
209       else if (cm > 255)
210         *out++ = cupsImageDensity[255];
211       else
212         *out++ = cupsImageDensity[cm];
213 
214       if (cy < 0)
215         *out++ = 0;
216       else if (cy > 255)
217         *out++ = cupsImageDensity[255];
218       else
219         *out++ = cupsImageDensity[cy];
220 
221       count --;
222     }
223   else
224     while (count > 0)
225     {
226       c = *in++;
227       m = *in++;
228       y = *in++;
229       k = *in++;
230 
231       c += k;
232       m += k;
233       y += k;
234 
235       if (c < 255)
236         *out++ = c;
237       else
238         *out++ = 255;
239 
240       if (m < 255)
241         *out++ = y;
242       else
243         *out++ = 255;
244 
245       if (y < 255)
246         *out++ = y;
247       else
248         *out++ = 255;
249 
250       count --;
251     }
252 }
253 
254 
255 /*
256  * 'cupsImageCMYKToCMYK()' - Convert CMYK colors to CMYK.
257  */
258 
259 void
cupsImageCMYKToCMYK(const cups_ib_t * in,cups_ib_t * out,int count)260 cupsImageCMYKToCMYK(
261     const cups_ib_t *in,		/* I - Input pixels */
262     cups_ib_t       *out,		/* I - Output pixels */
263     int             count)		/* I - Number of pixels */
264 {
265   int	c, m, y, k;			/* CMYK values */
266   int	cc, cm, cy;			/* Calibrated CMY values */
267 
268 
269   if (cupsImageHaveProfile)
270     while (count > 0)
271     {
272       c = *in++;
273       m = *in++;
274       y = *in++;
275       k = *in++;
276 
277       cc = (cupsImageMatrix[0][0][c] +
278             cupsImageMatrix[0][1][m] +
279 	    cupsImageMatrix[0][2][y]);
280       cm = (cupsImageMatrix[1][0][c] +
281             cupsImageMatrix[1][1][m] +
282 	    cupsImageMatrix[1][2][y]);
283       cy = (cupsImageMatrix[2][0][c] +
284             cupsImageMatrix[2][1][m] +
285 	    cupsImageMatrix[2][2][y]);
286 
287       if (cc < 0)
288         *out++ = 0;
289       else if (cc > 255)
290         *out++ = cupsImageDensity[255];
291       else
292         *out++ = cupsImageDensity[cc];
293 
294       if (cm < 0)
295         *out++ = 0;
296       else if (cm > 255)
297         *out++ = cupsImageDensity[255];
298       else
299         *out++ = cupsImageDensity[cm];
300 
301       if (cy < 0)
302         *out++ = 0;
303       else if (cy > 255)
304         *out++ = cupsImageDensity[255];
305       else
306         *out++ = cupsImageDensity[cy];
307 
308       *out++ = cupsImageDensity[k];
309 
310       count --;
311     }
312   else if (in != out)
313   {
314     while (count > 0)
315     {
316       *out++ = *in++;
317       *out++ = *in++;
318       *out++ = *in++;
319       *out++ = *in++;
320 
321       count --;
322     }
323   }
324 }
325 
326 
327 /*
328  * 'cupsImageCMYKToRGB()' - Convert CMYK colors to device-dependent RGB.
329  */
330 
331 void
cupsImageCMYKToRGB(const cups_ib_t * in,cups_ib_t * out,int count)332 cupsImageCMYKToRGB(
333     const cups_ib_t *in,		/* I - Input pixels */
334     cups_ib_t       *out,		/* I - Output pixels */
335     int             count)		/* I - Number of pixels */
336 {
337   int	c, m, y, k;			/* CMYK values */
338   int	cr, cg, cb;			/* Calibrated RGB values */
339 
340 
341   if (cupsImageHaveProfile)
342   {
343     while (count > 0)
344     {
345       c = *in++;
346       m = *in++;
347       y = *in++;
348       k = *in++;
349 
350       cr = cupsImageMatrix[0][0][c] +
351            cupsImageMatrix[0][1][m] +
352            cupsImageMatrix[0][2][y] + k;
353       cg = cupsImageMatrix[1][0][c] +
354            cupsImageMatrix[1][1][m] +
355 	   cupsImageMatrix[1][2][y] + k;
356       cb = cupsImageMatrix[2][0][c] +
357            cupsImageMatrix[2][1][m] +
358 	   cupsImageMatrix[2][2][y] + k;
359 
360       if (cr < 0)
361         *out++ = 255;
362       else if (cr > 255)
363         *out++ = 255 - cupsImageDensity[255];
364       else
365         *out++ = 255 - cupsImageDensity[cr];
366 
367       if (cg < 0)
368         *out++ = 255;
369       else if (cg > 255)
370         *out++ = 255 - cupsImageDensity[255];
371       else
372         *out++ = 255 - cupsImageDensity[cg];
373 
374       if (cb < 0)
375         *out++ = 255;
376       else if (cb > 255)
377         *out++ = 255 - cupsImageDensity[255];
378       else
379         *out++ = 255 - cupsImageDensity[cb];
380 
381       count --;
382     }
383   }
384   else
385   {
386     while (count > 0)
387     {
388       c = 255 - *in++;
389       m = 255 - *in++;
390       y = 255 - *in++;
391       k = *in++;
392 
393       c -= k;
394       m -= k;
395       y -= k;
396 
397       if (c > 0)
398 	*out++ = c;
399       else
400         *out++ = 0;
401 
402       if (m > 0)
403 	*out++ = m;
404       else
405         *out++ = 0;
406 
407       if (y > 0)
408 	*out++ = y;
409       else
410         *out++ = 0;
411 
412       if (cupsImageColorSpace == CUPS_CSPACE_CIELab ||
413           cupsImageColorSpace >= CUPS_CSPACE_ICC1)
414         rgb_to_lab(out - 3);
415       else if (cupsImageColorSpace == CUPS_CSPACE_CIEXYZ)
416         rgb_to_xyz(out - 3);
417 
418       count --;
419     }
420   }
421 }
422 
423 
424 /*
425  * 'cupsImageCMYKToWhite()' - Convert CMYK colors to luminance.
426  */
427 
428 void
cupsImageCMYKToWhite(const cups_ib_t * in,cups_ib_t * out,int count)429 cupsImageCMYKToWhite(
430     const cups_ib_t *in,		/* I - Input pixels */
431     cups_ib_t       *out,		/* I - Output pixels */
432     int             count)		/* I - Number of pixels */
433 {
434   int	w;				/* White value */
435 
436 
437   if (cupsImageHaveProfile)
438   {
439     while (count > 0)
440     {
441       w = 255 - (31 * in[0] + 61 * in[1] + 8 * in[2]) / 100 - in[3];
442 
443       if (w > 0)
444         *out++ = cupsImageDensity[w];
445       else
446         *out++ = cupsImageDensity[0];
447 
448       in += 4;
449       count --;
450     }
451   }
452   else
453   {
454     while (count > 0)
455     {
456       w = 255 - (31 * in[0] + 61 * in[1] + 8 * in[2]) / 100 - in[3];
457 
458       if (w > 0)
459         *out++ = w;
460       else
461         *out++ = 0;
462 
463       in += 4;
464       count --;
465     }
466   }
467 }
468 
469 
470 /*
471  * 'cupsImageLut()' - Adjust all pixel values with the given LUT.
472  */
473 
474 void
cupsImageLut(cups_ib_t * pixels,int count,const cups_ib_t * lut)475 cupsImageLut(cups_ib_t       *pixels,	/* IO - Input/output pixels */
476              int             count,	/* I  - Number of pixels/bytes to adjust */
477              const cups_ib_t *lut)	/* I  - Lookup table */
478 {
479   while (count > 0)
480   {
481     *pixels = lut[*pixels];
482     pixels ++;
483     count --;
484   }
485 }
486 
487 
488 /*
489  * 'cupsImageRGBAdjust()' - Adjust the hue and saturation of the given RGB colors.
490  */
491 
492 void
cupsImageRGBAdjust(cups_ib_t * pixels,int count,int saturation,int hue)493 cupsImageRGBAdjust(cups_ib_t *pixels,	/* IO - Input/output pixels */
494         	   int       count,	/* I - Number of pixels to adjust */
495         	   int       saturation,/* I - Color saturation (%) */
496         	   int       hue)	/* I - Color hue (degrees) */
497 {
498   int			i, j, k;	/* Looping vars */
499   float			mat[3][3];	/* Color adjustment matrix */
500   static int		last_sat = 100,	/* Last saturation used */
501 			last_hue = 0;	/* Last hue used */
502   static cups_clut_t	*lut = NULL;	/* Lookup table for matrix */
503 
504 
505   if (saturation != last_sat || hue != last_hue || !lut)
506   {
507    /*
508     * Build the color adjustment matrix...
509     */
510 
511     ident(mat);
512     saturate(mat, saturation * 0.01);
513     huerotate(mat, (float)hue);
514 
515    /*
516     * Allocate memory for the lookup table...
517     */
518 
519     if (lut == NULL)
520       lut = calloc(3, sizeof(cups_clut_t));
521 
522     if (lut == NULL)
523       return;
524 
525    /*
526     * Convert the matrix into a 3x3 array of lookup tables...
527     */
528 
529     for (i = 0; i < 3; i ++)
530       for (j = 0; j < 3; j ++)
531         for (k = 0; k < 256; k ++)
532           lut[i][j][k] = mat[i][j] * k + 0.5;
533 
534    /*
535     * Save the saturation and hue to compare later...
536     */
537 
538     last_sat = saturation;
539     last_hue = hue;
540   }
541 
542  /*
543   * Adjust each pixel in the given buffer.
544   */
545 
546   while (count > 0)
547   {
548     i = lut[0][0][pixels[0]] +
549         lut[1][0][pixels[1]] +
550         lut[2][0][pixels[2]];
551     if (i < 0)
552       pixels[0] = 0;
553     else if (i > 255)
554       pixels[0] = 255;
555     else
556       pixels[0] = i;
557 
558     i = lut[0][1][pixels[0]] +
559         lut[1][1][pixels[1]] +
560         lut[2][1][pixels[2]];
561     if (i < 0)
562       pixels[1] = 0;
563     else if (i > 255)
564       pixels[1] = 255;
565     else
566       pixels[1] = i;
567 
568     i = lut[0][2][pixels[0]] +
569         lut[1][2][pixels[1]] +
570         lut[2][2][pixels[2]];
571     if (i < 0)
572       pixels[2] = 0;
573     else if (i > 255)
574       pixels[2] = 255;
575     else
576       pixels[2] = i;
577 
578     count --;
579     pixels += 3;
580   }
581 }
582 
583 
584 /*
585  * 'cupsImageRGBToBlack()' - Convert RGB data to black.
586  */
587 
588 void
cupsImageRGBToBlack(const cups_ib_t * in,cups_ib_t * out,int count)589 cupsImageRGBToBlack(
590     const cups_ib_t *in,		/* I - Input pixels */
591     cups_ib_t       *out,		/* I - Output pixels */
592     int             count)		/* I - Number of pixels */
593 {
594   if (cupsImageHaveProfile)
595     while (count > 0)
596     {
597       *out++ = cupsImageDensity[255 - (31 * in[0] + 61 * in[1] + 8 * in[2]) / 100];
598       in += 3;
599       count --;
600     }
601   else
602     while (count > 0)
603     {
604       *out++ = 255 - (31 * in[0] + 61 * in[1] + 8 * in[2]) / 100;
605       in += 3;
606       count --;
607     }
608 }
609 
610 
611 /*
612  * 'cupsImageRGBToCMY()' - Convert RGB colors to CMY.
613  */
614 
615 void
cupsImageRGBToCMY(const cups_ib_t * in,cups_ib_t * out,int count)616 cupsImageRGBToCMY(
617     const cups_ib_t *in,		/* I - Input pixels */
618     cups_ib_t       *out,		/* I - Output pixels */
619     int             count)		/* I - Number of pixels */
620 {
621   int	c, m, y, k;			/* CMYK values */
622   int	cc, cm, cy;			/* Calibrated CMY values */
623 
624 
625   if (cupsImageHaveProfile)
626     while (count > 0)
627     {
628       c = 255 - *in++;
629       m = 255 - *in++;
630       y = 255 - *in++;
631       k = min(c, min(m, y));
632       c -= k;
633       m -= k;
634       y -= k;
635 
636       cc = cupsImageMatrix[0][0][c] +
637            cupsImageMatrix[0][1][m] +
638 	   cupsImageMatrix[0][2][y] + k;
639       cm = cupsImageMatrix[1][0][c] +
640            cupsImageMatrix[1][1][m] +
641 	   cupsImageMatrix[1][2][y] + k;
642       cy = cupsImageMatrix[2][0][c] +
643            cupsImageMatrix[2][1][m] +
644 	   cupsImageMatrix[2][2][y] + k;
645 
646       if (cc < 0)
647         *out++ = 0;
648       else if (cc > 255)
649         *out++ = cupsImageDensity[255];
650       else
651         *out++ = cupsImageDensity[cc];
652 
653       if (cm < 0)
654         *out++ = 0;
655       else if (cm > 255)
656         *out++ = cupsImageDensity[255];
657       else
658         *out++ = cupsImageDensity[cm];
659 
660       if (cy < 0)
661         *out++ = 0;
662       else if (cy > 255)
663         *out++ = cupsImageDensity[255];
664       else
665         *out++ = cupsImageDensity[cy];
666 
667       count --;
668     }
669   else
670     while (count > 0)
671     {
672       c    = 255 - in[0];
673       m    = 255 - in[1];
674       y    = 255 - in[2];
675       k    = min(c, min(m, y));
676 
677       *out++ = (255 - in[1] / 4) * (c - k) / 255 + k;
678       *out++ = (255 - in[2] / 4) * (m - k) / 255 + k;
679       *out++ = (255 - in[0] / 4) * (y - k) / 255 + k;
680       in += 3;
681       count --;
682     }
683 }
684 
685 
686 /*
687  * 'cupsImageRGBToCMYK()' - Convert RGB colors to CMYK.
688  */
689 
690 void
cupsImageRGBToCMYK(const cups_ib_t * in,cups_ib_t * out,int count)691 cupsImageRGBToCMYK(
692     const cups_ib_t *in,		/* I - Input pixels */
693     cups_ib_t       *out,		/* I - Output pixels */
694     int             count)		/* I - Number of pixels */
695 {
696   int	c, m, y, k,			/* CMYK values */
697 	km;				/* Maximum K value */
698   int	cc, cm, cy;			/* Calibrated CMY values */
699 
700 
701   if (cupsImageHaveProfile)
702     while (count > 0)
703     {
704       c = 255 - *in++;
705       m = 255 - *in++;
706       y = 255 - *in++;
707       k = min(c, min(m, y));
708 
709       if ((km = max(c, max(m, y))) > k)
710         k = k * k * k / (km * km);
711 
712       c -= k;
713       m -= k;
714       y -= k;
715 
716       cc = (cupsImageMatrix[0][0][c] +
717             cupsImageMatrix[0][1][m] +
718 	    cupsImageMatrix[0][2][y]);
719       cm = (cupsImageMatrix[1][0][c] +
720             cupsImageMatrix[1][1][m] +
721 	    cupsImageMatrix[1][2][y]);
722       cy = (cupsImageMatrix[2][0][c] +
723             cupsImageMatrix[2][1][m] +
724 	    cupsImageMatrix[2][2][y]);
725 
726       if (cc < 0)
727         *out++ = 0;
728       else if (cc > 255)
729         *out++ = cupsImageDensity[255];
730       else
731         *out++ = cupsImageDensity[cc];
732 
733       if (cm < 0)
734         *out++ = 0;
735       else if (cm > 255)
736         *out++ = cupsImageDensity[255];
737       else
738         *out++ = cupsImageDensity[cm];
739 
740       if (cy < 0)
741         *out++ = 0;
742       else if (cy > 255)
743         *out++ = cupsImageDensity[255];
744       else
745         *out++ = cupsImageDensity[cy];
746 
747       *out++ = cupsImageDensity[k];
748 
749       count --;
750     }
751   else
752     while (count > 0)
753     {
754       c = 255 - *in++;
755       m = 255 - *in++;
756       y = 255 - *in++;
757       k = min(c, min(m, y));
758 
759       if ((km = max(c, max(m, y))) > k)
760         k = k * k * k / (km * km);
761 
762       c -= k;
763       m -= k;
764       y -= k;
765 
766       *out++ = c;
767       *out++ = m;
768       *out++ = y;
769       *out++ = k;
770 
771       count --;
772     }
773 }
774 
775 
776 /*
777  * 'cupsImageRGBToRGB()' - Convert RGB colors to device-dependent RGB.
778  */
779 
780 void
cupsImageRGBToRGB(const cups_ib_t * in,cups_ib_t * out,int count)781 cupsImageRGBToRGB(
782     const cups_ib_t *in,		/* I - Input pixels */
783     cups_ib_t       *out,		/* I - Output pixels */
784     int             count)		/* I - Number of pixels */
785 {
786   int	c, m, y, k;			/* CMYK values */
787   int	cr, cg, cb;			/* Calibrated RGB values */
788 
789 
790   if (cupsImageHaveProfile)
791   {
792     while (count > 0)
793     {
794       c = 255 - *in++;
795       m = 255 - *in++;
796       y = 255 - *in++;
797       k = min(c, min(m, y));
798       c -= k;
799       m -= k;
800       y -= k;
801 
802       cr = cupsImageMatrix[0][0][c] +
803            cupsImageMatrix[0][1][m] +
804            cupsImageMatrix[0][2][y] + k;
805       cg = cupsImageMatrix[1][0][c] +
806            cupsImageMatrix[1][1][m] +
807 	   cupsImageMatrix[1][2][y] + k;
808       cb = cupsImageMatrix[2][0][c] +
809            cupsImageMatrix[2][1][m] +
810 	   cupsImageMatrix[2][2][y] + k;
811 
812       if (cr < 0)
813         *out++ = 255;
814       else if (cr > 255)
815         *out++ = 255 - cupsImageDensity[255];
816       else
817         *out++ = 255 - cupsImageDensity[cr];
818 
819       if (cg < 0)
820         *out++ = 255;
821       else if (cg > 255)
822         *out++ = 255 - cupsImageDensity[255];
823       else
824         *out++ = 255 - cupsImageDensity[cg];
825 
826       if (cb < 0)
827         *out++ = 255;
828       else if (cb > 255)
829         *out++ = 255 - cupsImageDensity[255];
830       else
831         *out++ = 255 - cupsImageDensity[cb];
832 
833       count --;
834     }
835   }
836   else
837   {
838     if (in != out)
839       memcpy(out, in, count * 3);
840 
841     if (cupsImageColorSpace == CUPS_CSPACE_CIELab ||
842         cupsImageColorSpace >= CUPS_CSPACE_ICC1)
843     {
844       while (count > 0)
845       {
846         rgb_to_lab(out);
847 
848 	out += 3;
849 	count --;
850       }
851     }
852     else if (cupsImageColorSpace == CUPS_CSPACE_CIEXYZ)
853     {
854       while (count > 0)
855       {
856         rgb_to_xyz(out);
857 
858 	out += 3;
859 	count --;
860       }
861     }
862   }
863 }
864 
865 
866 /*
867  * 'cupsImageRGBToWhite()' - Convert RGB colors to luminance.
868  */
869 
870 void
cupsImageRGBToWhite(const cups_ib_t * in,cups_ib_t * out,int count)871 cupsImageRGBToWhite(
872     const cups_ib_t *in,		/* I - Input pixels */
873     cups_ib_t       *out,		/* I - Output pixels */
874     int             count)		/* I - Number of pixels */
875 {
876   if (cupsImageHaveProfile)
877   {
878     while (count > 0)
879     {
880       *out++ = 255 - cupsImageDensity[255 - (31 * in[0] + 61 * in[1] + 8 * in[2]) / 100];
881       in += 3;
882       count --;
883     }
884   }
885   else
886   {
887     while (count > 0)
888     {
889       *out++ = (31 * in[0] + 61 * in[1] + 8 * in[2]) / 100;
890       in += 3;
891       count --;
892     }
893   }
894 }
895 
896 
897 /*
898  * 'cupsImageSetProfile()' - Set the device color profile.
899  */
900 
901 void
cupsImageSetProfile(float d,float g,float matrix[3][3])902 cupsImageSetProfile(float d,		/* I - Ink/marker density */
903                     float g,		/* I - Ink/marker gamma */
904                     float matrix[3][3])	/* I - Color transform matrix */
905 {
906   int	i, j, k;			/* Looping vars */
907   float	m;				/* Current matrix value */
908   int	*im;				/* Pointer into cupsImageMatrix */
909 
910 
911  /*
912   * Allocate memory for the profile data...
913   */
914 
915   if (cupsImageMatrix == NULL)
916     cupsImageMatrix = calloc(3, sizeof(cups_clut_t));
917 
918   if (cupsImageMatrix == NULL)
919     return;
920 
921   if (cupsImageDensity == NULL)
922     cupsImageDensity = calloc(256, sizeof(int));
923 
924   if (cupsImageDensity == NULL)
925     return;
926 
927  /*
928   * Populate the profile lookup tables...
929   */
930 
931   cupsImageHaveProfile  = 1;
932 
933   for (i = 0, im = cupsImageMatrix[0][0]; i < 3; i ++)
934     for (j = 0; j < 3; j ++)
935       for (k = 0, m = matrix[i][j]; k < 256; k ++)
936         *im++ = (int)(k * m + 0.5);
937 
938   for (k = 0, im = cupsImageDensity; k < 256; k ++)
939     *im++ = 255.0 * d * pow((float)k / 255.0, g) + 0.5;
940 }
941 
942 
943 /*
944  * 'cupsImageSetRasterColorSpace()' - Set the destination colorspace.
945  */
946 
947 void
cupsImageSetRasterColorSpace(cups_cspace_t cs)948 cupsImageSetRasterColorSpace(
949     cups_cspace_t cs)			/* I - Destination colorspace */
950 {
951  /*
952   * Set the destination colorspace...
953   */
954 
955   cupsImageColorSpace = cs;
956 
957  /*
958   * Don't use color profiles in colorimetric colorspaces...
959   */
960 
961   if (cs == CUPS_CSPACE_CIEXYZ ||
962       cs == CUPS_CSPACE_CIELab ||
963       cs >= CUPS_CSPACE_ICC1)
964     cupsImageHaveProfile = 0;
965 }
966 
967 
968 /*
969  * 'cupsImageWhiteToBlack()' - Convert luminance colors to black.
970  */
971 
972 void
cupsImageWhiteToBlack(const cups_ib_t * in,cups_ib_t * out,int count)973 cupsImageWhiteToBlack(
974     const cups_ib_t *in,		/* I - Input pixels */
975     cups_ib_t       *out,		/* I - Output pixels */
976     int             count)		/* I - Number of pixels */
977 {
978   if (cupsImageHaveProfile)
979     while (count > 0)
980     {
981       *out++ = cupsImageDensity[255 - *in++];
982       count --;
983     }
984   else
985     while (count > 0)
986     {
987       *out++ = 255 - *in++;
988       count --;
989     }
990 }
991 
992 
993 /*
994  * 'cupsImageWhiteToCMY()' - Convert luminance colors to CMY.
995  */
996 
997 void
cupsImageWhiteToCMY(const cups_ib_t * in,cups_ib_t * out,int count)998 cupsImageWhiteToCMY(
999     const cups_ib_t *in,		/* I - Input pixels */
1000     cups_ib_t       *out,		/* I - Output pixels */
1001     int             count)		/* I - Number of pixels */
1002 {
1003   if (cupsImageHaveProfile)
1004     while (count > 0)
1005     {
1006       out[0] = cupsImageDensity[255 - *in++];
1007       out[1] = out[0];
1008       out[2] = out[0];
1009       out += 3;
1010       count --;
1011     }
1012   else
1013     while (count > 0)
1014     {
1015       *out++ = 255 - *in;
1016       *out++ = 255 - *in;
1017       *out++ = 255 - *in++;
1018       count --;
1019     }
1020 }
1021 
1022 
1023 /*
1024  * 'cupsImageWhiteToCMYK()' - Convert luminance colors to CMYK.
1025  */
1026 
1027 void
cupsImageWhiteToCMYK(const cups_ib_t * in,cups_ib_t * out,int count)1028 cupsImageWhiteToCMYK(
1029     const cups_ib_t *in,		/* I - Input pixels */
1030     cups_ib_t       *out,		/* I - Output pixels */
1031     int             count)		/* I - Number of pixels */
1032 {
1033   if (cupsImageHaveProfile)
1034     while (count > 0)
1035     {
1036       *out++ = 0;
1037       *out++ = 0;
1038       *out++ = 0;
1039       *out++ = cupsImageDensity[255 - *in++];
1040       count --;
1041     }
1042   else
1043     while (count > 0)
1044     {
1045       *out++ = 0;
1046       *out++ = 0;
1047       *out++ = 0;
1048       *out++ = 255 - *in++;
1049       count --;
1050     }
1051 }
1052 
1053 
1054 /*
1055  * 'cupsImageWhiteToRGB()' - Convert luminance data to RGB.
1056  */
1057 
1058 void
cupsImageWhiteToRGB(const cups_ib_t * in,cups_ib_t * out,int count)1059 cupsImageWhiteToRGB(
1060     const cups_ib_t *in,		/* I - Input pixels */
1061     cups_ib_t       *out,		/* I - Output pixels */
1062     int             count)		/* I - Number of pixels */
1063 {
1064   if (cupsImageHaveProfile)
1065   {
1066     while (count > 0)
1067     {
1068       out[0] = 255 - cupsImageDensity[255 - *in++];
1069       out[1] = out[0];
1070       out[2] = out[0];
1071       out += 3;
1072       count --;
1073     }
1074   }
1075   else
1076   {
1077     while (count > 0)
1078     {
1079       *out++ = *in;
1080       *out++ = *in;
1081       *out++ = *in++;
1082 
1083       if (cupsImageColorSpace == CUPS_CSPACE_CIELab ||
1084           cupsImageColorSpace >= CUPS_CSPACE_ICC1)
1085         rgb_to_lab(out - 3);
1086       else if (cupsImageColorSpace == CUPS_CSPACE_CIEXYZ)
1087         rgb_to_xyz(out - 3);
1088 
1089       count --;
1090     }
1091   }
1092 }
1093 
1094 
1095 /*
1096  * 'cupsImageWhiteToWhite()' - Convert luminance colors to device-dependent
1097  *                             luminance.
1098  */
1099 
1100 void
cupsImageWhiteToWhite(const cups_ib_t * in,cups_ib_t * out,int count)1101 cupsImageWhiteToWhite(
1102     const cups_ib_t *in,		/* I - Input pixels */
1103     cups_ib_t       *out,		/* I - Output pixels */
1104     int             count)		/* I - Number of pixels */
1105 {
1106   if (cupsImageHaveProfile)
1107     while (count > 0)
1108     {
1109       *out++ = 255 - cupsImageDensity[255 - *in++];
1110       count --;
1111     }
1112   else if (in != out)
1113     memcpy(out, in, count);
1114 }
1115 
1116 
1117 /*
1118  * 'cielab()' - Map CIE Lab transformation...
1119  */
1120 
1121 static float				/* O - Adjusted color value */
cielab(float x,float xn)1122 cielab(float x,				/* I - Raw color value */
1123        float xn)			/* I - Whitepoint color value */
1124 {
1125   float x_xn;				/* Fraction of whitepoint */
1126 
1127 
1128   x_xn = x / xn;
1129 
1130   if (x_xn > 0.008856)
1131     return (cbrt(x_xn));
1132   else
1133     return (7.787 * x_xn + 16.0 / 116.0);
1134 }
1135 
1136 
1137 /*
1138  * 'huerotate()' - Rotate the hue, maintaining luminance.
1139  */
1140 
1141 static void
huerotate(float mat[3][3],float rot)1142 huerotate(float mat[3][3],		/* I - Matrix to append to */
1143           float rot)			/* I - Hue rotation in degrees */
1144 {
1145   float hmat[3][3];			/* Hue matrix */
1146   float lx, ly, lz;			/* Luminance vector */
1147   float xrs, xrc;			/* X rotation sine/cosine */
1148   float yrs, yrc;			/* Y rotation sine/cosine */
1149   float zrs, zrc;			/* Z rotation sine/cosine */
1150   float zsx, zsy;			/* Z shear x/y */
1151 
1152 
1153  /*
1154   * Load the identity matrix...
1155   */
1156 
1157   ident(hmat);
1158 
1159  /*
1160   * Rotate the grey vector into positive Z...
1161   */
1162 
1163   xrs = M_SQRT1_2;
1164   xrc = M_SQRT1_2;
1165   xrotate(hmat,xrs,xrc);
1166 
1167   yrs = -1.0 / sqrt(3.0);
1168   yrc = -M_SQRT2 * yrs;
1169   yrotate(hmat,yrs,yrc);
1170 
1171  /*
1172   * Shear the space to make the luminance plane horizontal...
1173   */
1174 
1175   xform(hmat, 0.3086, 0.6094, 0.0820, &lx, &ly, &lz);
1176   zsx = lx / lz;
1177   zsy = ly / lz;
1178   zshear(hmat, zsx, zsy);
1179 
1180  /*
1181   * Rotate the hue...
1182   */
1183 
1184   zrs = sin(rot * M_PI / 180.0);
1185   zrc = cos(rot * M_PI / 180.0);
1186 
1187   zrotate(hmat, zrs, zrc);
1188 
1189  /*
1190   * Unshear the space to put the luminance plane back...
1191   */
1192 
1193   zshear(hmat, -zsx, -zsy);
1194 
1195  /*
1196   * Rotate the grey vector back into place...
1197   */
1198 
1199   yrotate(hmat, -yrs, yrc);
1200   xrotate(hmat, -xrs, xrc);
1201 
1202  /*
1203   * Append it to the current matrix...
1204   */
1205 
1206   mult(hmat, mat, mat);
1207 }
1208 
1209 
1210 /*
1211  * 'ident()' - Make an identity matrix.
1212  */
1213 
1214 static void
ident(float mat[3][3])1215 ident(float mat[3][3])			/* I - Matrix to identify */
1216 {
1217   mat[0][0] = 1.0;
1218   mat[0][1] = 0.0;
1219   mat[0][2] = 0.0;
1220   mat[1][0] = 0.0;
1221   mat[1][1] = 1.0;
1222   mat[1][2] = 0.0;
1223   mat[2][0] = 0.0;
1224   mat[2][1] = 0.0;
1225   mat[2][2] = 1.0;
1226 }
1227 
1228 
1229 /*
1230  * 'mult()' - Multiply two matrices.
1231  */
1232 
1233 static void
mult(float a[3][3],float b[3][3],float c[3][3])1234 mult(float a[3][3],			/* I - First matrix */
1235      float b[3][3],			/* I - Second matrix */
1236      float c[3][3])			/* I - Destination matrix */
1237 {
1238   int	x, y;				/* Looping vars */
1239   float	temp[3][3];			/* Temporary matrix */
1240 
1241 
1242  /*
1243   * Multiply a and b, putting the result in temp...
1244   */
1245 
1246   for (y = 0; y < 3; y ++)
1247     for (x = 0; x < 3; x ++)
1248       temp[y][x] = b[y][0] * a[0][x] +
1249                    b[y][1] * a[1][x] +
1250                    b[y][2] * a[2][x];
1251 
1252  /*
1253   * Copy temp to c (that way c can be a pointer to a or b).
1254   */
1255 
1256   memcpy(c, temp, sizeof(temp));
1257 }
1258 
1259 
1260 /*
1261  * 'rgb_to_lab()' - Convert an RGB color to CIE Lab.
1262  */
1263 
1264 static void
rgb_to_lab(cups_ib_t * val)1265 rgb_to_lab(cups_ib_t *val)		/* IO - Color value */
1266 {
1267   float	r,				/* Red value */
1268 	g,				/* Green value */
1269 	b,				/* Blue value */
1270 	ciex,				/* CIE X value */
1271 	ciey,				/* CIE Y value */
1272 	ciez,				/* CIE Z value */
1273 	ciey_yn,			/* Normalized luminance */
1274 	ciel,				/* CIE L value */
1275 	ciea,				/* CIE a value */
1276 	cieb;				/* CIE b value */
1277 
1278 
1279  /*
1280   * Convert sRGB to linear RGB...
1281   */
1282 
1283   r = pow((val[0] / 255.0 + 0.055) / 1.055, 2.4);
1284   g = pow((val[1] / 255.0 + 0.055) / 1.055, 2.4);
1285   b = pow((val[2] / 255.0 + 0.055) / 1.055, 2.4);
1286 
1287  /*
1288   * Convert to CIE XYZ...
1289   */
1290 
1291   ciex = 0.412453 * r + 0.357580 * g + 0.180423 * b;
1292   ciey = 0.212671 * r + 0.715160 * g + 0.072169 * b;
1293   ciez = 0.019334 * r + 0.119193 * g + 0.950227 * b;
1294 
1295  /*
1296   * Normalize and convert to CIE Lab...
1297   */
1298 
1299   ciey_yn = ciey / D65_Y;
1300 
1301   if (ciey_yn > 0.008856)
1302     ciel = 116 * cbrt(ciey_yn) - 16;
1303   else
1304     ciel = 903.3 * ciey_yn;
1305 
1306 /*ciel = ciel;*/
1307   ciea = 500 * (cielab(ciex, D65_X) - cielab(ciey, D65_Y));
1308   cieb = 200 * (cielab(ciey, D65_Y) - cielab(ciez, D65_Z));
1309 
1310  /*
1311   * Scale the L value and bias the a and b values by 128 so that all
1312   * numbers are from 0 to 255.
1313   */
1314 
1315   ciel = ciel * 2.55 + 0.5;
1316   ciea += 128.5;
1317   cieb += 128.5;
1318 
1319  /*
1320   * Output 8-bit values...
1321   */
1322 
1323   if (ciel < 0.0)
1324     val[0] = 0;
1325   else if (ciel < 255.0)
1326     val[0] = (int)ciel;
1327   else
1328     val[0] = 255;
1329 
1330   if (ciea < 0.0)
1331     val[1] = 0;
1332   else if (ciea < 255.0)
1333     val[1] = (int)ciea;
1334   else
1335     val[1] = 255;
1336 
1337   if (cieb < 0.0)
1338     val[2] = 0;
1339   else if (cieb < 255.0)
1340     val[2] = (int)cieb;
1341   else
1342     val[2] = 255;
1343 }
1344 
1345 
1346 /*
1347  * 'rgb_to_xyz()' - Convert an RGB color to CIE XYZ.
1348  */
1349 
1350 static void
rgb_to_xyz(cups_ib_t * val)1351 rgb_to_xyz(cups_ib_t *val)		/* IO - Color value */
1352 {
1353   float	r,				/* Red value */
1354 	g,				/* Green value */
1355 	b,				/* Blue value */
1356 	ciex,				/* CIE X value */
1357 	ciey,				/* CIE Y value */
1358 	ciez;				/* CIE Z value */
1359 
1360 
1361  /*
1362   * Convert sRGB to linear RGB...
1363   */
1364 
1365   r = pow((val[0] / 255.0 + 0.055) / 1.055, 2.4);
1366   g = pow((val[1] / 255.0 + 0.055) / 1.055, 2.4);
1367   b = pow((val[2] / 255.0 + 0.055) / 1.055, 2.4);
1368 
1369  /*
1370   * Convert to CIE XYZ...
1371   */
1372 
1373   ciex = 0.412453 * r + 0.357580 * g + 0.180423 * b;
1374   ciey = 0.212671 * r + 0.715160 * g + 0.072169 * b;
1375   ciez = 0.019334 * r + 0.119193 * g + 0.950227 * b;
1376 
1377  /*
1378   * Encode as 8-bit XYZ...
1379   */
1380 
1381   if (ciex < 0.0f)
1382     val[0] = 0;
1383   else if (ciex < 1.1f)
1384     val[0] = (int)(231.8181f * ciex + 0.5);
1385   else
1386     val[0] = 255;
1387 
1388   if (ciey < 0.0f)
1389     val[1] = 0;
1390   else if (ciey < 1.1f)
1391     val[1] = (int)(231.8181f * ciey + 0.5);
1392   else
1393     val[1] = 255;
1394 
1395   if (ciez < 0.0f)
1396     val[2] = 0;
1397   else if (ciez < 1.1f)
1398     val[2] = (int)(231.8181f * ciez + 0.5);
1399   else
1400     val[2] = 255;
1401 }
1402 
1403 
1404 /*
1405  * 'saturate()' - Make a saturation matrix.
1406  */
1407 
1408 static void
saturate(float mat[3][3],float sat)1409 saturate(float mat[3][3],		/* I - Matrix to append to */
1410          float sat)			/* I - Desired color saturation */
1411 {
1412   float	smat[3][3];			/* Saturation matrix */
1413 
1414 
1415   smat[0][0] = (1.0 - sat) * 0.3086 + sat;
1416   smat[0][1] = (1.0 - sat) * 0.3086;
1417   smat[0][2] = (1.0 - sat) * 0.3086;
1418   smat[1][0] = (1.0 - sat) * 0.6094;
1419   smat[1][1] = (1.0 - sat) * 0.6094 + sat;
1420   smat[1][2] = (1.0 - sat) * 0.6094;
1421   smat[2][0] = (1.0 - sat) * 0.0820;
1422   smat[2][1] = (1.0 - sat) * 0.0820;
1423   smat[2][2] = (1.0 - sat) * 0.0820 + sat;
1424 
1425   mult(smat, mat, mat);
1426 }
1427 
1428 
1429 /*
1430  * 'xform()' - Transform a 3D point using a matrix...
1431  */
1432 
1433 static void
xform(float mat[3][3],float x,float y,float z,float * tx,float * ty,float * tz)1434 xform(float mat[3][3],			/* I - Matrix */
1435       float x,				/* I - Input X coordinate */
1436       float y,				/* I - Input Y coordinate */
1437       float z,				/* I - Input Z coordinate */
1438       float *tx,			/* O - Output X coordinate */
1439       float *ty,			/* O - Output Y coordinate */
1440       float *tz)			/* O - Output Z coordinate */
1441 {
1442   *tx = x * mat[0][0] + y * mat[1][0] + z * mat[2][0];
1443   *ty = x * mat[0][1] + y * mat[1][1] + z * mat[2][1];
1444   *tz = x * mat[0][2] + y * mat[1][2] + z * mat[2][2];
1445 }
1446 
1447 
1448 /*
1449  * 'xrotate()' - Rotate about the x (red) axis...
1450  */
1451 
1452 static void
xrotate(float mat[3][3],float rs,float rc)1453 xrotate(float mat[3][3],		/* I - Matrix */
1454         float rs,			/* I - Rotation angle sine */
1455         float rc)			/* I - Rotation angle cosine */
1456 {
1457   float rmat[3][3];			/* I - Rotation matrix */
1458 
1459 
1460   rmat[0][0] = 1.0;
1461   rmat[0][1] = 0.0;
1462   rmat[0][2] = 0.0;
1463 
1464   rmat[1][0] = 0.0;
1465   rmat[1][1] = rc;
1466   rmat[1][2] = rs;
1467 
1468   rmat[2][0] = 0.0;
1469   rmat[2][1] = -rs;
1470   rmat[2][2] = rc;
1471 
1472   mult(rmat, mat, mat);
1473 }
1474 
1475 
1476 /*
1477  * 'yrotate()' - Rotate about the y (green) axis...
1478  */
1479 
1480 static void
yrotate(float mat[3][3],float rs,float rc)1481 yrotate(float mat[3][3],		/* I - Matrix */
1482         float rs,			/* I - Rotation angle sine */
1483         float rc)			/* I - Rotation angle cosine */
1484 {
1485   float rmat[3][3];			/* I - Rotation matrix */
1486 
1487 
1488   rmat[0][0] = rc;
1489   rmat[0][1] = 0.0;
1490   rmat[0][2] = -rs;
1491 
1492   rmat[1][0] = 0.0;
1493   rmat[1][1] = 1.0;
1494   rmat[1][2] = 0.0;
1495 
1496   rmat[2][0] = rs;
1497   rmat[2][1] = 0.0;
1498   rmat[2][2] = rc;
1499 
1500   mult(rmat,mat,mat);
1501 }
1502 
1503 
1504 /*
1505  * 'zrotate()' - Rotate about the z (blue) axis...
1506  */
1507 
1508 static void
zrotate(float mat[3][3],float rs,float rc)1509 zrotate(float mat[3][3],		/* I - Matrix */
1510         float rs,			/* I - Rotation angle sine */
1511         float rc)			/* I - Rotation angle cosine */
1512 {
1513   float rmat[3][3];			/* I - Rotation matrix */
1514 
1515 
1516   rmat[0][0] = rc;
1517   rmat[0][1] = rs;
1518   rmat[0][2] = 0.0;
1519 
1520   rmat[1][0] = -rs;
1521   rmat[1][1] = rc;
1522   rmat[1][2] = 0.0;
1523 
1524   rmat[2][0] = 0.0;
1525   rmat[2][1] = 0.0;
1526   rmat[2][2] = 1.0;
1527 
1528   mult(rmat,mat,mat);
1529 }
1530 
1531 
1532 /*
1533  * 'zshear()' - Shear z using x and y...
1534  */
1535 
1536 static void
zshear(float mat[3][3],float dx,float dy)1537 zshear(float mat[3][3],			/* I - Matrix */
1538        float dx,			/* I - X shear */
1539        float dy)			/* I - Y shear */
1540 {
1541   float smat[3][3];			/* Shear matrix */
1542 
1543 
1544   smat[0][0] = 1.0;
1545   smat[0][1] = 0.0;
1546   smat[0][2] = dx;
1547 
1548   smat[1][0] = 0.0;
1549   smat[1][1] = 1.0;
1550   smat[1][2] = dy;
1551 
1552   smat[2][0] = 0.0;
1553   smat[2][1] = 0.0;
1554   smat[2][2] = 1.0;
1555 
1556   mult(smat, mat, mat);
1557 }
1558 
1559