File Coverage

blib/lib/ICC/Shared.pm
Criterion Covered Total %
statement 91 802 11.3
branch 8 384 2.0
condition 1 168 0.6
subroutine 28 107 26.1
pod 61 68 89.7
total 189 1529 12.3


line stmt bran cond sub pod time code
1             package ICC::Shared;
2              
3 41     41   101821 use strict;
  41         78  
  41         925  
4 41     41   158 use Carp;
  41         57  
  41         1592  
5 41     41   233 use File::Path;
  41         104  
  41         1672  
6 41     41   215 use File::Spec;
  41         62  
  41         851  
7 41     41   165 use List::Util;
  41         72  
  41         1994  
8 41     41   18128 use Math::Matrix;
  41         153561  
  41         1103  
9 41     41   2202 use POSIX qw(:math_h :float_h);
  41         26559  
  41         255  
10 41     41   65395 use Scalar::Util;
  41         98  
  41         1557  
11 41     41   22156 use Storable;
  41         111821  
  41         1873  
12 41     41   19619 use YAML::Tiny;
  41         192427  
  41         2028  
13 41     41   246 use Exporter qw(import);
  41         67  
  41         1637  
14              
15             our $VERSION = 0.58;
16              
17             # revised 2019-11-16
18             #
19             # Copyright © 2004-2020 by William B. Birkett
20              
21             # this module contains common methods, constants and functions
22              
23             # enable static variables
24 41     41   210 use feature 'state';
  41         68  
  41         4504  
25              
26             # constants
27 41     41   235 use constant D50 => [96.42, 100, 82.49]; # ICC D50 XYZ
  41         66  
  41         3231  
28 41     41   229 use constant d50 => [0.9642, 1.0, 0.8249]; # ICC D50 XYZNumber
  41         63  
  41         2221  
29 41     41   195 use constant d50P => [map {$_ * 32768/65535} 0.9642, 1.0, 0.8249]; # ICC D50 16-bit XYZ PCS
  41         69  
  41         72  
  123         2394  
30 41     41   205 use constant PI => atan2(0, -1); # pi
  41         64  
  41         1832  
31 41     41   193 use constant radian => 180/PI; # radian (in degrees)
  41         65  
  41         1886  
32 41     41   240 use constant ln10 => log(10); # natural log of 10
  41         60  
  41         2568  
33 41     41   224 use constant sin16 => sin(16/radian); # sin(16˚)
  41         693  
  41         2249  
34 41     41   246 use constant cos16 => cos(16/radian); # cos(16˚)
  41         62  
  41         5490  
35              
36             BEGIN {
37            
38             # set constant vector elements to be read-only
39 41     41   113 for (@{(D50)}) {Internals::SvREADONLY($_, 1)}
  41         98  
  123         257  
40 41         87 for (@{(d50)}) {Internals::SvREADONLY($_, 1)}
  41         99  
  123         235  
41 41         72 for (@{(d50P)}) {Internals::SvREADONLY($_, 1)}
  41         106  
  123         427133  
42            
43             }
44              
45             # list of constants and functions to export
46             our @EXPORT = qw(D50 d50 d50P PI radian ln10 xyz2Lab Lab2xyz xyz2Lxyz Lxyz2xyz Lab2Lxyz Lxyz2Lab
47             XYZ2Lab Lab2XYZ XYZ2Lxyz Lxyz2XYZ XYZ2xyz xyz2XYZ XYZ2xyY xyY2XYZ x2L L2x dLdx dxdL XYZ2W xyz2dwv
48             xyz2Lab_jac Lab2xyz_jac dEab dEcmc dE94 dE00 dE99 dH dLCh dLab CCT CCT2 bbrad bbxy bbuv XYZ2ucs xy2ucs daylight
49             linear linear_matrix cspline cspline_matrix lagrange lagrange_matrix dotProduct crossProduct
50             flatten clip_struct round s15f162v v2s15f16 makeProfileFolder getICCPath filterPath setFile parse_tokens
51             is_vector is_num_vector is_matrix is_num_matrix);
52              
53             # add POSIX functions/constants
54             push(@EXPORT, @{$POSIX::EXPORT_TAGS{'math_h'}});
55             push(@EXPORT, @{$POSIX::EXPORT_TAGS{'float_h'}});
56              
57             #----------- common methods -----------
58              
59             # copy an object (using Storable)
60             # parameters: ([number_of_copies])
61             # returns: (new_object_references)
62             sub copy {
63              
64             # get parameters
65 0     0 1 0 my ($self, $n) = @_;
66              
67             # set default to one copy
68 0 0       0 $n = 1 if (! defined($n));
69              
70             # return empty if zero copies
71 0 0       0 return() if ($n == 0);
72              
73             # verify number of copies is a positive integer
74 0 0 0     0 (! ref($n) && $n == int($n) && $n > 0) or croak('number of copies not a positive integer');
      0        
75              
76             # verify parameter is a reference
77 0 0       0 (ref($self)) or croak("can't copy class $self");
78              
79             # return copies
80 0 0       0 return(wantarray ? map {Storable::dclone($self)} (1 .. $n) : Storable::dclone($self));
  0         0  
81              
82             }
83              
84             # print object contents
85             # format parameter is an array structure
86             # parameter: ([format])
87             # returns: (string)
88             sub dump {
89              
90             # get parameters
91 0     0 1 0 my ($self, $format) = @_;
92              
93             # verify object has 'sdump' method
94 0 0       0 ($self->can('sdump')) or croak("object lacks 'sdump' method");
95              
96             # get string from 'sdump'
97 0         0 my $s = $self->sdump($format);
98              
99             # print string
100 0         0 print $s, "\n";
101              
102             # return string
103 0         0 return($s);
104              
105             }
106              
107             #------- color encoding functions -------
108              
109             # convert xyz to L*a*b*
110             # parameters: (input_array -or- input_structure)
111             # returns: (output_array -or- output_structure)
112             sub xyz2Lab {
113              
114             # push subroutine reference and number of parameters
115 0     0 1 0 push(@_, \&_xyz2Lab);
116              
117             # convert input
118 0         0 &_convert3;
119              
120             }
121              
122             # convert L*a*b* to xyz
123             # parameters: (input_array -or- input_structure)
124             # returns: (output_array -or- output_structure)
125             sub Lab2xyz {
126              
127             # push subroutine reference and number of parameters
128 0     0 1 0 push(@_, \&_Lab2xyz);
129              
130             # convert input
131 0         0 &_convert3;
132              
133             }
134              
135             # convert xyz to LxLyLz
136             # parameters: (input_array -or- input_structure)
137             # returns: (output_array -or- output_structure)
138             sub xyz2Lxyz {
139              
140             # push subroutine reference and number of parameters
141 0     0 1 0 push(@_, \&_xyz2Lxyz);
142              
143             # convert input
144 0         0 &_convert3;
145              
146             }
147              
148             # convert LxLyLz to xyz
149             # parameters: (input_array -or- input_structure)
150             # returns: (output_array -or- output_structure)
151             sub Lxyz2xyz {
152              
153             # push subroutine reference and number of parameters
154 0     0 1 0 push(@_, \&_Lxyz2xyz);
155              
156             # convert input
157 0         0 &_convert3;
158              
159             }
160              
161             # convert L*a*b* to LxLyLz
162             # parameters: (input_array -or- input_structure)
163             # returns: (output_array -or- output_structure)
164             sub Lab2Lxyz {
165              
166             # push subroutine reference and number of parameters
167 0     0 1 0 push(@_, \&_Lab2Lxyz);
168              
169             # convert input
170 0         0 &_convert3;
171              
172             }
173              
174             # convert LxLyLz to L*a*b*
175             # parameters: (input_array -or- input_structure)
176             # returns: (output_array -or- output_structure)
177             sub Lxyz2Lab {
178              
179             # push subroutine reference and number of parameters
180 0     0 1 0 push(@_, \&_Lxyz2Lab);
181              
182             # convert input
183 0         0 &_convert3;
184              
185             }
186              
187             # convert XYZ to L*a*b*
188             # default illuminant is D50
189             # parameters: (input_array -or- input_structure, [white_point_vector])
190             # returns: (output_array -or- output_structure)
191             sub XYZ2Lab {
192              
193             # push subroutine reference and number of parameters
194 0     0 1 0 push(@_, \&_XYZ2Lab);
195              
196             # convert input
197 0         0 &_convert4;
198              
199             }
200              
201             # convert L*a*b* to XYZ
202             # default illuminant is D50
203             # parameters: (input_array -or- input_structure, [white_point_vector])
204             # returns: (output_array -or- output_structure)
205             sub Lab2XYZ {
206              
207             # push subroutine reference and number of parameters
208 0     0 1 0 push(@_, \&_Lab2XYZ);
209              
210             # convert input
211 0         0 &_convert4;
212              
213             }
214              
215             # convert XYZ to LxLyLz
216             # default illuminant is D50
217             # parameters: (input_array -or- input_structure, [white_point_vector])
218             # returns: (output_array -or- output_structure)
219             sub XYZ2Lxyz {
220              
221             # push subroutine reference and number of parameters
222 0     0 1 0 push(@_, \&_XYZ2Lxyz);
223              
224             # convert input
225 0         0 &_convert4;
226              
227             }
228              
229             # convert LxLyLz to XYZ
230             # default illuminant is D50
231             # parameters: (input_array -or- input_structure, [white_point_vector])
232             # returns: (output_array -or- output_structure)
233             sub Lxyz2XYZ {
234              
235             # push subroutine reference and number of parameters
236 0     0 1 0 push(@_, \&_Lxyz2XYZ);
237              
238             # convert input
239 0         0 &_convert4;
240              
241             }
242              
243             # convert XYZ to xyz
244             # default illuminant is D50
245             # parameters: (input_array -or- input_structure, [white_point_vector])
246             # returns: (output_array -or- output_structure)
247             sub XYZ2xyz {
248              
249             # push subroutine reference and number of parameters
250 0     0 1 0 push(@_, \&_XYZ2xyz);
251              
252             # convert input
253 0         0 &_convert4;
254              
255             }
256              
257             # convert xyz to XYZ
258             # default illuminant is D50
259             # parameters: (input_array -or- input_structure, [white_point_vector])
260             # returns: (output_array -or- output_structure)
261             sub xyz2XYZ {
262              
263             # push subroutine reference and number of parameters
264 0     0 1 0 push(@_, \&_xyz2XYZ);
265              
266             # convert input
267 0         0 &_convert4;
268              
269             }
270              
271             # # convert XYZ to xyY
272             # parameters: (input_array -or- input_structure)
273             # returns: (output_array -or- output_structure)
274             sub XYZ2xyY {
275              
276             # push subroutine reference and number of parameters
277 0     0 1 0 push(@_, \&_XYZ2xyY);
278              
279             # convert input
280 0         0 &_convert3;
281              
282             }
283              
284             # convert xyY to XYZ
285             # parameters: (input_array -or- input_structure)
286             # returns: (output_array -or- output_structure)
287             sub xyY2XYZ {
288              
289             # push subroutine reference and number of parameters
290 0     0 1 0 push(@_, \&_xyY2XYZ);
291              
292             # convert input
293 0         0 &_convert3;
294              
295             }
296              
297             # convert x to L
298             # parameter: (x)
299             # returns: (L)
300             sub x2L {
301              
302             # return L
303 30 50   30 1 66 return(($_[0] > 216/24389) ? 116 * $_[0]**(1/3) - 16 : $_[0] * 24389/27);
304              
305             }
306              
307             # convert L to x
308             # parameter: (L)
309             # returns: (x)
310             sub L2x {
311              
312             # return x
313 120 50   120 1 315 return(($_[0] > 8) ? (($_[0] + 16)/116)**3 : $_[0] * 27/24389);
314              
315             }
316              
317             # compute dL/dx
318             # parameter: (x)
319             # returns: (dL/dx)
320             sub dLdx {
321              
322             # get x
323 0     0 1 0 my $x = shift();
324              
325             # return dL/dx
326 0 0       0 return(($x > 216/24389) ? 116 * $x**(-2/3)/3 : 24389/27);
327              
328             }
329              
330             # compute dx/dL
331             # parameter: (L)
332             # returns: (dx/dL)
333             sub dxdL {
334              
335             # get L
336 0     0 1 0 my $L = shift();
337              
338             # return dx/dL
339 0 0       0 return(($L > 8) ? 3 * (($L + 16)/116)**2/116 : 27/24389);
340              
341             }
342              
343             # compute xyz to L*a*b* Jacobian matrix
344             # parameters: (x, y, z)
345             # returns: (Jacobian_matrix)
346             sub xyz2Lab_jac {
347              
348             # get parameters
349 0     0 1 0 my ($x, $y, $z) = @_;
350              
351             # compute ∂Lx/∂x, ∂Ly/∂y, ∂Lz/∂z values
352 0         0 my $dLx = dLdx($x);
353 0         0 my $dLy = dLdx($y);
354 0         0 my $dLz = dLdx($z);
355              
356             # return Jacobian matrix
357 0         0 return(Math::Matrix->new(
358             [0, $dLy, 0],
359             [$dLx * 500/116, -$dLy * 500/116, 0],
360             [0, $dLy * 200/116, -$dLz * 200/116]
361             ));
362            
363             }
364              
365             # compute L*a*b* to xyz Jacobian matrix
366             # parameters: (L, a, b)
367             # returns: (Jacobian_matrix)
368             sub Lab2xyz_jac {
369              
370             # get parameters
371 0     0 1 0 my ($L, $a, $b) = @_;
372              
373             # compute ∂x/∂Lx, ∂y/∂Ly, ∂z/∂Lz values
374 0         0 my $dx = dxdL($L + 116 * $a/500);
375 0         0 my $dy = dxdL($L);
376 0         0 my $dz = dxdL($L - 116 * $b/200);
377              
378             # return Jacobian matrix
379 0         0 return(Math::Matrix->new(
380             [$dx, $dx * 116/500, 0],
381             [$dy, 0, 0],
382             [$dz, 0, -$dz * 116/200]
383             ));
384            
385             }
386              
387             # convert XYZ to CIE Whiteness
388             # parameters: (X, Y, Z, ref_to_WP_vector)
389             # returns: (W)
390             sub XYZ2W {
391              
392             # get parameters
393 0     0 1 0 my ($X, $Y, $Z, $wtpt) = @_;
394              
395             # compute chromaticity values
396 0         0 my @xyYs = _XYZ2xyY($X, $Y, $Z); # sample
397 0         0 my @xyYw = _XYZ2xyY(@{$wtpt}); # white point
  0         0  
398              
399             # return Whiteness
400 0         0 return($Y + 800 * ($xyYw[0] - $xyYs[0]) + 1700 * ($xyYw[1] - $xyYs[1]));
401              
402             }
403              
404             # convert xyz to density-weighted value
405             # parameters: (x, y, z)
406             # returns: (dw_value)
407             sub xyz2dwv {
408              
409             # get parameters
410 0     0 1 0 my ($x, $y, $z) = @_;
411              
412             # return 0 if any parameter <= 0
413 0 0 0     0 return(0) if ($x <= 0 || $y <= 0 || $x <= 0);
      0        
414              
415             # compute weight values
416 0         0 my $wx = -log($x/($x + $y + $z));
417 0         0 my $wy = -log($y/($x + $y + $z));
418 0         0 my $wz = -log($z/($x + $y + $z));
419              
420             # return density-weighted value
421 0         0 return(($wx * $x + $wy * $y + $wz * $z)/($wx + $wy + $wz));
422              
423             }
424              
425             #-------- color difference functions --------
426              
427             # compute ∆E*ab color difference
428             # parameters: (L*a*b*_1, L*a*b*_2)
429             # returns: (∆E*ab_value)
430             sub dEab {
431              
432             # return ∆E*ab
433 0     0 1 0 return(sqrt(($_[0] - $_[3])**2 + ($_[1] - $_[4])**2 + ($_[2] - $_[5])**2));
434              
435             }
436              
437             # compute CMC color difference
438             # note: this function is not commutative
439             # L*a*b*_1 is the reference, L*a*b*_2 is the sample
440             # optional parameters are l (lightness) and c (chroma)
441             # default values of optional parameters are 2:1 (l:c)
442             # note: this function assumes D65, 10• colorimetry
443             # parameters: (L*a*b*_1, L*a*b*_2, [l, c])
444             # returns: (∆E*cmc_value)
445             sub dEcmc {
446              
447             # local variables
448 0     0 1 0 my ($dL, $C1, $C2, $dC, $dH2, $h1);
449 0         0 my ($F, $T, $Sl, $Sc, $Sh, $l, $c);
450              
451             # compute ∆L*
452 0         0 $dL = $_[0] - $_[3];
453              
454             # compute C1*
455 0         0 $C1 = sqrt($_[1]**2 + $_[2]**2);
456              
457             # compute C2*
458 0         0 $C2 = sqrt($_[4]**2 + $_[5]**2);
459              
460             # compute ∆C*
461 0         0 $dC = $C1 - $C2;
462              
463             # compute ∆H*^2
464 0         0 $dH2 = ($_[1] - $_[4])**2 + ($_[2] - $_[5])**2 - $dC**2;
465              
466             # compute h1
467 0 0 0     0 $h1 = ($_[2] || $_[1]) ? radian * atan2($_[2], $_[1]) : 0;
468              
469             # adjust h1 if negative
470 0 0       0 $h1 += 360 if ($h1 < 0);
471              
472             # compute F
473 0         0 $F = sqrt($C1**4/($C1**4 + 1900));
474              
475             # compute T
476 0 0 0     0 $T = ($h1 >= 164 && $h1 <= 345) ? 0.56 + abs(0.2 * cos(($h1 + 168)/radian)) : 0.36 + abs(0.4 * cos(($h1 + 35)/radian));
477              
478             # compute Sl
479 0 0       0 $Sl = $_[0] < 16 ? 0.511 : 0.040975 * $_[0]/(1 + 0.01765 * $_[0]);
480              
481             # compute Sc
482 0         0 $Sc = 0.0638 * $C1/(1 + 0.0131 * $C1) + 0.638;
483              
484             # compute Sh
485 0         0 $Sh = $Sc * ($F * $T + 1 - $F);
486              
487             # get l
488 0 0       0 $l = defined($_[6]) ? $_[6] : 2;
489              
490             # get c
491 0 0       0 $c = defined($_[7]) ? $_[7] : 1;
492              
493             # return ∆Ecmc
494 0         0 return(sqrt(($dL/($l * $Sl))**2 + ($dC/($c * $Sc))**2 + $dH2/($Sh**2)));
495              
496             }
497              
498             # compute ∆E*94 (graphic arts) color difference
499             # parameters: (L*a*b*_1, L*a*b*_2)
500             # returns: (∆E*94_value)
501             sub dE94 {
502              
503             # local variables
504 0     0 1 0 my ($dL, $C1, $C2, $dC, $dH2, $dH, $Cm);
505              
506             # compute ∆L*
507 0         0 $dL = $_[0] - $_[3];
508              
509             # compute C1*
510 0         0 $C1 = sqrt($_[1]**2 + $_[2]**2);
511              
512             # compute C2*
513 0         0 $C2 = sqrt($_[4]**2 + $_[5]**2);
514              
515             # compute ∆C*
516 0         0 $dC = $C1 - $C2;
517              
518             # compute ∆H*^2
519 0         0 $dH2 = ($_[1] - $_[4])**2 + ($_[2] - $_[5])**2 - $dC**2;
520              
521             # compute geometric mean C*ab (makes function commutative)
522 0         0 $Cm = sqrt($C1 * $C2);
523              
524             # return ∆E*94
525 0         0 return(sqrt($dL**2 + ($dC/(1 + 0.045 * $Cm))**2 + ($dH2/(1 + 0.015 * $Cm)**2)));
526              
527             }
528              
529             # compute ∆E*00 color difference
530             # parameters: (L*a*b*_1, L*a*b*_2)
531             # returns: (∆E*00_value)
532             sub dE00 {
533              
534             # local variables
535 0     0 1 0 my ($dL, $C1, $C2, $Lm, $Cm, $amult, $a1p, $a2p, $C1p, $C2p, $Cmp);
536 0         0 my ($dCp, $h1p, $h2p, $dhp, $dHp, $Hp, $T, $Sl, $Sc, $Sh, $Rt);
537              
538             # compute ∆L*
539 0         0 $dL = $_[3] - $_[0];
540              
541             # compute mean L*
542 0         0 $Lm = ($_[3] + $_[0])/2;
543              
544             # compute C1*
545 0         0 $C1 = sqrt($_[1]**2 + $_[2]**2);
546              
547             # compute C2*
548 0         0 $C2 = sqrt($_[4]**2 + $_[5]**2);
549              
550             # compute mean C*
551 0         0 $Cm = ($C1 + $C2)/2;
552              
553             # compute a' multiplier
554 0         0 $amult = (1 + (1 - sqrt($Cm**7/($Cm**7 + 6103515625)))/2);
555              
556             # compute a1'
557 0         0 $a1p = $_[1] * $amult;
558              
559             # compute a2'
560 0         0 $a2p = $_[4] * $amult;
561              
562             # compute C1'
563 0         0 $C1p = sqrt($a1p**2 + $_[2]**2);
564              
565             # compute C2'
566 0         0 $C2p = sqrt($a2p**2 + $_[5]**2);
567            
568             # compute mean C'
569 0         0 $Cmp = ($C1p + $C2p)/2;
570            
571             # compute ∆C'
572 0         0 $dCp = $C2p - $C1p;
573              
574             # compute h1'
575 0 0 0     0 $h1p = ($_[2] || $a1p) ? radian * atan2($_[2], $a1p) : 0;
576              
577             # adjust h1' if negative
578 0 0       0 $h1p += 360 if ($h1p < 0);
579              
580             # compute h2'
581 0 0 0     0 $h2p = ($_[5] || $a2p) ? radian * atan2($_[5], $a2p) : 0;
582              
583             # adjust h2' if negative
584 0 0       0 $h2p += 360 if ($h2p < 0);
585              
586             # abs(h1' - h2') > 180
587 0 0       0 if (abs($h1p - $h2p) > 180) {
588            
589             # if h2' > h1'
590 0 0       0 if ($h2p > $h1p) {
591            
592             # compute ∆h'
593 0         0 $dhp = $h2p - $h1p - 360;
594            
595             } else {
596            
597             # compute ∆h'
598 0         0 $dhp = $h2p - $h1p + 360;
599            
600             }
601            
602             } else {
603            
604             # compute ∆h'
605 0         0 $dhp = $h2p - $h1p;
606            
607             }
608              
609             # compute ∆H'
610 0         0 $dHp = 2 * sqrt($C1p * $C2p) * sin($dhp/(radian * 2));
611              
612             # if C1' or C2' is zero
613 0 0 0     0 if ($C1p == 0 || $C2p == 0) {
614            
615             # compute H'
616 0         0 $Hp = $h1p + $h2p;
617            
618             } else {
619            
620             # if abs(h1' - h2') > 180
621 0 0       0 if (abs($h1p - $h2p) > 180) {
622            
623             # compute H'
624 0         0 $Hp = ($h1p + $h2p + 360)/2;
625            
626             } else {
627            
628             # compute H'
629 0         0 $Hp = ($h1p + $h2p)/2;
630            
631             }
632            
633             }
634              
635             # compute T
636 0         0 $T = 1 - 0.17 * cos(($Hp - 30)/radian) + 0.24 * cos((2 * $Hp)/radian) + 0.32 * cos((3 * $Hp + 6)/radian) - 0.20 * cos((4 * $Hp - 63)/radian);
637              
638             # compute Sl
639 0         0 $Sl = 1 + (0.015 * ($Lm - 50)**2)/sqrt(20 + ($Lm - 50)**2);
640              
641             # compute Sc
642 0         0 $Sc = 1 + 0.045 * $Cmp;
643              
644             # compute Sh
645 0         0 $Sh = 1 + 0.015 * $Cmp * $T;
646              
647             # compute Rt
648 0         0 $Rt = -2 * sqrt($Cmp**7/($Cmp**7 + 6103515625)) * sin((60 * exp(-(($Hp - 275)/25)**2))/radian);
649              
650             # return ∆E*00 difference
651 0         0 return(sqrt(($dL/$Sl)**2 + ($dCp/$Sc)**2 + ($dHp/$Sh)**2 + $Rt * ($dCp/$Sc) * ($dHp/$Sh)));
652              
653             }
654              
655             # compute DIN 99 color difference
656             # optional parameters are Ke (lightness) and Kch (chroma)
657             # default values of optional parameters are Ke = 1, Kch = 1
658             # parameters: (L*a*b*_1, L*a*b*_2, [Ke, Kch])
659             # returns: (DIN_99_value)
660             sub dE99 {
661              
662             # compute DIN 99 Lab for sample 1
663 0     0 1 0 my @Lab1 = _Lab2DIN99(@_[0 .. 2, 6, 7]);
664              
665             # compute DIN 99 Lab for sample 2
666 0         0 my @Lab2 = _Lab2DIN99(@_[3 .. 7]);
667              
668             # return DIN 99 color difference
669 0         0 return(sqrt(($Lab1[0] - $Lab2[0])**2 + ($Lab1[1] - $Lab2[1])**2 +($Lab1[2] - $Lab2[2])**2));
670              
671             }
672              
673             # compute ∆H color difference
674             # parameters: (L*a*b*_1, L*a*b*_2)
675             # returns: (∆H)
676             sub dH {
677              
678             # compute chroma values
679 0     0 0 0 my $C0 = sqrt($_[1]**2 + $_[2]**2);
680 0         0 my $C1 = sqrt($_[4]**2 + $_[5]**2);
681              
682             # return ∆H
683 0         0 return(sqrt(abs(($_[1] - $_[4])**2 + ($_[2] - $_[5])**2 - ($C1 - $C0)**2)));
684              
685             }
686              
687             # compute ∆L*, ∆Ch (aka ∆F) color differences
688             # as used to determine G7 grayscale compliance
689             # parameters: (L*a*b*_1, L*a*b*_2)
690             # returns: (∆L*, ∆Ch)
691             sub dLCh {
692              
693             # return ∆L*, ∆Ch
694 0     0 1 0 return($_[0] - $_[3], sqrt(($_[1] - $_[4])**2 + ($_[2] - $_[5])**2));
695              
696             }
697              
698             # compute ∆L*, ∆a*, ∆b* color differences
699             # parameters: (L*a*b*_1, L*a*b*_2)
700             # returns: (∆L*, ∆a*, ∆b*)
701             sub dLab {
702              
703             # return ∆L*, ∆a*, ∆b*
704 0     0 1 0 return($_[0] - $_[3], $_[1] - $_[4], $_[2] - $_[5]);
705              
706             }
707              
708             #--------- illuminant functions ---------
709              
710             # correlated color temperature (CCT)
711             # minimizes the UCS u,v error value
712             # the error value should be less than 5E-2
713             # parameters: (x, y)
714             # returns: (CCT, error_value)
715             sub CCT {
716              
717             # get parameters
718 0     0 1 0 my ($x, $y) = @_;
719              
720             # local variables
721 0         0 my ($ut, $vt, $T, $dT, $u, $v, $err, $ux, $vx, $errx, $derv, $derv0);
722              
723             # compute target u,v values
724 0         0 ($ut, $vt) = xy2ucs($x, $y);
725              
726             # compute CCT using McCamy's approximation
727 0         0 $T = CCT2($x, $y);
728              
729             # compute current u,v values
730 0         0 ($u, $v) = bbuv($T);
731              
732             # compute current error
733 0         0 $err = sqrt(($u - $ut)**2 + ($v - $vt)**2);
734              
735             # compute delta values
736 0         0 ($ux, $vx) = bbuv($T + 1E-3);
737              
738             # compute delta error
739 0         0 $errx = sqrt(($ux - $ut)**2 + ($vx - $vt)**2);
740              
741             # compute derr/dT
742 0         0 $derv = ($errx - $err)/1E-3;
743              
744             # initialize delta T
745 0 0       0 $dT = $derv > 0 ? -2 : 2;
746              
747             # optimization loop
748 0         0 for (0 .. 30) {
749            
750             # adjust T value
751 0         0 $T += $dT;
752            
753             # compute current u,v values
754 0         0 ($u, $v) = bbuv($T);
755            
756             # compute current error
757 0         0 $err = sqrt(($u - $ut)**2 + ($v - $vt)**2);
758            
759             # compute delta values
760 0         0 ($ux, $vx) = bbuv($T + 1E-3);
761            
762             # compute delta error
763 0         0 $errx = sqrt(($ux - $ut)**2 + ($vx - $vt)**2);
764            
765             # save previous derr/dT values
766 0         0 $derv0 = $derv;
767            
768             # compute new derr/dT
769 0         0 $derv = ($errx - $err)/1E-3;
770            
771             # quit loop if derr/dT < 1E-9
772 0 0       0 last if (abs($derv) < 1E-9);
773            
774             # adjust delta T if sign of derivative changes
775 0 0       0 $dT /= -2 if (($derv > 0) ^ ($derv0 > 0));
776            
777             }
778              
779             # return CCT
780 0         0 return($T, $err);
781              
782             }
783              
784             # correlated color temperature (CCT)
785             # using McCamy's approximation
786             # parameters: (x, y)
787             # returns: (CCT)
788             sub CCT2 {
789              
790             # get parameters
791 0     0 1 0 my ($x, $y) = @_;
792              
793             # compute n
794 0         0 my $n = ($x - 0.3320)/($y - 0.1858);
795              
796             # return CCT
797 0         0 return(-449 * $n**3 + 3525 * $n**2 - 6823.3 * $n + 5520.33);
798              
799             }
800              
801             # black body radiance (Planck's law)
802             # using constants and formula per CIE 15
803             # wavelength in nm, temperature in degrees Kelvin
804             # parameters: (wavelength, temperature)
805             # returns: (radiance)
806             sub bbrad {
807              
808             # get parameters
809 0     0 1 0 my ($lambda, $T) = @_;
810              
811             # CIE constants
812 0         0 my $c1 = 3.741771E-16; # 2πhc²
813 0         0 my $c2 = 1.4388E-2; # hc/kB
814              
815             # convert wavelength to meters
816 0         0 $lambda *= 1E-9;
817              
818             # return radiance value
819 0         0 return($c1/(PI * $lambda**5 * (exp($c2/($lambda * $T)) - 1)));
820              
821             }
822              
823             # compute chromaticity values of black body radiator
824             # parameter: (temperature)
825             # returns: (x, y)
826             sub bbxy {
827              
828             # get temperature
829 0     0 1 0 my $T = shift();
830              
831             # local variables
832 0         0 my ($b, $X, $Y, $Z, $d);
833              
834             # load CIE color matching functions (YAML format)
835 0         0 state $cmf = YAML::Tiny->read(getICCPath('Data/CIE_cmfs_360-830_x_1.yml'))->[0];
836              
837             # for each wavelength (360 - 830 nm)
838 0         0 for my $i (0 .. 470) {
839            
840             # compute black body reflectance
841 0         0 $b->[$i] = bbrad($i + 360, $T);
842            
843             }
844              
845             # compute colorimetry
846 0         0 $X = dotProduct($cmf->{'CIE1931x'}, $b);
847 0         0 $Y = dotProduct($cmf->{'CIE1931y'}, $b);
848 0         0 $Z = dotProduct($cmf->{'CIE1931z'}, $b);
849              
850             # compute denominator
851 0 0       0 ($d = $X + $Y + $Z) or croak('X + Y + Z = 0 when computing chromaticity');
852              
853             # return x,y
854 0         0 return($X/$d, $Y/$d);
855              
856             }
857              
858             # compute UCS 1960 values of black body radiator
859             # parameter: (temperature)
860             # returns: (u, v)
861             sub bbuv {
862              
863             # get temperature
864 0     0 1 0 my $T = shift();
865              
866             # local variables
867 0         0 my ($b, $X, $Y, $Z, $d);
868              
869             # load CIE color matching functions (YAML format)
870 0         0 state $cmf = YAML::Tiny->read(getICCPath('Data/CIE_cmfs_360-830_x_1.yml'))->[0];
871              
872             # for each wavelength (360 - 830 nm)
873 0         0 for my $i (0 .. 470) {
874            
875             # compute black body reflectance
876 0         0 $b->[$i] = bbrad($i + 360, $T);
877            
878             }
879              
880             # compute colorimetry
881 0         0 $X = dotProduct($cmf->{'CIE1931x'}, $b);
882 0         0 $Y = dotProduct($cmf->{'CIE1931y'}, $b);
883 0         0 $Z = dotProduct($cmf->{'CIE1931z'}, $b);
884              
885             # compute denominator
886 0 0       0 ($d = $X + 15 * $Y + 3 * $Z) or croak('X + 15Y + 3Z = 0 when computing UCS values');
887              
888             # return u,v
889 0         0 return(4 * $X/$d, 6 * $Y/$d);
890              
891             }
892              
893             # convert XYZ to UCS 1960
894             # used to compute color differences for CCT
895             # parameters: (X, Y, Z)
896             # returns: (u, v)
897             sub XYZ2ucs {
898              
899             # get parameters
900 0     0 1 0 my ($X, $Y, $Z) = @_;
901              
902             # compute denominator
903 0 0       0 (my $d = $X + 15 * $Y + 3 * $Z) or croak('X + 15Y + 3Z = 0 when computing UCS values');
904              
905             # return u,v
906 0         0 return(4 * $X/$d, 6 * $Y/$d);
907              
908             }
909              
910             # convert chromaticity to UCS 1960
911             # used to compute color differences for CCT
912             # parameters: (x, y)
913             # returns: (u, v)
914             sub xy2ucs {
915              
916             # get parameters
917 0     0 1 0 my ($x, $y) = @_;
918              
919             # denominator
920 0 0       0 (my $d = 12 * $y - 2 * $x + 3) or croak('12Y - 2X + 3 = 0 when computing UCS values');
921              
922             # return u,v
923 0         0 return(4 * $x/$d, 6 * $y/$d);
924              
925             }
926              
927             # compute daylight SPD
928             # parameter: (color_temperature)
929             # returns: (range, SPD_vector)
930             sub daylight {
931              
932             # get parameter
933 0     0 1 0 my $cct = shift();
934              
935             # local variables
936 0         0 my ($range, $spd, $xD, $yD, $M1, $M2);
937            
938             # verify color temperature
939 0 0 0     0 ($cct >= 4000 && $cct <= 25000) or croak('CCT must be between 4000˚K and 25000˚K');
940              
941             # load CIE daylight Eigenvectors (YAML format)
942 0         0 state $eigen = YAML::Tiny->read(getICCPath('Data/CIE_daylight_300-830_x_5.yml'))->[0];
943              
944             # set range
945 0         0 $range = [300, 830, 5];
946            
947             # if CCT > 7000
948 0 0       0 if ($cct > 7000) {
949            
950             # compute x value
951 0         0 $xD = -2.0064e9/$cct**3 + 1.9018e6/$cct**2 + 0.24748e3/$cct + 0.23704;
952            
953             } else {
954            
955             # compute x value
956 0         0 $xD = -4.6070e9/$cct**3 + 2.9678e6/$cct**2 + 0.09911e3/$cct + 0.244063;
957            
958             }
959              
960             # compute y value
961 0         0 $yD = -3.000 * $xD**2 + 2.870 * $xD - 0.275;
962              
963             # compute M1
964 0         0 $M1 = (-1.3515 - 1.7703 * $xD + 5.9114 * $yD)/(0.0241 + 0.2562 * $xD - 0.7341 * $yD);
965              
966             # compute M2
967 0         0 $M2 = (0.0300 - 31.4424 * $xD + 30.0717 * $yD)/(0.0241 + 0.2562 * $xD - 0.7341 * $yD);
968              
969             # for each wavelength
970 0         0 for my $i (0 .. 106) {
971            
972             # compute spectral power
973 0         0 $spd->[$i] = $eigen->{'CIE_S0'}->[$i] + $M1 * $eigen->{'CIE_S1'}->[$i] + $M2 * $eigen->{'CIE_S2'}->[$i];
974            
975             }
976              
977             # return
978 0         0 return($range, $spd);
979              
980             }
981              
982             #--------- interpolation functions ---------
983              
984             # linear interpolation function
985             # interpolates and/or extrapolates equally spaced data
986             # input/output range structure: [start_nm, end_nm, increment]
987             # optional extrapolation method is 'copy' or 'linear', none returns zeros
988             # parameters: (input_vector, input_range, output_range, [extrapolation_method])
989             # returns: (output_vector)
990             sub linear {
991              
992             # get parameters
993 0     0 1 0 my ($vector_in, $range_in, $range_out, $ext) = @_;
994              
995             # local variables
996 0         0 my ($ix, $ox, $w, $vector_out, $t, $low);
997              
998             # verify input vector size
999 0 0       0 (($ix = $#{$vector_in}) > 0) or croak('input vector must contain two or more elements');
  0         0  
1000              
1001             # verify input range
1002 0 0 0     0 (abs($ix * $range_in->[2] - $range_in->[1] + $range_in->[0]) < 1E-12 && $range_in->[2] > 0) or croak('invalid input range');
1003              
1004             # compute output upper index from range
1005 0 0       0 (($ox = round(($range_out->[1] - $range_out->[0])/$range_out->[2])) > 0) or croak('output upper index must be > 0');
1006              
1007             # verify output range
1008 0 0 0     0 (abs($ox * $range_out->[2] - $range_out->[1] + $range_out->[0]) < 1E-12 && $range_out->[2] > 0) or croak('invalid output range');
1009              
1010             # for each output vector element
1011 0         0 for my $i (0 .. $ox) {
1012            
1013             # compute wavelength
1014 0         0 $w = $range_out->[0] + $i * $range_out->[2];
1015            
1016             # if wavelength < start of source
1017 0 0       0 if ($w < $range_in->[0]) {
    0          
    0          
1018            
1019             # if extrapolation defined
1020 0 0       0 if (defined($ext)) {
1021            
1022             # if linear extrapolation
1023 0 0       0 if ($ext eq 'linear') {
1024            
1025             # compute ratio
1026 0         0 $t = ($w - $range_in->[0])/$range_in->[2];
1027            
1028             # extrapolate
1029 0         0 $vector_out->[$i] = (1 - $t) * $vector_in->[0] + $t * $vector_in->[1];
1030            
1031             # if copy extrapolation
1032             } elsif ('copy') {
1033            
1034             # set element to first knot
1035 0         0 $vector_out->[$i] = $vector_in->[0];
1036            
1037             } else {
1038            
1039             # error
1040             croak('invalid extrapolation type');
1041            
1042             }
1043            
1044             } else {
1045            
1046             # set element to 0
1047 0         0 $vector_out->[$i] = 0;
1048            
1049             }
1050            
1051             # if wavelength > end of source
1052             } elsif ($w > $range_in->[1]) {
1053            
1054             # if extrapolation defined
1055 0 0       0 if (defined($ext)) {
1056            
1057             # if linear extrapolation
1058 0 0       0 if ($ext eq 'linear') {
1059            
1060             # compute ratio
1061 0         0 $t = ($range_in->[1] - $w)/$range_in->[2];
1062            
1063             # extrapolate
1064 0         0 $vector_out->[$i] = (1 - $t) * $vector_in->[-1] + $t * $vector_in->[-2];
1065            
1066             # if copy extrapolation
1067             } elsif ('copy') {
1068            
1069             # set element to last knot
1070 0         0 $vector_out->[$i] = $vector_in->[-1];
1071            
1072             } else {
1073            
1074             # error
1075             croak('invalid extrapolation type');
1076            
1077             }
1078            
1079             } else {
1080            
1081             # set element to 0
1082 0         0 $vector_out->[$i] = 0;
1083            
1084             }
1085            
1086             # if wavelength == end of source
1087             } elsif ($w == $range_in->[1]) {
1088            
1089             # set element to last knot
1090 0         0 $vector_out->[$i] = $vector_in->[-1];
1091            
1092             } else {
1093            
1094             # compute ratio and lower index
1095 0         0 ($t, $low) = POSIX::modf(($w - $range_in->[0])/$range_in->[2]);
1096            
1097             # if ratio is non-zero
1098 0 0       0 if ($t) {
1099            
1100             # interpolate value
1101 0         0 $vector_out->[$i] = (1 - $t) * $vector_in->[$low] + $t * $vector_in->[$low + 1];
1102            
1103             } else {
1104            
1105             # use low source value
1106 0         0 $vector_out->[$i] = $vector_in->[$low];
1107            
1108             }
1109            
1110             }
1111            
1112             }
1113              
1114             # return
1115 0         0 return($vector_out);
1116              
1117             }
1118              
1119             # compute linear interpolation matrix
1120             # input/output range structure: [start_nm, end_nm, increment]
1121             # optional extrapolation method is 'copy' or 'linear', none returns zeros
1122             # parameters: (input_range, output_range, [extrapolation_method])
1123             # returns: (interpolation_matrix)
1124             sub linear_matrix {
1125              
1126             # get parameters
1127 0     0 1 0 my ($range_in, $range_out, $ext) = @_;
1128              
1129             # local variables
1130 0         0 my ($ix, $ox, $mat, $w, $low, $t);
1131              
1132             # check if ICC::Support::Lapack module is loaded
1133 0         0 state $lapack = defined($INC{'ICC/Support/Lapack.pm'});
1134              
1135             # compute input vector size from range
1136 0         0 $ix = round(($range_in->[1] - $range_in->[0])/$range_in->[2]);
1137              
1138             # verify input range
1139 0 0 0     0 ($ix > 0 && abs($ix * $range_in->[2] - $range_in->[1] + $range_in->[0]) < 1E-12 && $range_in->[2] > 0) or croak('invalid input range');
      0        
1140              
1141             # compute output vector size from range
1142 0         0 $ox = round(($range_out->[1] - $range_out->[0])/$range_out->[2]);
1143              
1144             # verify output range
1145 0 0 0     0 ($ox > 0 && abs($ox * $range_out->[2] - $range_out->[1] + $range_out->[0]) < 1E-12 && $range_out->[2] > 0) or croak('invalid output range');
      0        
1146              
1147             # if ICC::Support::Lapack module is loaded
1148 0 0       0 if ($lapack) {
1149            
1150             # make matrix of zeros
1151 0         0 $mat = ICC::Support::Lapack::zeros($ox + 1, $ix + 1);
1152            
1153             } else {
1154            
1155             # make matrix of zeros
1156 0         0 $mat = [map {[(0) x ($ix + 1)]} (0 .. $ox)];
  0         0  
1157            
1158             }
1159              
1160             # for each output
1161 0         0 for my $i (0 .. $ox) {
1162            
1163             # compute wavelength
1164 0         0 $w = $range_out->[0] + $i * $range_out->[2];
1165            
1166             # if wavelength < start of source
1167 0 0       0 if ($w < $range_in->[0]) {
    0          
    0          
1168            
1169             # if extrapolation defined
1170 0 0       0 if (defined($ext)) {
1171            
1172             # if linear extrapolation
1173 0 0       0 if ($ext eq 'linear') {
1174            
1175             # compute ratio
1176 0         0 $t = ($w - $range_in->[0])/$range_in->[2];
1177            
1178             # set elements to ratio
1179 0         0 $mat->[$i][0] = 1 - $t;
1180 0         0 $mat->[$i][1] = $t;
1181            
1182             # if copy extrapolation
1183             } elsif ('copy') {
1184            
1185             # set element to first knot
1186 0         0 $mat->[$i][0] = 1;
1187            
1188             } else {
1189            
1190             # error
1191             croak('invalid extrapolation type');
1192            
1193             }
1194            
1195             }
1196            
1197             # if wavelength > end of source
1198             } elsif ($w > $range_in->[1]) {
1199            
1200             # if extrapolation defined
1201 0 0       0 if (defined($ext)) {
1202            
1203             # if linear extrapolation
1204 0 0       0 if ($ext eq 'linear') {
1205            
1206             # compute ratio
1207 0         0 $t = ($range_in->[1] - $w)/$range_in->[2];
1208            
1209             # set elements to ratio
1210 0         0 $mat->[$i][$ix - 1] = $t;
1211 0         0 $mat->[$i][$ix] = 1 - $t;
1212            
1213             # if copy extrapolation
1214             } elsif ('copy') {
1215            
1216             # set element to last knot
1217 0         0 $mat->[$i][$ix] = 1;
1218            
1219             } else {
1220            
1221             # error
1222             croak('invalid extrapolation type');
1223            
1224             }
1225            
1226             }
1227            
1228             # if wavelength == end of source
1229             } elsif ($w == $range_in->[1]) {
1230            
1231             # set element to last knot
1232 0         0 $mat->[$i][$ix] = 1;
1233            
1234             } else {
1235            
1236             # compute ratio and lower index
1237 0         0 ($t, $low) = POSIX::modf(($w - $range_in->[0])/$range_in->[2]);
1238            
1239             # set elements to ratio
1240 0         0 $mat->[$i][$low + 1] = $t;
1241 0         0 $mat->[$i][$low] = 1 - $t;
1242            
1243             }
1244            
1245             }
1246              
1247             # return
1248 0         0 return(bless($mat, 'Math::Matrix'));
1249              
1250             }
1251              
1252             # cubic spline interpolation function
1253             # interpolates and/or extrapolates equally spaced data
1254             # input/output range structure: [start_nm, end_nm, increment]
1255             # optional extrapolation method is 'copy' or 'linear', none returns zeros
1256             # parameters: (input_vector, input_range, output_range, [extrapolation_method])
1257             # returns: (output_vector)
1258             sub cspline {
1259              
1260             # get parameters
1261 0     0 1 0 my ($vector_in, $range_in, $range_out, $ext) = @_;
1262              
1263             # local variables
1264 0         0 my ($ix, $ox, $mat, $derv, $info, $w, $vector_out);
1265 0         0 my ($low, $t, $tc, $h00, $h01, $h10, $h11);
1266              
1267             # verify input vector size
1268 0 0       0 (($ix = $#{$vector_in}) > 0) or croak('input vector must contain two or more elements');
  0         0  
1269              
1270             # verify input range
1271 0 0 0     0 (abs($ix * $range_in->[2] - $range_in->[1] + $range_in->[0]) < 1E-12 && $range_in->[2] > 0) or croak('invalid input range');
1272              
1273             # compute output upper index from range
1274 0 0       0 (($ox = round(($range_out->[1] - $range_out->[0])/$range_out->[2])) > 0) or croak('output upper index must be > 0');
1275              
1276             # verify output range
1277 0 0 0     0 (abs($ox * $range_out->[2] - $range_out->[1] + $range_out->[0]) < 1E-12 && $range_out->[2] > 0) or croak('invalid output range');
1278              
1279             # check if ICC::Support::Lapack module is loaded
1280 0         0 state $lapack = defined($INC{'ICC/Support/Lapack.pm'});
1281              
1282             # if ICC::Support::Lapack module is loaded
1283 0 0       0 if ($lapack) {
1284            
1285             # for each input element
1286 0         0 for my $i (0 .. $ix) {
1287            
1288             # compute rhs (3 * (y[i+1] - y[i-1]))
1289 0 0       0 $mat->[$i][0] = 3 * ($vector_in->[$i + 1 > $ix ? $ix : $i + 1] - $vector_in->[$i - 1 < 0 ? 0 : $i - 1]);
    0          
1290            
1291             }
1292            
1293             # solve for derivative matrix
1294 0         0 ($info, $derv) = ICC::Support::Lapack::trisolve([(1) x $ix], [2, (4) x ($ix - 1), 2], [(1) x $ix], $mat);
1295            
1296             # otherwise, use Math::Matrix package (slow)
1297             } else {
1298            
1299             # make tri-diagonal matrix object
1300 0         0 $mat = Math::Matrix->tridiagonal([2, (4) x ($ix - 1), 2]);
1301            
1302             # for each input element
1303 0         0 for my $i (0 .. $ix) {
1304            
1305             # append rhs (3 * (y[i+1] - y[i-1]))
1306 0 0       0 $mat->[$i][$ix + 1] = 3 * ($vector_in->[$i + 1 > $ix ? $ix : $i + 1] - $vector_in->[$i - 1 < 0 ? 0 : $i - 1]);
    0          
1307            
1308             }
1309            
1310             # solve for derivative matrix
1311 0         0 $derv = $mat->solve();
1312            
1313             }
1314              
1315             # for each output vector element
1316 0         0 for my $i (0 .. $ox) {
1317            
1318             # compute wavelength
1319 0         0 $w = $range_out->[0] + $i * $range_out->[2];
1320            
1321             # if wavelength < start of source
1322 0 0       0 if ($w < $range_in->[0]) {
    0          
    0          
1323            
1324             # if extrapolation defined
1325 0 0       0 if (defined($ext)) {
1326            
1327             # if linear extrapolation
1328 0 0       0 if ($ext eq 'linear') {
1329            
1330             # set element to linear value
1331 0         0 $vector_out->[$i] = $vector_in->[0] + $derv->[0][0] * ($w - $range_in->[0])/$range_in->[2];
1332            
1333             # if copy extrapolation
1334             } elsif ('copy') {
1335            
1336             # set element to first knot
1337 0         0 $vector_out->[$i] = $vector_in->[0];
1338            
1339             } else {
1340            
1341             # error
1342             croak('invalid extrapolation type');
1343            
1344             }
1345            
1346             } else {
1347            
1348             # set element to 0
1349 0         0 $vector_out->[$i] = 0;
1350            
1351             }
1352            
1353             # if wavelength > end of source
1354             } elsif ($w > $range_in->[1]) {
1355            
1356             # if extrapolation defined
1357 0 0       0 if (defined($ext)) {
1358            
1359             # if linear extrapolation
1360 0 0       0 if ($ext eq 'linear') {
1361            
1362             # set element to linear value
1363 0         0 $vector_out->[$i] = $vector_in->[-1] + $derv->[-1][0] * ($w - $range_in->[1])/$range_in->[2];
1364            
1365             # if copy extrapolation
1366             } elsif ('copy') {
1367            
1368             # set element to last knot
1369 0         0 $vector_out->[$i] = $vector_in->[-1];
1370            
1371             } else {
1372            
1373             # error
1374             croak('invalid extrapolation type');
1375            
1376             }
1377            
1378             } else {
1379            
1380             # set element to 0
1381 0         0 $vector_out->[$i] = 0;
1382            
1383             }
1384            
1385             # if wavelength == end of source
1386             } elsif ($w == $range_in->[1]) {
1387            
1388             # set element to last knot
1389 0         0 $vector_out->[$i] = $vector_in->[-1];
1390            
1391             } else {
1392            
1393             # compute ratio and lower index
1394 0         0 ($t, $low) = POSIX::modf(($w - $range_in->[0])/$range_in->[2]);
1395            
1396             # if ratio is non-zero
1397 0 0       0 if ($t) {
1398            
1399             # if ICC::Support::Lapack module is loaded
1400 0 0       0 if ($lapack) {
1401            
1402             # interpolate value
1403 0         0 $vector_out->[$i] = ICC::Support::Lapack::hermite($t, $vector_in->[$low], $vector_in->[$low + 1], $derv->[$low][0], $derv->[$low + 1][0]);
1404            
1405             } else {
1406            
1407             # compute Hermite coefficients
1408 0         0 $tc = 1 - $t;
1409 0         0 $h00 = (1 + 2 * $t) * $tc * $tc;
1410 0         0 $h01 = 1 - $h00;
1411 0         0 $h10 = $t * $tc * $tc;
1412 0         0 $h11 = -$t * $t * $tc;
1413            
1414             # interpolate value
1415 0         0 $vector_out->[$i] = $h00 * $vector_in->[$low] + $h01 * $vector_in->[$low + 1] + $h10 * $derv->[$low][0] + $h11 * $derv->[$low + 1][0];
1416            
1417             }
1418            
1419             } else {
1420            
1421             # use lower source value
1422 0         0 $vector_out->[$i] = $vector_in->[$low];
1423            
1424             }
1425            
1426             }
1427            
1428             }
1429              
1430             # return
1431 0         0 return($vector_out);
1432              
1433             }
1434              
1435             # compute cubic spline interpolation matrix
1436             # input/output range structure: [start_nm, end_nm, increment]
1437             # optional extrapolation method is 'copy' or 'linear', none returns zeros
1438             # parameters: (input_range, output_range, [extrapolation_method])
1439             # returns: (interpolation_matrix)
1440             sub cspline_matrix {
1441              
1442             # get parameters
1443 0     0 1 0 my ($range_in, $range_out, $ext) = @_;
1444              
1445             # local variables
1446 0         0 my ($ix, $ox, $rhs, $info, $derv, $mat);
1447 0         0 my ($w, $low, $t, $tc, $h00, $h01, $h10, $h11);
1448              
1449             # check if ICC::Support::Lapack module is loaded
1450 0         0 state $lapack = defined($INC{'ICC/Support/Lapack.pm'});
1451              
1452             # compute input vector size from range
1453 0         0 $ix = round(($range_in->[1] - $range_in->[0])/$range_in->[2]);
1454              
1455             # verify input range
1456 0 0 0     0 ($ix > 0 && abs($ix * $range_in->[2] - $range_in->[1] + $range_in->[0]) < 1E-12 && $range_in->[2] > 0) or croak('invalid input range');
      0        
1457              
1458             # compute output vector size from range
1459 0         0 $ox = round(($range_out->[1] - $range_out->[0])/$range_out->[2]);
1460              
1461             # verify output range
1462 0 0 0     0 ($ox > 0 && abs($ox * $range_out->[2] - $range_out->[1] + $range_out->[0]) < 1E-12 && $range_out->[2] > 0) or croak('invalid output range');
      0        
1463              
1464             # if ICC::Support::Lapack module is loaded
1465 0 0       0 if ($lapack) {
1466            
1467             # make rhs matrix (filled with zeros)
1468 0         0 $rhs = ICC::Support::Lapack::zeros($ix + 1);
1469            
1470             # for each row
1471 0         0 for my $i (1 .. $ix - 1) {
1472            
1473             # set diagonal values
1474 0         0 $rhs->[$i - 1][$i] = 3;
1475 0         0 $rhs->[$i + 1][$i] = -3;
1476            
1477             }
1478            
1479             # set endpoint values
1480 0         0 $rhs->[0][0] = -3;
1481 0         0 $rhs->[1][0] = -3;
1482 0         0 $rhs->[$ix][$ix] = 3;
1483 0         0 $rhs->[$ix - 1][$ix] = 3;
1484            
1485             # solve for derivative matrix
1486 0         0 ($info, $derv) = ICC::Support::Lapack::trisolve([(1) x $ix], [2, (4) x ($ix - 1), 2], [(1) x $ix], $rhs);
1487            
1488             # make matrix of zeros
1489 0         0 $mat = ICC::Support::Lapack::zeros($ox + 1, $ix + 1);
1490            
1491             # otherwise, use Math::Matrix package (slow)
1492             } else {
1493            
1494             # make rhs matrix (fill with zeros)
1495 0         0 $rhs = bless([map {[(0) x ($ix + 1)]} (0 .. $ix)], 'Math::Matrix');
  0         0  
1496            
1497             # for each row
1498 0         0 for my $i (1 .. $ix - 1) {
1499            
1500             # set diagonal values
1501 0         0 $rhs->[$i - 1][$i] = 3;
1502 0         0 $rhs->[$i + 1][$i] = -3;
1503            
1504             }
1505            
1506             # set endpoint values
1507 0         0 $rhs->[0][0] = -3;
1508 0         0 $rhs->[1][0] = -3;
1509 0         0 $rhs->[$ix][$ix] = 3;
1510 0         0 $rhs->[$ix - 1][$ix] = 3;
1511            
1512             # solve for derivative matrix
1513 0         0 $derv = Math::Matrix->tridiagonal([2, (4) x ($ix - 1), 2])->concat($rhs)->solve();
1514            
1515             # make matrix of zeros
1516 0         0 $mat = [map {[(0) x ($ix + 1)]} (0 .. $ox)];
  0         0  
1517            
1518             }
1519              
1520             # for each output
1521 0         0 for my $i (0 .. $ox) {
1522            
1523             # compute wavelength
1524 0         0 $w = $range_out->[0] + $i * $range_out->[2];
1525            
1526             # if wavelength < start of source
1527 0 0       0 if ($w < $range_in->[0]) {
    0          
    0          
1528            
1529             # if extrapolation defined
1530 0 0       0 if (defined($ext)) {
1531            
1532             # if linear extrapolation
1533 0 0       0 if ($ext eq 'linear') {
1534            
1535             # compute ratio
1536 0         0 $t = ($w - $range_in->[0])/$range_in->[2];
1537            
1538             # for each input
1539 0         0 for my $j (0 .. $ix) {
1540            
1541             # set element to linear value
1542 0 0       0 $mat->[$i][$j] = ($j == 0 ? 1 : 0) + $derv->[0][$j] * $t;
1543            
1544             }
1545            
1546             # if copy extrapolation
1547             } elsif ('copy') {
1548            
1549             # set element to first knot
1550 0         0 $mat->[$i][0] = 1;
1551            
1552             } else {
1553            
1554             # error
1555             croak('invalid extrapolation type');
1556            
1557             }
1558            
1559             }
1560            
1561             # if wavelength > end of source
1562             } elsif ($w > $range_in->[1]) {
1563            
1564             # if extrapolation defined
1565 0 0       0 if (defined($ext)) {
1566            
1567             # if linear extrapolation
1568 0 0       0 if ($ext eq 'linear') {
1569            
1570             # compute ratio
1571 0         0 $t = ($w - $range_in->[1])/$range_in->[2];
1572            
1573             # for each input
1574 0         0 for my $j (0 .. $ix) {
1575            
1576             # set element to linear value
1577 0 0       0 $mat->[$i][$j] = ($j == $ix ? 1 : 0) + $derv->[-1][$j] * $t;
1578            
1579             }
1580            
1581             # if copy extrapolation
1582             } elsif ('copy') {
1583            
1584             # set element to last knot
1585 0         0 $mat->[$i][$ix] = 1;
1586            
1587             } else {
1588            
1589             # error
1590             croak('invalid extrapolation type');
1591            
1592             }
1593            
1594             }
1595            
1596             # if wavelength == end of source
1597             } elsif ($w == $range_in->[1]) {
1598            
1599             # set element to last knot
1600 0         0 $mat->[$i][$ix] = 1;
1601            
1602             } else {
1603            
1604             # compute ratio and lower index
1605 0         0 ($t, $low) = POSIX::modf(($w - $range_in->[0])/$range_in->[2]);
1606            
1607             # if ratio is non-zero
1608 0 0       0 if ($t) {
1609            
1610             # if ICC::Support::Lapack module is loaded
1611 0 0       0 if ($lapack) {
1612            
1613             # for each input
1614 0         0 for my $j (0 .. $ix) {
1615            
1616             # interpolate value
1617 0 0       0 $mat->[$i][$j] = ICC::Support::Lapack::hermite($t, ($low == $j ? 1 : 0), ($low + 1 == $j ? 1 : 0), $derv->[$low][$j], $derv->[$low + 1][$j]);
    0          
1618            
1619             }
1620            
1621             } else {
1622            
1623             # for each input
1624 0         0 for my $j (0 .. $ix) {
1625            
1626             # compute Hermite coefficients
1627 0         0 $tc = 1 - $t;
1628 0         0 $h00 = (1 + 2 * $t) * $tc * $tc;
1629 0         0 $h01 = 1 - $h00;
1630 0         0 $h10 = $t * $tc * $tc;
1631 0         0 $h11 = -$t * $t * $tc;
1632            
1633             # interpolate partial derivative value
1634 0 0       0 $mat->[$i][$j] = $h00 * ($low == $j ? 1 : 0) + $h01 * ($low + 1 == $j ? 1 : 0) + $h10 * $derv->[$low][$j] + $h11 * $derv->[$low + 1][$j];
    0          
1635            
1636             }
1637            
1638             }
1639            
1640             } else {
1641            
1642             # set to one
1643 0         0 $mat->[$i][$low] = 1;
1644            
1645             }
1646            
1647             }
1648            
1649             }
1650              
1651             # return
1652 0         0 return(bless($mat, 'Math::Matrix'));
1653              
1654             }
1655              
1656             # Lagrange interpolation function (ASTM E 2022)
1657             # interpolates and/or extrapolates equally spaced data
1658             # input/output range structure: [start_nm, end_nm, increment]
1659             # optional extrapolation method is 'copy' or 'linear', none returns zeros
1660             # parameters: (input_vector, input_range, output_range, [extrapolation_method])
1661             # returns: (output_vector)
1662             sub lagrange {
1663              
1664             # get parameters
1665 0     0 1 0 my ($vector_in, $range_in, $range_out, $ext) = @_;
1666              
1667             # local variables
1668 0         0 my ($ix, $ox, $w, $vector_out, $t, $low);
1669              
1670             # verify input vector size
1671 0 0       0 (($ix = $#{$vector_in}) > 0) or croak('input vector must contain two or more elements');
  0         0  
1672              
1673             # verify input range
1674 0 0 0     0 (abs($ix * $range_in->[2] - $range_in->[1] + $range_in->[0]) < 1E-12 && $range_in->[2] > 0) or croak('invalid input range');
1675              
1676             # compute output upper index from range
1677 0 0       0 (($ox = round(($range_out->[1] - $range_out->[0])/$range_out->[2])) > 0) or croak('output upper index must be > 0');
1678              
1679             # verify output range
1680 0 0 0     0 (abs($ox * $range_out->[2] - $range_out->[1] + $range_out->[0]) < 1E-12 && $range_out->[2] > 0) or croak('invalid output range');
1681              
1682             # for each output vector element
1683 0         0 for my $i (0 .. $ox) {
1684            
1685             # compute wavelength
1686 0         0 $w = $range_out->[0] + $i * $range_out->[2];
1687            
1688             # if wavelength < start of source
1689 0 0       0 if ($w < $range_in->[0]) {
    0          
    0          
1690            
1691             # if extrapolation defined
1692 0 0       0 if (defined($ext)) {
1693            
1694             # if linear extrapolation
1695 0 0       0 if ($ext eq 'linear') {
1696            
1697             # compute ratio
1698 0         0 $t = ($range_in->[0] - $w)/$range_in->[2];
1699            
1700             # compute extrapolated value
1701 0         0 $vector_out->[$i] = (1 + 3 * $t/2) * $vector_in->[0];
1702 0         0 $vector_out->[$i] += -2 * $t * $vector_in->[1];
1703 0         0 $vector_out->[$i] += $t * $vector_in->[2]/2;
1704            
1705             # if copy extrapolation
1706             } elsif ('copy') {
1707            
1708             # set element to first knot
1709 0         0 $vector_out->[$i] = $vector_in->[0];
1710            
1711             } else {
1712            
1713             # error
1714             croak('invalid extrapolation type');
1715            
1716             }
1717            
1718             } else {
1719            
1720             # set element to 0
1721 0         0 $vector_out->[$i] = 0;
1722            
1723             }
1724            
1725             # if wavelength > end of source
1726             } elsif ($w > $range_in->[1]) {
1727            
1728             # if extrapolation defined
1729 0 0       0 if (defined($ext)) {
1730            
1731             # if linear extrapolation
1732 0 0       0 if ($ext eq 'linear') {
1733            
1734             # compute ratio
1735 0         0 $t = ($w - $range_in->[1])/$range_in->[2];
1736            
1737             # compute extrapolated value
1738 0         0 $vector_out->[$i] = (1 + 3 * $t/2) * $vector_in->[$ix];
1739 0         0 $vector_out->[$i] += -2 * $t * $vector_in->[$ix - 1];
1740 0         0 $vector_out->[$i] += $t * $vector_in->[$ix - 2]/2;
1741            
1742             # if copy extrapolation
1743             } elsif ('copy') {
1744            
1745             # set element to last knot
1746 0         0 $vector_out->[$i] = $vector_in->[-1];
1747            
1748             } else {
1749            
1750             # error
1751             croak('invalid extrapolation type');
1752            
1753             }
1754            
1755             } else {
1756            
1757             # set element to 0
1758 0         0 $vector_out->[$i] = 0;
1759            
1760             }
1761            
1762             # if wavelength == end of source
1763             } elsif ($w == $range_in->[1]) {
1764            
1765             # set element to last knot
1766 0         0 $vector_out->[$i] = $vector_in->[-1];
1767            
1768             } else {
1769            
1770             # compute ratio and lower index
1771 0         0 ($t, $low) = POSIX::modf(($w - $range_in->[0])/$range_in->[2]);
1772            
1773             # if ratio is zero
1774 0 0       0 if ($t == 0) {
    0          
    0          
1775            
1776             # set element to low knot
1777 0         0 $vector_out->[$i] = $vector_in->[$low];
1778            
1779             # quadratic Lagrange interpolation (lower)
1780             } elsif ($low == 0) {
1781            
1782             # compute Lagrange quadratic interpolated value
1783 0         0 $vector_out->[$i] = ($t - 1) * ($t - 2) * $vector_in->[0]/2;
1784 0         0 $vector_out->[$i] += -$t * ($t - 2) * $vector_in->[1];
1785 0         0 $vector_out->[$i] += ($t - 1) * $t * $vector_in->[2]/2;
1786            
1787             # quadratic Lagrange interpolation (upper)
1788             } elsif ($low == $ix - 1) {
1789            
1790             # complement ratio
1791 0         0 $t = 1 - $t;
1792            
1793             # compute Lagrange quadratic interpolated value
1794 0         0 $vector_out->[$i] = ($t - 1) * ($t - 2) * $vector_in->[$ix]/2;
1795 0         0 $vector_out->[$i] += -$t * ($t - 2) * $vector_in->[$ix - 1];
1796 0         0 $vector_out->[$i] += ($t - 1) * $t * $vector_in->[$ix - 2]/2;
1797            
1798             # cubic Lagrange interpolation
1799             } else {
1800            
1801             # increment ratio
1802 0         0 $t = 1 + $t;
1803            
1804             # compute Lagrange cubic interpolated value
1805 0         0 $vector_out->[$i] = ($t - 1) * ($t - 2) * ($t - 3) * $vector_in->[$low - 1]/-6;
1806 0         0 $vector_out->[$i] += $t * ($t - 2) * ($t - 3) * $vector_in->[$low]/2;
1807 0         0 $vector_out->[$i] += ($t - 1) * $t * ($t - 3) * $vector_in->[$low + 1]/-2;
1808 0         0 $vector_out->[$i] += ($t - 1) * ($t - 2) * $t * $vector_in->[$low + 2]/6;
1809            
1810             }
1811            
1812             }
1813            
1814             }
1815              
1816             # return
1817 0         0 return($vector_out);
1818              
1819             }
1820              
1821             # compute Lagrange interpolation matrix (ASTM E 2022)
1822             # input/output range structure: [start_nm, end_nm, increment]
1823             # optional extrapolation method is 'copy' or 'linear', none returns zeros
1824             # parameters: (input_range, output_range, [extrapolation_method])
1825             # returns: (interpolation_matrix)
1826             sub lagrange_matrix {
1827              
1828             # get parameters
1829 0     0 1 0 my ($range_in, $range_out, $ext) = @_;
1830              
1831             # local variables
1832 0         0 my ($ix, $ox, $mat, $w, $low, $t);
1833              
1834             # check if ICC::Support::Lapack module is loaded
1835 0         0 state $lapack = defined($INC{'ICC/Support/Lapack.pm'});
1836              
1837             # compute input vector size from range
1838 0         0 $ix = round(($range_in->[1] - $range_in->[0])/$range_in->[2]);
1839              
1840             # verify input range
1841 0 0 0     0 ($ix > 0 && abs($ix * $range_in->[2] - $range_in->[1] + $range_in->[0]) < 1E-12 && $range_in->[2] > 0) or croak('invalid input range');
      0        
1842              
1843             # compute output vector size from range
1844 0         0 $ox = round(($range_out->[1] - $range_out->[0])/$range_out->[2]);
1845              
1846             # verify output range
1847 0 0 0     0 ($ox > 0 && abs($ox * $range_out->[2] - $range_out->[1] + $range_out->[0]) < 1E-12 && $range_out->[2] > 0) or croak('invalid output range');
      0        
1848              
1849             # if ICC::Support::Lapack module is loaded
1850 0 0       0 if ($lapack) {
1851            
1852             # make matrix of zeros
1853 0         0 $mat = ICC::Support::Lapack::zeros($ox + 1, $ix + 1);
1854            
1855             } else {
1856            
1857             # make matrix of zeros
1858 0         0 $mat = [map {[(0) x ($ix + 1)]} (0 .. $ox)];
  0         0  
1859            
1860             }
1861              
1862             # for each output
1863 0         0 for my $i (0 .. $ox) {
1864            
1865             # compute wavelength
1866 0         0 $w = $range_out->[0] + $i * $range_out->[2];
1867            
1868             # if wavelength < start of source
1869 0 0       0 if ($w < $range_in->[0]) {
    0          
    0          
1870            
1871             # if extrapolation defined
1872 0 0       0 if (defined($ext)) {
1873            
1874             # if linear extrapolation
1875 0 0       0 if ($ext eq 'linear') {
1876            
1877             # compute ratio
1878 0         0 $t = ($range_in->[0] - $w)/$range_in->[2];
1879            
1880             # set matrix elements
1881 0         0 $mat->[$i][0] = 1 + 3 * $t/2;
1882 0         0 $mat->[$i][1] = -2 * $t;
1883 0         0 $mat->[$i][2] = $t/2;
1884            
1885             # if copy extrapolation
1886             } elsif ('copy') {
1887            
1888             # set element to first knot
1889 0         0 $mat->[$i][0] = 1;
1890            
1891             } else {
1892            
1893             # error
1894             croak('invalid extrapolation type');
1895            
1896             }
1897            
1898             }
1899            
1900             # if wavelength > end of source
1901             } elsif ($w > $range_in->[1]) {
1902            
1903             # if extrapolation defined
1904 0 0       0 if (defined($ext)) {
1905            
1906             # if linear extrapolation
1907 0 0       0 if ($ext eq 'linear') {
1908            
1909             # compute ratio
1910 0         0 $t = ($w - $range_in->[1])/$range_in->[2];
1911            
1912             # set matrix elements
1913 0         0 $mat->[$i][$ix] = 1 + 3 * $t/2;
1914 0         0 $mat->[$i][$ix - 1] = -2 * $t;
1915 0         0 $mat->[$i][$ix - 2] = $t/2;
1916            
1917             # if copy extrapolation
1918             } elsif ('copy') {
1919            
1920             # set element to last knot
1921 0         0 $mat->[$i][$ix] = 1;
1922            
1923             } else {
1924            
1925             # error
1926             croak('invalid extrapolation type');
1927            
1928             }
1929            
1930             }
1931            
1932             # if wavelength == end of source
1933             } elsif ($w == $range_in->[1]) {
1934            
1935             # set element to last knot
1936 0         0 $mat->[$i][$ix] = 1;
1937            
1938             } else {
1939            
1940             # compute ratio and lower index
1941 0         0 ($t, $low) = POSIX::modf(($w - $range_in->[0])/$range_in->[2]);
1942            
1943             # if ratio is zero
1944 0 0       0 if ($t == 0) {
    0          
    0          
1945            
1946             # set element to low knot
1947 0         0 $mat->[$i][$low] = 1;
1948            
1949             # quadratic Lagrange interpolation (lower)
1950             } elsif ($low == 0) {
1951            
1952             # compute Lagrange quadratic interpolation coefficients
1953 0         0 $mat->[$i][0] = ($t - 1) * ($t - 2)/2;
1954 0         0 $mat->[$i][1] = -$t * ($t - 2);
1955 0         0 $mat->[$i][2] = ($t - 1) * $t/2;
1956            
1957             # quadratic Lagrange interpolation (upper)
1958             } elsif ($low == $ix - 1) {
1959            
1960             # complement ratio
1961 0         0 $t = 1 - $t;
1962            
1963             # compute Lagrange quadratic interpolation coefficients
1964 0         0 $mat->[$i][$ix] = ($t - 1) * ($t - 2)/2;
1965 0         0 $mat->[$i][$ix - 1] = -$t * ($t - 2);
1966 0         0 $mat->[$i][$ix - 2] = ($t - 1) * $t/2;
1967            
1968             # cubic Lagrange interpolation
1969             } else {
1970            
1971             # increment ratio
1972 0         0 $t = 1 + $t;
1973            
1974             # compute Lagrange cubic interpolation coefficients
1975 0         0 $mat->[$i][$low - 1] = ($t - 1) * ($t - 2) * ($t - 3)/-6;
1976 0         0 $mat->[$i][$low] = $t * ($t - 2) * ($t - 3)/2;
1977 0         0 $mat->[$i][$low + 1] = ($t - 1) * $t * ($t - 3)/-2;
1978 0         0 $mat->[$i][$low + 2] = ($t - 1) * ($t - 2) * $t/6;
1979            
1980             }
1981            
1982             }
1983            
1984             }
1985              
1986             # return
1987 0         0 return(bless($mat, 'Math::Matrix'));
1988              
1989             }
1990              
1991             #--------- vector functions ---------
1992              
1993             # compute vector dot product
1994             # parameters: (vector_1, vector_2)
1995             # vectors must have equal dimensions
1996             # returns: (dot_product)
1997             sub dotProduct {
1998              
1999             # check if ICC::Support::Lapack module is loaded
2000 0     0 1 0 state $lapack = defined($INC{'ICC/Support/Lapack.pm'});
2001              
2002             # if ICC::Support::Lapack module is loaded
2003 0 0       0 if ($lapack) {
2004            
2005             # return dot product
2006 0         0 return(ICC::Support::Lapack::dot(@_));
2007            
2008             } else {
2009            
2010             # verify vector sizes match
2011 0 0       0 ($#{$_[0]} == $#{$_[1]}) or croak('vector size mismatch');
  0         0  
  0         0  
2012            
2013             # initialize sum
2014 0         0 my $sum = 0;
2015            
2016             # for each pair of array elements
2017 0         0 for (0 .. $#{$_[0]}) {
  0         0  
2018            
2019             # accumulate product
2020 0         0 $sum += $_[0][$_] * $_[1][$_];
2021            
2022             }
2023            
2024             # return dot product
2025 0         0 return($sum);
2026            
2027             }
2028            
2029             }
2030              
2031             # compute vector cross product
2032             # parameters: (vector_1, vector_2)
2033             # vectors must be three-dimensional
2034             # returns: (cross_product_vector)
2035             sub crossProduct {
2036              
2037             # verify vectors are three-dimensional
2038 0 0 0 0 1 0 ($#{$_[0]} == 2 && $#{$_[1]} == 2) or croak('vectors not three-dimensional');
  0         0  
  0         0  
2039              
2040             # return cross product
2041 0         0 return([$_[0][1] * $_[1][2] - $_[0][2] * $_[1][1], $_[0][2] * $_[1][0] - $_[0][0] * $_[1][2], $_[0][0] * $_[1][1] - $_[0][1] * $_[1][0]]);
2042              
2043             }
2044              
2045             #--------- utility functions ---------
2046              
2047             # flatten slice structure
2048             # structure may contain scalars, arrays or Math::Matrix objects
2049             # the structure is converted to an array of scalars
2050             # parameters: (ref_to_slice_structure)
2051             # returns: (reference_to_flattened_slice)
2052             sub flatten {
2053              
2054             # local variables
2055 612     612 1 665 my (@out);
2056              
2057             # while input array
2058 612         963 while (@_) {
2059            
2060             # shift input value
2061 3945         4021 my $v = shift();
2062            
2063             # if undefined
2064 3945 50 33     5725 if (! defined($v)) {
    100          
    50          
2065            
2066             # warn
2067 0         0 carp('slice contains undefined value');
2068            
2069             # if scalar
2070             } elsif (! ref($v)) {
2071            
2072             # push on output
2073 3672         5632 push(@out, $v);
2074            
2075             # if array reference or Math::Matrix object
2076             } elsif (ref($v) eq 'ARRAY' || UNIVERSAL::isa($v, 'Math::Matrix')) {
2077            
2078             # dereference and push on input
2079 273         293 push(@_, @{$v});
  273         646  
2080            
2081             } else {
2082            
2083             # error
2084 0         0 croak('illegal slice structure');
2085            
2086             }
2087            
2088             }
2089              
2090             # return reference to flattened slice
2091 612         1295 return(\@out);
2092              
2093             }
2094              
2095             # recursive clipping function
2096             # structure may contain scalars, arrays or Math::Matrix objects
2097             # clips scalar elements < 0.0 or > 1.0
2098             # parameter: (reference_to_structure)
2099             sub clip_struct {
2100              
2101             # if array reference or Math::Matrix object
2102 0 0 0 0 1 0 if (ref($_[0]) eq 'ARRAY' || UNIVERSAL::isa($_[0], 'Math::Matrix')) {
2103            
2104             # for each array element
2105 0         0 for (@{$_[0]}) {
  0         0  
2106            
2107             # if a reference
2108 0 0       0 if (ref()) {
2109            
2110             # call myself
2111 0         0 clip_struct($_)
2112            
2113             } else {
2114            
2115             # clip scalar
2116 0 0       0 $_ = 0.0 if ($_ < 0.0);
2117 0 0       0 $_ = 1.0 if ($_ > 1.0);
2118            
2119             }
2120            
2121             }
2122            
2123             } else {
2124            
2125             # error
2126 0         0 croak('error clipping structure');
2127            
2128             }
2129            
2130             }
2131              
2132             # round off number to nearest integer
2133             # parameter: (numeric_value)
2134             # returns: (integer_value)
2135             sub round {
2136              
2137             # return rounded value
2138 0 0   0 1 0 $_[0] > 0 ? int($_[0] + 0.5) : int($_[0] - 0.5)
2139              
2140             }
2141              
2142             # convert array of s15Fixed16Number values to numeric
2143             # parameters: (input_array)
2144             # returns: (converted_array)
2145             sub s15f162v {
2146              
2147             # get parameters
2148 15     15 1 24 my (@in) = @_;
2149              
2150             # return converted array
2151 15 50       25 return(map {($_ & 0x80000000) ? $_/65536 - 65536 : $_/65536} @in);
  45         128  
2152              
2153             }
2154              
2155             # convert array of numeric values to s15Fixed16Number
2156             # parameters: (input_array)
2157             # returns: (converted_array)
2158             sub v2s15f16 {
2159              
2160             # get parameters
2161 11     11 1 22 my (@in) = @_;
2162              
2163             # return converted array
2164 11 50       18 return(map {($_ < 0) ? ($_ + 65536) * 65536 : $_ * 65536} @in);
  45         107  
2165              
2166             }
2167              
2168             # make profile folder
2169             # makes '~profiles' folder, and an alias
2170             # default alias folder is '~/Library/ColorSync/Profiles/'
2171             # customer and/or job will be undefined if not in path
2172             # parameters: (file/folder_path, [alias_folder_path])
2173             # returns: (profiles_folder_path, directory_segs, customer, job)
2174             sub makeProfileFolder {
2175              
2176             # get parameters
2177 0     0 1 0 my ($path, $alias) = @_;
2178              
2179             # local variables
2180 0         0 my ($vol, $dir, $file, @dsegs, $dref);
2181 0         0 my (@jobs, $i, $cust, $job, $folder, $sym);
2182              
2183             # verify parameter is a valid path
2184 0 0       0 (-e $path) or croak('invalid path parameter');
2185              
2186             # split the absolute path (note use of 'no_file' flag)
2187 0         0 ($vol, $dir, $file) = File::Spec->splitpath(File::Spec->rel2abs($path), -d $path);
2188              
2189             # split directory into segments
2190 0         0 @dsegs = File::Spec->splitdir($dir);
2191              
2192             # remove any trailing path delimiters
2193 0   0     0 while (@dsegs > 1 && $dsegs[-1] eq '') {pop(@dsegs)};
  0         0  
2194            
2195             # make copy of directory segments array
2196 0         0 $dref = [@dsegs];
2197              
2198             # find indices of 'jobs' directory segments (if any)
2199 0         0 @jobs = grep {lc($dsegs[$_]) eq 'jobs'} (0 .. $#dsegs);
  0         0  
2200              
2201             # if customer directory segment exists
2202 0 0 0     0 if (@jobs && ($i = $jobs[0] + 1) <= $#dsegs) {
2203            
2204             # set customer
2205 0         0 $cust = $dsegs[$i];
2206            
2207             # replace spaces with underscores
2208 0         0 $cust =~ s/ /_/g;
2209            
2210             }
2211              
2212             # if job number directory segment exists
2213 0 0 0     0 if (@jobs && ($i = $jobs[0] + 2) <= $#dsegs) {
2214            
2215             # match job number
2216 0 0       0 if ($dsegs[$i] =~ m/^(\d+)/) {
2217            
2218             # set job number
2219 0         0 $job = $1;
2220            
2221             } else {
2222            
2223             # set job number to segment
2224 0         0 $job = $dsegs[$i];
2225            
2226             # replace spacer with underscores
2227 0         0 $job =~ s/ /_/g;
2228            
2229             }
2230            
2231             # truncate directory segments array
2232 0         0 splice(@dsegs, $i + 1);
2233            
2234             }
2235              
2236             # add '~profiles' segment
2237 0         0 push(@dsegs, '~profiles');
2238              
2239             # join directory segments
2240 0         0 $dir = File::Spec->catdir(@dsegs);
2241              
2242             # make profile folder path
2243 0         0 $folder = File::Spec->catpath($vol, $dir, '');
2244              
2245             # make profile folder
2246 0 0       0 mkdir($folder) if (! -d $folder);
2247              
2248             # make profile folder alias name (note colons appear as '/' in OS X paths)
2249 0         0 $sym = 'alias_to_' . join(':', @dsegs);
2250              
2251             # if alias parameter provided
2252 0 0       0 if (defined($alias)) {
2253            
2254             # verify a valid directory
2255 0 0       0 (-d $alias) or croak('invalid alias parameter');
2256            
2257             # append alias folder name
2258 0         0 $alias .= "/$sym";
2259            
2260             } else {
2261            
2262             # append alias folder name
2263 0         0 $alias = "$ENV{'HOME'}/Library/ColorSync/Profiles/$sym";
2264            
2265             }
2266              
2267             # make profile folder alias (Unix symbolic link)
2268 0 0       0 symlink($folder, $alias) if (! -d $alias);
2269              
2270             # return
2271 0         0 return($folder, $dref, $cust, $job);
2272              
2273             }
2274              
2275             # get 'ICC' directory path
2276             # optional parameter is appended to path
2277             # parameter: ([file/folder])
2278             # returns: (path)
2279             sub getICCPath {
2280              
2281             # get path to this module
2282 0     0 1 0 my $path = $INC{'ICC/Shared.pm'};
2283              
2284             # if file/folder name supplied
2285 0 0       0 if (@_) {
2286            
2287             # replace module file name
2288 0         0 $path =~ s/Shared.pm$/$_[0]/;
2289            
2290             } else {
2291            
2292             # remove module file name
2293 0         0 $path =~ s/Shared.pm$//;
2294            
2295             }
2296              
2297             # return path (may be invalid)
2298 0         0 return($path);
2299              
2300             }
2301              
2302             # filter file path
2303             # replaces '~' with user home directory
2304             # removes '\' characters (Unix) or replaces '/' with '\' (Windows)
2305             # parameter: (file_path)
2306             sub filterPath {
2307              
2308             # return if path undefined or a reference
2309 0 0 0 0 1 0 return() if (! defined($_[0]) || ref($_[0]));
2310              
2311             # if Windows
2312 0 0       0 if ($^O eq 'MSWin32') {
2313            
2314             # replace '~'
2315 0         0 $_[0] =~ s/^~/$ENV{'USERPROFILE'}/;
2316            
2317             # replace '/' with '\'
2318 0         0 $_[0] =~ s/\//\\/g;
2319            
2320             } else {
2321            
2322             # replace '~'
2323 0         0 $_[0] =~ s/^~/$ENV{'HOME'}/;
2324              
2325             # remove '\'
2326 0         0 $_[0] =~ s/\\//g;
2327            
2328             }
2329              
2330             # return
2331 0         0 return();
2332              
2333             }
2334              
2335             # set macOS file creator and type
2336             # calls the 'SetFile' utility, if present
2337             # parameters: (file_path, creator, type)
2338             sub setFile {
2339              
2340             # get parameters
2341 0     0 1 0 my ($path, $creator, $type) = @_;
2342              
2343             # replace '~'
2344 0         0 $path =~ s/^~/$ENV{'HOME'}/;
2345              
2346             # escape spaces and special characters
2347 0         0 $path =~ s/([^\w\-\/+.@])/\\$1/g;
2348              
2349             # if SetFile in /Library/Developer directory
2350 0 0       0 if (-f '/Library/Developer/CommandLineTools/usr/bin/SetFile') {
2351            
2352             # make system call
2353 0         0 qx(/Library/Developer/CommandLineTools/usr/bin/SetFile -c $creator -t $type $path);
2354            
2355             }
2356              
2357             }
2358              
2359             # parse token string
2360             # string is split into tokens
2361             # token parameters are processed with 'eval'
2362             # returned array contains tokens and their parameters
2363             # parameters: (token_string)
2364             # returns: (ref_to_array)
2365             sub parse_tokens {
2366              
2367             # get parameter
2368 0     0 1 0 my $str = shift();
2369              
2370             # local variables
2371 0         0 my ($nest, $j, $char, @px, $sub, @pars, @eval);
2372              
2373             # return if parameter undefined
2374 0 0       0 return([]) if (! defined($str));
2375              
2376             # if parameter not a scalar
2377 0 0       0 if (ref($str)) {
2378            
2379             # warn
2380 0         0 carp("invalid 'parse_token' parameter\n");
2381            
2382             # return
2383 0         0 return([]);
2384            
2385             }
2386              
2387             # initialize loop
2388 0         0 $nest = 0;
2389 0         0 $j = 0;
2390              
2391             # for each character
2392 0         0 for my $i (0 .. (length($str) - 1)) {
2393            
2394             # get character
2395 0         0 $char = substr($str, $i, 1);
2396            
2397             # if a '('
2398 0 0       0 if ($char eq '(') {
    0          
2399            
2400             # if first nesting level
2401 0 0       0 if ($nest == 0) {
2402            
2403             # save index (start parameter group)
2404 0         0 $px[$j]->[0] = $i;
2405            
2406             }
2407            
2408             # increment nesting level
2409 0         0 $nest++;
2410            
2411             # if a ')'
2412             } elsif ($char eq ')') {
2413            
2414             # if second nesting level
2415 0 0       0 if ($nest == 1) {
2416            
2417             # save index (end parameter group)
2418 0         0 $px[$j]->[1] = $i;
2419            
2420             # increment parameter index
2421 0         0 $j++;
2422            
2423             }
2424            
2425             # decrement nesting level
2426 0         0 $nest--;
2427            
2428             }
2429            
2430             }
2431              
2432             # check final nesting level (should be 0)
2433 0 0       0 carp("missing ')' in parameters\n") if ($nest > 0);
2434 0 0       0 carp("missing '(' in parameters\n") if ($nest < 0);
2435              
2436             # for each parameter group
2437 0         0 for my $i (0 .. $#px) {
2438            
2439             # get parameter substring
2440 0         0 $sub = substr($str, $px[$i]->[0], $px[$i]->[1] - $px[$i]->[0] + 1);
2441            
2442             # evaluate parameter substring
2443 0         0 @eval = eval($sub);
2444            
2445             # split text within (), if 'eval' failed
2446 0 0 0     0 @eval = split(/[\s,]+/, substr($str, $px[$i]->[0] + 1, $px[$i]->[1] - $px[$i]->[0] - 1)) if ((! @eval) && ($px[$i]->[1] - $px[$i]->[0] > 1));
2447            
2448             # save parameters
2449 0         0 $pars[$i] = [@eval];
2450            
2451             }
2452              
2453             # for each parameter group (reverse order)
2454 0         0 for my $i (reverse(0 .. $#px)) {
2455            
2456             # replace parameter group with ' ($i) '
2457 0         0 substr($str, $px[$i]->[0], $px[$i]->[1] - $px[$i]->[0] + 1) = " ($i) ";
2458            
2459             }
2460              
2461             # split string into tokens, then replace '($i)' with 'eval'd parameters
2462 0 0       0 return([map {m/\((\d+)\)/ ? $pars[$1] : $_} split(/[\s,]+/, $str)]);
  0         0  
2463              
2464             }
2465              
2466             #--------- testing parameters ---------
2467              
2468             # test if a vector
2469             # elements may be numeric or 'undef'
2470             sub is_vector {
2471              
2472             # return test result
2473 0   0 0 1 0 return(ref($_[0]) eq 'ARRAY' && @{$_[0]} == grep {! defined || Scalar::Util::looks_like_number($_)} @{$_[0]});
2474              
2475             }
2476              
2477             # test if a numeric vector
2478             # elements must all be numeric
2479             sub is_num_vector {
2480              
2481             # return test result
2482 0   0 0 1 0 return(ref($_[0]) eq 'ARRAY' && @{$_[0]} == grep {Scalar::Util::looks_like_number($_)} @{$_[0]});
2483              
2484             }
2485              
2486             # test if a matrix
2487             # elements may be numeric or 'undef'
2488             sub is_matrix {
2489              
2490             # return test result
2491 0   0 0 1 0 return((ref($_[0]) eq 'ARRAY' || ref($_[0]) eq 'Math::Matrix') && @{$_[0]} == grep {is_vector($_)} @{$_[0]});
2492              
2493             }
2494              
2495             # test if a numeric matrix
2496             # elements must all be numeric
2497             sub is_num_matrix {
2498              
2499             # return test result
2500 0   0 0 1 0 return((ref($_[0]) eq 'ARRAY' || ref($_[0]) eq 'Math::Matrix') && @{$_[0]} == grep {is_num_vector($_)} @{$_[0]});
2501              
2502             }
2503              
2504             #--------- internal functions ---------
2505              
2506             # convert structure (3 parameters)
2507             # parameters: (input_array -or- input_structure, subroutine_reference)
2508             # returns: (output_array -or- output_structure)
2509             sub _convert3 {
2510              
2511             # get subroutine code reference
2512 0     0   0 my ($sub) = pop();
2513              
2514             # if valid array input (3 scalars)
2515 0 0 0     0 if (@_ == 3 && 3 == grep {! ref()} @_) {
  0 0 0     0  
      0        
2516            
2517             # call subroutine (passing @_)
2518 0         0 &$sub;
2519            
2520             # if valid structure input (array_ref -or- Math::Matrix object
2521             } elsif (@_ == 1 && (ref($_[0]) eq 'ARRAY' || UNIVERSAL::isa($_[0], 'Math::Matrix'))) {
2522            
2523             # process structure
2524 0         0 _crawl($_[0], my $out = [], $sub);
2525            
2526             # bless output array if input a Math::Matrix object
2527 0 0       0 bless($out, 'Math::Matrix') if (UNIVERSAL::isa($_[0], 'Math::Matrix'));
2528            
2529             # return output
2530 0         0 return($out);
2531            
2532             } else {
2533            
2534             # error
2535 0         0 croak('invalid input structure');
2536            
2537             }
2538            
2539             }
2540              
2541             # convert structure (4 parameters)
2542             # parameters: (input_array -or- input_structure, subroutine_reference)
2543             # returns: (output_array -or- output_structure)
2544             sub _convert4 {
2545              
2546             # get subroutine code reference
2547 0     0   0 my ($sub) = pop(@_);
2548              
2549             # get white point vector (input_structure)
2550 0 0       0 my $wtpt = (@_ == 2) ? pop() : D50;
2551            
2552             # push D50 white point vector, if none (input_array)
2553 0 0       0 push(@_, D50) if (@_ == 3);
2554              
2555             # if valid array input (3 scalars and white point vector)
2556 0 0 0     0 if (@_ == 4 && ref($_[3]) eq 'ARRAY' && 3 == grep {! ref()} @_) {
  0 0 0     0  
      0        
      0        
2557            
2558             # call subroutine (passing @_)
2559 0         0 &$sub;
2560            
2561             # if valid structure input (array_ref -or- Math::Matrix object)
2562             } elsif (@_ == 1 && (ref($_[0]) eq 'ARRAY' || UNIVERSAL::isa($_[0], 'Math::Matrix'))) {
2563            
2564             # process structure
2565 0         0 _crawl($_[0], my $out = [], $sub, $wtpt);
2566            
2567             # bless output array if input a Math::Matrix object
2568 0 0       0 bless($out, 'Math::Matrix') if (UNIVERSAL::isa($_[0], 'Math::Matrix'));
2569            
2570             # return output
2571 0         0 return($out);
2572            
2573             } else {
2574            
2575             # error
2576 0         0 croak('invalid input structure');
2577            
2578             }
2579            
2580             }
2581              
2582             # process structure
2583             # recursive subroutine
2584             # parameters: (input_structure, output_structure, subroutine_reference, [white_point_vector])
2585             sub _crawl {
2586              
2587             # get parameters
2588 0     0   0 my ($in, $out, $sub, $wtpt) = @_;
2589              
2590             # if input a reference to an array of 3 scalars
2591 0 0 0     0 if (@{$in} == 3 && 3 == grep {! ref()} @{$in}) {
  0 0       0  
  0         0  
  0         0  
2592            
2593             # call subroutine
2594 0         0 @{$out} = &$sub(@{$in}, $wtpt);
  0         0  
  0         0  
2595            
2596             # if input a reference to an array of array references
2597 0         0 } elsif (@{$in} == grep {ref() eq 'ARRAY'} @{$in}) {
  0         0  
  0         0  
2598            
2599             # for each array reference
2600 0         0 for my $i (0 .. $#{$in}) {
  0         0  
2601            
2602             # process next level
2603 0         0 _crawl($in->[$i], $out->[$i] = [], $sub, $wtpt);
2604            
2605             }
2606            
2607             } else {
2608            
2609             # error
2610 0         0 croak('invalid input structure');
2611            
2612             }
2613            
2614             }
2615              
2616             # convert xyz to L*a*b*
2617             # parameters: (x, y, z)
2618             # returns: (L, a, b)
2619             sub _xyz2Lab {
2620              
2621             # get parameters
2622 10     10   15 my ($x, $y, $z) = @_;
2623              
2624             # compute L* value
2625 10         14 my $L = x2L($y);
2626              
2627             # return L*a*b* values
2628 10         14 return($L, 500 * (x2L($x) - $L)/116, 200 * ($L - x2L($z))/116);
2629              
2630             }
2631              
2632             # convert L*a*b* to xyz
2633             # parameters: (L, a, b)
2634             # returns: (x, y, z)
2635             sub _Lab2xyz {
2636              
2637             # get parameters
2638 40     40   64 my ($L, $a, $b) = @_;
2639              
2640             # return xyz values
2641 40         123 return (L2x($L + 116 * $a/500), L2x($L), L2x($L - 116 * $b/200));
2642              
2643             }
2644              
2645             # convert xyz to LxLyLz
2646             # parameters: (x, y, z)
2647             # returns: (Lx, Ly, Lz)
2648             sub _xyz2Lxyz {
2649              
2650             # get parameters
2651 0     0     my ($x, $y, $z) = @_;
2652              
2653             # return LxLyLz values
2654 0           return(x2L($x), x2L($y), x2L($z));
2655              
2656             }
2657              
2658             # convert LxLyLz to xyz
2659             # parameters: (Lx, Ly, Lz)
2660             # returns: (x, y, z)
2661             sub _Lxyz2xyz {
2662              
2663             # get parameters
2664 0     0     my ($Lx, $Ly, $Lz) = @_;
2665              
2666             # return xyz values
2667 0           return (L2x($Lx), L2x($Ly), L2x($Lz));
2668              
2669             }
2670              
2671             # convert L*a*b* to LxLyLz
2672             # parameters: (L, a, b)
2673             # returns: (Lx, Ly, Lz)
2674             sub _Lab2Lxyz {
2675              
2676             # get parameters
2677 0     0     my ($L, $a, $b) = @_;
2678              
2679             # return LxLyLz values
2680 0           return($L + 116 * $a/500, $L, $L - 116 * $b/200);
2681              
2682             }
2683              
2684             # convert LxLyLz to L*a*b*
2685             # parameters: (Lx, Ly, Lz)
2686             # returns: (L, a, b)
2687             sub _Lxyz2Lab {
2688              
2689             # get parameters
2690 0     0     my ($Lx, $Ly, $Lz) = @_;
2691              
2692             # return L*a*b* values
2693 0           return($Ly, 500 * ($Lx - $Ly)/116, 200 * ($Ly - $Lz)/116);
2694              
2695             }
2696              
2697             # convert XYZ to L*a*b*
2698             # parameters: (X, Y, Z, ref_to_WP_vector)
2699             # returns: (L, a, b)
2700             sub _XYZ2Lab {
2701              
2702             # get parameters
2703 0     0     my ($x, $y, $z, $WPxyz) = @_;
2704              
2705             # compute L* value
2706 0           my $L = x2L($y/$WPxyz->[1]);
2707              
2708             # return L*a*b* values
2709 0           return($L, 500 * (x2L($x/$WPxyz->[0]) - $L)/116, 200 * ($L - x2L($z/$WPxyz->[2]))/116);
2710              
2711             }
2712              
2713             # convert L*a*b* to XYZ
2714             # parameters: (L, a, b, ref_to_WP_vector)
2715             # returns: (X, Y, Z)
2716             sub _Lab2XYZ {
2717              
2718             # get parameters
2719 0     0     my ($L, $a, $b, $WPxyz) = @_;
2720              
2721             # return XYZ values
2722 0           return (L2x($L + 116 * $a/500) * $WPxyz->[0], L2x($L) * $WPxyz->[1], L2x($L - 116 * $b/200) * $WPxyz->[2]);
2723              
2724             }
2725              
2726             # convert XYZ to LxLyLz
2727             # parameters: (X, Y, Z, ref_to_WP_vector)
2728             # returns: (Lx, Ly, Lz)
2729             sub _XYZ2Lxyz {
2730              
2731             # get parameters
2732 0     0     my ($x, $y, $z, $WPxyz) = @_;
2733              
2734             # return values
2735 0           return(x2L($x/$WPxyz->[0]), x2L($y/$WPxyz->[1]), x2L($z/$WPxyz->[2]));
2736              
2737             }
2738              
2739             # convert LxLyLz to XYZ
2740             # parameters: (Lx, Ly, Lz, ref_to_WP_vector)
2741             # returns: (X, Y, Z)
2742             sub _Lxyz2XYZ {
2743              
2744             # get parameters
2745 0     0     my ($Lx, $Ly, $Lz, $WPxyz) = @_;
2746              
2747             # return values
2748 0           return (L2x($Lx) * $WPxyz->[0], L2x($Ly) * $WPxyz->[1], L2x($Lz) * $WPxyz->[2]);
2749              
2750             }
2751              
2752             # convert XYZ to xyz
2753             # parameters: (X, Y, Z, ref_to_WP_vector)
2754             # returns: (x, y, z)
2755             sub _XYZ2xyz {
2756              
2757             # get parameters
2758 0     0     my ($x, $y, $z, $WPxyz) = @_;
2759              
2760             # return values
2761 0           return($x/$WPxyz->[0], $y/$WPxyz->[1], $z/$WPxyz->[2]);
2762              
2763             }
2764              
2765             # convert xyz to XYZ
2766             # parameters: (x, y, z, ref_to_WP_vector)
2767             # returns: (X, Y, Z)
2768             sub _xyz2XYZ {
2769              
2770             # get parameters
2771 0     0     my ($x, $y, $z, $WPxyz) = @_;
2772              
2773             # return values
2774 0           return($x * $WPxyz->[0], $y * $WPxyz->[1], $z * $WPxyz->[2]);
2775              
2776             }
2777              
2778             # convert XYZ to xyY
2779             # parameters: (X, Y, Z)
2780             # returns: (x, y, Y)
2781             sub _XYZ2xyY {
2782              
2783             # get parameters
2784 0     0     my ($X, $Y, $Z) = @_;
2785              
2786             # compute denominator
2787 0 0         (my $d = $X + $Y + $Z) or croak('X + Y + Z = 0 when computing chromaticity');
2788              
2789             # return values
2790 0           return($X/$d, $Y/$d, $Y);
2791              
2792             }
2793              
2794             # convert xyY to XYZ
2795             # parameters: (x, y, Y)
2796             # returns: (X, Y, Z)
2797             sub _xyY2XYZ {
2798              
2799             # get parameters
2800 0     0     my ($x, $y, $Y) = @_;
2801              
2802             # if y is zero
2803 0 0         if ($y == 0) {
2804            
2805             # error
2806 0           croak('cannot compute XYZ when y = 0');
2807            
2808             } else {
2809            
2810             # return values
2811 0           return($Y * $x/$y, $Y, $Y * (1 - $x - $y)/$y);
2812            
2813             }
2814              
2815             }
2816              
2817             # convert L*a*b* values to DIN 99 Lab values
2818             # optional parameters are Ke (lightness) and Kch (chroma)
2819             # default values of optional parameters are Ke = 1, Kch = 1
2820             # parameters: (L*a*b*, [Ke, Kch])
2821             # returns: (DIN_99_Lab)
2822             sub _Lab2DIN99 {
2823              
2824             # local variables
2825 0     0     my ($e, $f, $G, $h99, $C99, $Ke, $Kch);
2826              
2827             # compute redness
2828 0           $e = cos16 * $_[1] + sin16 * $_[2];
2829              
2830             # compute yellowness
2831 0           $f = 0.7 * (cos16 * $_[2] - sin16 * $_[1]);
2832              
2833             # compute chroma
2834 0           $G = sqrt($e**2 + $f**2);
2835              
2836             # compute DIN 99 hue angle
2837 0           $h99 = atan2($f, $e);
2838              
2839             # get Ke
2840 0 0         $Ke = defined($_[3]) ? $_[3] : 1;
2841              
2842             # get Kch
2843 0 0         $Kch = defined($_[4]) ? $_[4] : 1;
2844              
2845             # compute DIN 99 chroma
2846 0           $C99 = (log(1 + 0.045 * $G))/(0.045 * $Kch * $Ke);
2847              
2848             # return DIN 99 Lab
2849 0           return(105.509 * (log(1 + 0.0158 * $_[0])) * $Ke, $C99 * cos($h99), $C99 * sin($h99));
2850              
2851             }
2852              
2853             #--------- additional Math::Matrix methods ---------
2854              
2855             package Math::Matrix;
2856              
2857             # print object contents to string
2858             # format parameter is an array structure
2859             # parameter: ([format])
2860             # returns: (string)
2861             sub sdump {
2862              
2863             # get parameters
2864 0     0 0   my ($self, $p) = @_;
2865              
2866             # local variables
2867 0           my ($s, $fmt);
2868              
2869             # resolve parameter to a scalar
2870 0 0         $p = defined($p) ? ref($p) eq 'ARRAY' ? $p->[0] : $p : '10.5f';
    0          
2871              
2872             # set string to object ID
2873 0           $s = sprintf("'%s' object, (0x%x)\n", ref($self), Scalar::Util::refaddr($self));
2874              
2875             # append string
2876 0           $s .= "matrix values\n";
2877              
2878             # for each row
2879 0           for my $i (0 .. $#{$self}) {
  0            
2880            
2881             # make number format from parameter
2882 0           $fmt = " %$p" x @{$self->[$i]};
  0            
2883            
2884             # append matrix row
2885 0           $s .= sprintf("$fmt\n", @{$self->[$i]});
  0            
2886            
2887             }
2888              
2889             # return
2890 0           return($s);
2891              
2892             }
2893              
2894             # print object contents
2895             # format parameter is an array structure
2896             # parameter: ([format])
2897             # returns: (string)
2898             sub dump {
2899              
2900             # get parameters
2901 0     0 0   my ($self, $format) = @_;
2902              
2903             # get string from 'sdump'
2904 0           my $s = $self->sdump($format);
2905              
2906             # print string
2907 0           print $s, "\n";
2908              
2909             # return string
2910 0           return($s);
2911              
2912             }
2913              
2914             # clip matrix elements
2915             # default limits are [0, 1]
2916             # limit is disabled with value of 'none'
2917             # parameters: ([lower_limit], [upper_limit])
2918             # returns: (modified_matrix_object)
2919             sub clip {
2920              
2921             # get parameters
2922 0     0 0   my ($self, $low, $up) = @_;
2923              
2924             # set defaults, as needed
2925 0 0         $low = 0.0 if (! defined($low));
2926 0 0         $low = undef if ($low eq 'none');
2927 0 0         $up = 1.0 if (! defined($up));
2928 0 0         $up = undef if ($up eq 'none');
2929              
2930             # verify limits
2931 0 0 0       (! defined($low) || Scalar::Util::looks_like_number($low)) or croak('invalid lower clipping limit');
2932 0 0 0       (! defined($up) || Scalar::Util::looks_like_number($up)) or croak('invalid upper clipping limit');
2933 0 0 0       (! defined($low) || ! defined($up) || $up >= $low) or croak('invalid clipping limits');
      0        
2934              
2935             # for each row
2936 0           for my $i (0 .. $#{$self}) {
  0            
2937            
2938             # for each column
2939 0           for my $j (0 .. $#{$self->[$i]}) {
  0            
2940            
2941             # if matrix element is numeric
2942 0 0         if (Scalar::Util::looks_like_number($self->[$i][$j])) {
2943            
2944             # apply limits to matrix element
2945 0 0 0       $self->[$i][$j] = $low if (defined($low) && $self->[$i][$j] < $low);
2946 0 0 0       $self->[$i][$j] = $up if (defined($up) && $self->[$i][$j] > $up);
2947            
2948             }
2949            
2950             }
2951            
2952             }
2953              
2954             # return
2955 0           return($self);
2956              
2957             }
2958              
2959             # exponentiate matrix elements
2960             # exponent may be a scalar or a vector
2961             # no testing, 'inf' or 'nan' values are possible
2962             # parameter: (exponent)
2963             # returns: (new_matrix_object)
2964             sub power {
2965              
2966             # get parameters
2967 0     0 0   my ($self, $x) = @_;
2968              
2969             # local variables
2970 0           my ($result);
2971              
2972             # check if ICC::Support::Lapack module is loaded
2973 0           state $lapack = defined($INC{'ICC/Support/Lapack.pm'});
2974              
2975             # if ICC::Support::Lapack module is loaded
2976 0 0         if ($lapack) {
2977            
2978             # compute using Lapack module
2979 0           $result = ICC::Support::Lapack::power($self, $x);
2980            
2981             } else {
2982            
2983             # if exponent a scalar
2984 0 0         if (! ref($x)) {
    0          
2985            
2986             # verify exponent is a number
2987 0 0         (Scalar::Util::looks_like_number($x)) or croak('exponent must be a real number');
2988            
2989             # for each row
2990 0           for my $i (0 .. $#{$self}) {
  0            
2991            
2992             # for each column
2993 0           for my $j (0 .. $#{$self->[0]}) {
  0            
2994            
2995             # compute element
2996 0           $result->[$i][$j] = $self->[$i][$j]**$x;
2997            
2998             }
2999            
3000             }
3001            
3002             # if exponent a vector
3003             } elsif (ref($x) eq 'ARRAY') {
3004            
3005             # verify exponents are numbers
3006 0 0         (@{$x} == grep {Scalar::Util::looks_like_number($_)} @{$x}) or croak('exponents must be real numbers');
  0            
  0            
  0            
3007            
3008             # verify vector size same as column size
3009 0 0         ($#{$x} == $#{$self->[0]}) or croak('exponent vector wrong size');
  0            
  0            
3010            
3011             # for each row
3012 0           for my $i (0 .. $#{$self}) {
  0            
3013            
3014             # for each column
3015 0           for my $j (0 .. $#{$self->[0]}) {
  0            
3016            
3017             # compute element, avoiding math overflow
3018 0           $result->[$i][$j] = $self->[$i][$j]**$x->[$j];
3019            
3020             }
3021            
3022             }
3023            
3024             } else {
3025            
3026             # error
3027 0           croak('exponent must be a scalar or vector');
3028            
3029             }
3030            
3031             }
3032              
3033             # return new object
3034 0           return(bless($result, 'Math::Matrix'));
3035              
3036             }
3037              
3038             # convert matrix working space from xyz to XYZ
3039             # default XYZ white point is D50
3040             # parameter: ([XYZ_white_point])
3041             # returns: (new_matrix_object)
3042             sub xyz2XYZ {
3043              
3044             # get parameters
3045 0     0 0   my ($self, $wtpt) = @_;
3046              
3047             # local variables
3048 0           my ($mat);
3049              
3050             # verify a 3x3 matrix
3051 0 0 0       (@{$self} == 3 && @{$self->[0]} == 3) or croak('must be a 3x3 matrix');
  0            
  0            
3052              
3053             # set white point to D50 if undefined
3054 0 0         $wtpt = [96.42, 100, 82.49] if (! defined($wtpt));
3055              
3056             # clone matrix
3057 0           $mat = Storable::dclone($self);
3058              
3059             # for each row
3060 0           for my $i (0 .. 2) {
3061            
3062             # for each column
3063 0           for my $j (0 .. 2) {
3064            
3065             # modify non-diagonal element
3066 0 0         $mat->[$i][$j] *= $wtpt->[$i]/$wtpt->[$j] if ($i != $j);
3067            
3068             }
3069            
3070             }
3071              
3072             # return
3073 0           return($mat);
3074              
3075             }
3076              
3077             # convert matrix working space from XYZ to xyz
3078             # default XYZ white point is D50
3079             # parameter: ([XYZ_white_point])
3080             # returns: (new_matrix_object)
3081             sub XYZ2xyz {
3082              
3083             # get parameters
3084 0     0 0   my ($self, $wtpt) = @_;
3085              
3086             # local variables
3087 0           my ($mat);
3088              
3089             # verify a 3x3 matrix
3090 0 0 0       (@{$self} == 3 && @{$self->[0]} == 3) or croak('must be a 3x3 matrix');
  0            
  0            
3091              
3092             # set white point to D50 if undefined
3093 0 0         $wtpt = [96.42, 100, 82.49] if (! defined($wtpt));
3094              
3095             # clone matrix
3096 0           $mat = Storable::dclone($self);
3097              
3098             # for each row
3099 0           for my $i (0 .. 2) {
3100            
3101             # for each column
3102 0           for my $j (0 .. 2) {
3103            
3104             # modify non-diagonal element
3105 0 0         $mat->[$i][$j] *= $wtpt->[$j]/$wtpt->[$i] if ($i != $j);
3106            
3107             }
3108            
3109             }
3110              
3111             # return
3112 0           return($mat);
3113              
3114             }
3115              
3116             1;