line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
1
|
|
|
|
|
|
|
package ICC::Shared; |
2
|
|
|
|
|
|
|
|
3
|
41
|
|
|
41
|
|
122359
|
use strict; |
|
41
|
|
|
|
|
147
|
|
|
41
|
|
|
|
|
1128
|
|
4
|
41
|
|
|
41
|
|
182
|
use Carp; |
|
41
|
|
|
|
|
65
|
|
|
41
|
|
|
|
|
2305
|
|
5
|
41
|
|
|
41
|
|
280
|
use File::Path; |
|
41
|
|
|
|
|
86
|
|
|
41
|
|
|
|
|
1910
|
|
6
|
41
|
|
|
41
|
|
205
|
use File::Spec; |
|
41
|
|
|
|
|
62
|
|
|
41
|
|
|
|
|
936
|
|
7
|
41
|
|
|
41
|
|
191
|
use List::Util; |
|
41
|
|
|
|
|
77
|
|
|
41
|
|
|
|
|
2235
|
|
8
|
41
|
|
|
41
|
|
21564
|
use Math::Matrix; |
|
41
|
|
|
|
|
181905
|
|
|
41
|
|
|
|
|
1318
|
|
9
|
41
|
|
|
41
|
|
2808
|
use POSIX qw(:math_h :float_h); |
|
41
|
|
|
|
|
32914
|
|
|
41
|
|
|
|
|
303
|
|
10
|
41
|
|
|
41
|
|
76368
|
use Scalar::Util; |
|
41
|
|
|
|
|
97
|
|
|
41
|
|
|
|
|
1824
|
|
11
|
41
|
|
|
41
|
|
26201
|
use Storable; |
|
41
|
|
|
|
|
127951
|
|
|
41
|
|
|
|
|
2314
|
|
12
|
41
|
|
|
41
|
|
22868
|
use YAML::Tiny; |
|
41
|
|
|
|
|
245720
|
|
|
41
|
|
|
|
|
2487
|
|
13
|
41
|
|
|
41
|
|
307
|
use Exporter qw(import); |
|
41
|
|
|
|
|
75
|
|
|
41
|
|
|
|
|
1987
|
|
14
|
|
|
|
|
|
|
|
15
|
|
|
|
|
|
|
our $VERSION = 0.58; |
16
|
|
|
|
|
|
|
|
17
|
|
|
|
|
|
|
# revised 2019-11-16 |
18
|
|
|
|
|
|
|
# |
19
|
|
|
|
|
|
|
# Copyright © 2004-2019 by William B. Birkett |
20
|
|
|
|
|
|
|
|
21
|
|
|
|
|
|
|
# this module contains common methods, constants and functions |
22
|
|
|
|
|
|
|
|
23
|
|
|
|
|
|
|
# enable static variables |
24
|
41
|
|
|
41
|
|
249
|
use feature 'state'; |
|
41
|
|
|
|
|
75
|
|
|
41
|
|
|
|
|
5209
|
|
25
|
|
|
|
|
|
|
|
26
|
|
|
|
|
|
|
# constants |
27
|
41
|
|
|
41
|
|
260
|
use constant D50 => [96.42, 100, 82.49]; # ICC D50 XYZ |
|
41
|
|
|
|
|
68
|
|
|
41
|
|
|
|
|
3734
|
|
28
|
41
|
|
|
41
|
|
277
|
use constant d50 => [0.9642, 1.0, 0.8249]; # ICC D50 XYZNumber |
|
41
|
|
|
|
|
66
|
|
|
41
|
|
|
|
|
2634
|
|
29
|
41
|
|
|
41
|
|
220
|
use constant d50P => [map {$_ * 32768/65535} 0.9642, 1.0, 0.8249]; # ICC D50 16-bit XYZ PCS |
|
41
|
|
|
|
|
73
|
|
|
41
|
|
|
|
|
76
|
|
|
123
|
|
|
|
|
2805
|
|
30
|
41
|
|
|
41
|
|
234
|
use constant PI => atan2(0, -1); # pi |
|
41
|
|
|
|
|
67
|
|
|
41
|
|
|
|
|
1928
|
|
31
|
41
|
|
|
41
|
|
209
|
use constant radian => 180/PI; # radian (in degrees) |
|
41
|
|
|
|
|
66
|
|
|
41
|
|
|
|
|
2175
|
|
32
|
41
|
|
|
41
|
|
217
|
use constant ln10 => log(10); # natural log of 10 |
|
41
|
|
|
|
|
936
|
|
|
41
|
|
|
|
|
3044
|
|
33
|
41
|
|
|
41
|
|
248
|
use constant sin16 => sin(16/radian); # sin(16˚) |
|
41
|
|
|
|
|
83
|
|
|
41
|
|
|
|
|
2520
|
|
34
|
41
|
|
|
41
|
|
263
|
use constant cos16 => cos(16/radian); # cos(16˚) |
|
41
|
|
|
|
|
93
|
|
|
41
|
|
|
|
|
6704
|
|
35
|
|
|
|
|
|
|
|
36
|
|
|
|
|
|
|
BEGIN { |
37
|
|
|
|
|
|
|
|
38
|
|
|
|
|
|
|
# set constant vector elements to be read-only |
39
|
41
|
|
|
41
|
|
131
|
for (@{(D50)}) {Internals::SvREADONLY($_, 1)} |
|
41
|
|
|
|
|
112
|
|
|
123
|
|
|
|
|
255
|
|
40
|
41
|
|
|
|
|
60
|
for (@{(d50)}) {Internals::SvREADONLY($_, 1)} |
|
41
|
|
|
|
|
89
|
|
|
123
|
|
|
|
|
208
|
|
41
|
41
|
|
|
|
|
55
|
for (@{(d50P)}) {Internals::SvREADONLY($_, 1)} |
|
41
|
|
|
|
|
70
|
|
|
123
|
|
|
|
|
512029
|
|
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
|
70
|
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
|
336
|
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
|
714
|
my (@out); |
2056
|
|
|
|
|
|
|
|
2057
|
|
|
|
|
|
|
# while input array |
2058
|
612
|
|
|
|
|
1101
|
while (@_) { |
2059
|
|
|
|
|
|
|
|
2060
|
|
|
|
|
|
|
# shift input value |
2061
|
3945
|
|
|
|
|
3993
|
my $v = shift(); |
2062
|
|
|
|
|
|
|
|
2063
|
|
|
|
|
|
|
# if undefined |
2064
|
3945
|
50
|
33
|
|
|
6010
|
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
|
|
|
|
|
5889
|
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
|
|
|
|
|
303
|
push(@_, @{$v}); |
|
273
|
|
|
|
|
663
|
|
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
|
|
|
|
|
1399
|
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
|
28
|
my (@in) = @_; |
2149
|
|
|
|
|
|
|
|
2150
|
|
|
|
|
|
|
# return converted array |
2151
|
15
|
50
|
|
|
|
26
|
return(map {($_ & 0x80000000) ? $_/65536 - 65536 : $_/65536} @in); |
|
45
|
|
|
|
|
149
|
|
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
|
30
|
my (@in) = @_; |
2162
|
|
|
|
|
|
|
|
2163
|
|
|
|
|
|
|
# return converted array |
2164
|
11
|
50
|
|
|
|
18
|
return(map {($_ < 0) ? ($_ + 65536) * 65536 : $_ * 65536} @in); |
|
45
|
|
|
|
|
148
|
|
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
|
|
14
|
my ($x, $y, $z) = @_; |
2623
|
|
|
|
|
|
|
|
2624
|
|
|
|
|
|
|
# compute L* value |
2625
|
10
|
|
|
|
|
13
|
my $L = x2L($y); |
2626
|
|
|
|
|
|
|
|
2627
|
|
|
|
|
|
|
# return L*a*b* values |
2628
|
10
|
|
|
|
|
15
|
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
|
|
63
|
my ($L, $a, $b) = @_; |
2639
|
|
|
|
|
|
|
|
2640
|
|
|
|
|
|
|
# return xyz values |
2641
|
40
|
|
|
|
|
97
|
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; |