File Coverage

blib/lib/ICC/Support/geo2.pm
Criterion Covered Total %
statement 14 106 13.2
branch 1 48 2.0
condition 1 48 2.0
subroutine 4 10 40.0
pod 1 5 20.0
total 21 217 9.6


line stmt bran cond sub pod time code
1             package ICC::Support::geo2;
2              
3 2     2   80608 use strict;
  2         10  
  2         51  
4 2     2   10 use Carp;
  2         4  
  2         118  
5              
6             our $VERSION = 0.11;
7              
8             # revised 2016-05-17
9             #
10             # Copyright © 2004-2020 by William B. Birkett
11              
12             # inherit from Shared
13 2     2   378 use parent qw(ICC::Shared);
  2         265  
  2         10  
14              
15             # create new object
16             # hash keys are:
17             # 'points' value is an array reference containing two 3-D point coordinate arrays
18             # parameters: ([ref_to_parameter_hash])
19             # returns: (ref_to_object)
20             sub new {
21              
22             # get object class
23 1     1 0 693 my $class = shift();
24              
25             # create empty geo2 object
26 1         4 my $self = [
27             {}, # object header
28             [] # points array
29             ];
30              
31             # if one parameter, a hash reference
32 1 50 33     4 if (@_ == 1 && ref($_[0]) eq 'HASH') {
33            
34             # create new object from parameter hash
35 0         0 _new_from_hash($self, @_);
36            
37             }
38              
39             # bless object
40 1         3 bless($self, $class);
41              
42             # return object reference
43 1         2 return($self);
44              
45             }
46              
47             # compute distance and offset
48             # optional hash keys:
49             # 'limit0' limits offset to values >= 0
50             # 'limit1' limits offset to values <= 1
51             # parameters: (ref_to_input_array, [hash])
52             # returns: (distance, offset)
53             sub transform {
54              
55             # get parameters
56 0     0 0   my ($self, $in, $hash) = @_;
57              
58             # local variables
59 0           my ($v01, $v21, $vx, $s, $t);
60              
61             # compute (x0 - x1) and (x2 - x1)
62 0           $v01 = [$in->[0] - $self->[1][0][0], $in->[1] - $self->[1][0][1], $in->[2] - $self->[1][0][2]];
63 0           $v21 = [$self->[1][1][0] - $self->[1][0][0], $self->[1][1][1] - $self->[1][0][1], $self->[1][1][2] - $self->[1][0][2]];
64              
65             # compute |(x2 - x1)|^2
66 0 0         if ($s = $v21->[0]**2 + $v21->[1]**2 + $v21->[2]**2) {
67            
68             # compute offset
69 0           $t = ICC::Shared::dotProduct($v01, $v21)/$s;
70            
71             # if offset limited >= 0 and t < 0
72 0 0 0       if ($hash->{'limit0'} && $t < 0) {
    0 0        
73            
74             # return distance and offset
75 0           return(sqrt(($in->[0] - $self->[1][0][0])**2 + ($in->[1] - $self->[1][0][1])**2 + ($in->[2] - $self->[1][0][2])**2), 0);
76            
77             # if offset limited <= 0 and t > 1
78             } elsif ($hash->{'limit1'} && $t > 1) {
79            
80             # return distance and offset
81 0           return(sqrt(($in->[0] - $self->[1][1][0])**2 + ($in->[1] - $self->[1][1][1])**2 + ($in->[2] - $self->[1][1][2])**2), 1);
82            
83             } else {
84            
85             # compute (x0 - x1) x (x2 - x1)
86 0           $vx = ICC::Shared::crossProduct($v01, $v21);
87            
88             # return distance and offset
89 0           return(sqrt(($vx->[0]**2 + $vx->[1]**2 + $vx->[2]**2)/$s), $t);
90            
91             }
92            
93             # identical endpoints
94             } else {
95            
96             # return distance and offset
97 0           return(sqrt(($in->[0] - $self->[1][0][0])**2 + ($in->[1] - $self->[1][0][1])**2 + ($in->[2] - $self->[1][0][2])**2), 0);
98            
99             }
100            
101             }
102              
103             # compute Jacobian matrix
104             # optional hash keys:
105             # 'limit0' limits offset to values >= 0
106             # 'limit1' limits offset to values <= 1
107             # parameters: (ref_to_input_array, [hash])
108             # returns: (Jacobian_matrix, [distance, offset])
109             sub jacobian {
110              
111             # get parameters
112 0     0 0   my ($self, $in, $hash) = @_;
113              
114             # local variables
115 0           my ($v01, $v21, $vx, $vr, $s, $t, $d, $jac, $wx);
116              
117             # compute (x0 - x1) and (x2 - x1)
118 0           $v01 = [$in->[0] - $self->[1][0][0], $in->[1] - $self->[1][0][1], $in->[2] - $self->[1][0][2]];
119 0           $v21 = [$self->[1][1][0] - $self->[1][0][0], $self->[1][1][1] - $self->[1][0][1], $self->[1][1][2] - $self->[1][0][2]];
120              
121             # compute |(x2 - x1)|^2
122 0 0         if ($s = $v21->[0]**2 + $v21->[1]**2 + $v21->[2]**2) {
123            
124             # compute offset
125 0           $t = ICC::Shared::dotProduct($v01, $v21)/$s;
126            
127             # if offset limited >= 0 and t < 0
128 0 0 0       if ($hash->{'limit0'} && $t < 0) {
    0 0        
129            
130             # set offset
131 0           $t = 0;
132            
133             # compute Jacobian vector and radius
134 0           ($jac->[0], $d) = _radjac($self->[1][0], $in);
135            
136             # complete Jacobian
137 0           $jac->[1] = [0, 0, 0];
138            
139             # if offset limited <= 0 and t > 1
140             } elsif ($hash->{'limit1'} && $t > 1) {
141            
142             # set offset
143 0           $t = 1;
144            
145             # compute Jacobian vector and radius
146 0           ($jac->[0], $d) = _radjac($self->[1][1], $in);
147            
148             # complete Jacobian
149 0           $jac->[1] = [0, 0, 0];
150            
151             } else {
152            
153             # compute (x0 - x1) x (x2 - x1)
154 0           $vx = ICC::Shared::crossProduct($v01, $v21);
155            
156             # compute distance
157 0           $d = sqrt(($vx->[0]**2 + $vx->[1]**2 + $vx->[2]**2)/$s);
158            
159             # compute offset partial derivatives
160 0           $jac->[1] = [$v21->[0]/$s, $v21->[1]/$s, $v21->[2]/$s];
161            
162             # compute cross product matrix
163 0           $wx = [
164             [0, $jac->[1][2]/$d, -$jac->[1][1]/$d],
165             [-$jac->[1][2]/$d, 0, $jac->[1][0]/$d],
166             [$jac->[1][1]/$d, -$jac->[1][0]/$d, 0]
167             ];
168            
169             # compute distance partial derivatives
170 0           $jac->[0] = ICC::Support::Lapack::vec_xplus($wx, $vx, {'trans' => 'T'});
171            
172             }
173            
174             # identical endpoints
175             } else {
176            
177             # set offset
178 0           $t = 0;
179            
180             # compute Jacobian vector and radius
181 0           ($jac->[0], $d) = _radjac($self->[1][0], $in);
182            
183             # complete Jacobian
184 0           $jac->[1] = [0, 0, 0];
185            
186             }
187            
188             # bless Jacobian as Math::Matrix object
189 0           bless($jac, 'Math::Matrix');
190              
191             # if array wanted
192 0 0         if (wantarray) {
193            
194             # return Jacobian, distance and offset
195 0           return($jac, $d, $t);
196            
197             } else {
198            
199             # return Jacobian
200 0           return($jac);
201            
202             }
203            
204             }
205              
206             # get/set points array
207             # parameters: ([ref_to_points_array])
208             # returns: (ref_to_points_array)
209             sub points {
210              
211             # get parameters
212 0     0 0   my ($self, $points) = @_;
213              
214             # if parameter supplied
215 0 0         if (defined($points)) {
216            
217             # verify a 2-D array
218 0 0 0       (ref($points) eq 'ARRAY' && @{$points} == grep {ref() eq 'ARRAY'} @{$points}) or croak('\'points\' parameter not a 2-D array');
  0            
  0            
  0            
219            
220             # verify array has 2 rows
221 0 0         (@{$points} == 2) or croak('\'points\' parameter must contain 2 points');
  0            
222            
223             # verify point 0 contains 3 coordinates
224 0 0 0       (@{$points->[0]} == 3 && 3 == grep {Scalar::Util::looks_like_number($_)} @{$points->[0]}) or croak('\'points\' parameter has invalid point 0');
  0            
  0            
  0            
225            
226             # verify point 1 contains 3 coordinates
227 0 0 0       (@{$points->[1]} == 3 && 3 == grep {Scalar::Util::looks_like_number($_)} @{$points->[1]}) or croak('\'points\' parameter has invalid point 1');
  0            
  0            
  0            
228            
229             # verify points are unique
230 0 0 0       ($points->[0][0] != $points->[1][0] || $points->[0][1] != $points->[1][1] || $points->[0][2] != $points->[1][2]) or carp('\'points\' parameter contains identical points');
      0        
231            
232             # copy points array
233 0           $self->[1] = Storable::dclone($points);
234            
235             }
236              
237             # return end point array reference
238 0           return($self->[1]);
239              
240             }
241              
242             # print object contents to string
243             # format is an array structure
244             # parameter: ([format])
245             # returns: (string)
246             sub sdump {
247              
248             # get parameters
249 0     0 1   my ($self, $p) = @_;
250              
251             # local variables
252 0           my ($s, $fmt);
253              
254             # resolve parameter to an array reference
255 0 0         $p = defined($p) ? ref($p) eq 'ARRAY' ? $p : [$p] : [];
    0          
256              
257             # get format string
258 0 0 0       $fmt = defined($p->[0]) && ! ref($p->[0]) ? $p->[0] : 'undef';
259              
260             # set string to object ID
261 0           $s = sprintf("'%s' object, (0x%x)\n", ref($self), $self);
262              
263             # return
264 0           return($s);
265              
266             }
267              
268             # compute radial Jacobian
269             # parameters: (point_vector, ref_to_input_vector)
270             # returns: (Jacobian_vector, radius)
271             sub _radjac {
272              
273             # get parameters
274 0     0     my ($point, $in) = @_;
275              
276             # local variables
277 0           my ($jac, $r);
278              
279             # for each dimension
280 0           for my $i (0 .. 2) {
281            
282             # compute Jacobian element
283 0           $jac->[$i] = ($in->[$i] - $point->[$i]);
284            
285             }
286              
287             # compute radius
288 0           $r = sqrt($jac->[0]**2 + $jac->[1]**2 + $jac->[2]**2);
289              
290             # if radius is zero
291 0 0         if ($r == 0) {
292            
293             # set Jacobian to all ones
294 0           $jac = [1, 1, 1];
295            
296             } else {
297            
298             # for each dimension
299 0           for my $i (0 .. 2) {
300            
301             # divide by radius
302 0           $jac->[$i] /= $r;
303            
304             }
305            
306             }
307              
308             # return
309 0           return($jac, $r);
310              
311             }
312              
313             # populate object from parameter hash
314             # parameters: (object_reference, ref_to_parameter_hash)
315             sub _new_from_hash {
316              
317             # get parameters
318 0     0     my ($self, $hash) = @_;
319              
320             # local variables
321 0           my ($points);
322              
323             # get points array
324 0 0         if ($points = $hash->{'points'}) {
325            
326             # verify a 2-D array
327 0 0 0       (ref($points) eq 'ARRAY' && @{$points} == grep {ref() eq 'ARRAY'} @{$points}) or croak('\'points\' parameter not a 2-D array');
  0            
  0            
  0            
328            
329             # verify array has 2 rows
330 0 0         (@{$points} == 2) or croak('\'points\' parameter must contain 2 points');
  0            
331            
332             # verify point 0 contains 3 coordinates
333 0 0 0       (@{$points->[0]} == 3 && 3 == grep {Scalar::Util::looks_like_number($_)} @{$points->[0]}) or croak('\'points\' parameter has invalid point 0');
  0            
  0            
  0            
334            
335             # verify point 1 contains 3 coordinates
336 0 0 0       (@{$points->[1]} == 3 && 3 == grep {Scalar::Util::looks_like_number($_)} @{$points->[1]}) or croak('\'points\' parameter has invalid point 1');
  0            
  0            
  0            
337            
338             # verify points are unique
339 0 0 0       ($points->[0][0] != $points->[1][0] || $points->[0][1] != $points->[1][1] || $points->[0][2] != $points->[1][2]) or carp('\'points\' parameter contains identical points');
      0        
340            
341             # copy points array
342 0           $self->[1] = Storable::dclone($points);
343            
344             }
345              
346             }
347              
348             1;