1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
|
SUBROUTINE slUPCD ( DISCO, X, Y )
*+
* - - - - - -
* U P C D
* - - - - - -
*
* Remove pincushion/barrel distortion from a distorted [x,y] to give
* tangent-plane [x,y].
*
* Given:
* DISCO d pincushion/barrel distortion coefficient
* X,Y d distorted coordinates
*
* Returned:
* X,Y d tangent-plane coordinates
*
* Notes:
*
* 1) The distortion is of the form RP = R*(1+C*R^2), where R is
* the radial distance from the tangent point, C is the DISCO
* argument, and RP is the radial distance in the presence of
* the distortion.
*
* 2) For pincushion distortion, C is +ve; for barrel distortion,
* C is -ve.
*
* 3) For X,Y in "radians" - units of one projection radius,
* which in the case of a photograph is the focal length of
* the camera - the following DISCO values apply:
*
* Geometry DISCO
*
* astrograph 0.0
* Schmidt -0.3333
* AAT PF doublet +147.069
* AAT PF triplet +178.585
* AAT f/8 +21.20
* JKT f/8 +13.32
*
* 4) The present routine is a rigorous inverse of the companion
* routine slPCD. The expression for RP in Note 1 is rewritten
* in the form x^3+a*x+b=0 and solved by standard techniques.
*
* 5) Cases where the cubic has multiple real roots can sometimes
* occur, corresponding to extreme instances of barrel distortion
* where up to three different undistorted [X,Y]s all produce the
* same distorted [X,Y]. However, only one solution is returned,
* the one that produces the smallest change in [X,Y].
*
* P.T.Wallace Starlink 3 September 2000
*
* Copyright (C) 2000 Rutherford Appleton Laboratory
*
* License:
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program (see SLA_CONDITIONS); if not, write to the
* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
* Boston, MA 02110-1301 USA
*
* Copyright (C) 1995 Association of Universities for Research in Astronomy Inc.
*-
IMPLICIT NONE
DOUBLE PRECISION DISCO,X,Y
DOUBLE PRECISION THIRD
PARAMETER (THIRD=1D0/3D0)
DOUBLE PRECISION D2PI
PARAMETER (D2PI=6.283185307179586476925286766559D0)
DOUBLE PRECISION RP,Q,R,D,W,S,T,F,C,T3,F1,F2,F3,W1,W2,W3
* Distance of the point from the origin.
RP = SQRT(X*X+Y*Y)
* If zero, or if no distortion, no action is necessary.
IF (RP.NE.0D0.AND.DISCO.NE.0D0) THEN
* Begin algebraic solution.
Q = 1D0/(3D0*DISCO)
R = RP/(2D0*DISCO)
W = Q*Q*Q+R*R
* Continue if one real root, or three of which only one is positive.
IF (W.GE.0D0) THEN
D = SQRT(W)
W = R+D
S = SIGN(ABS(W)**THIRD,W)
W = R-D
T = SIGN((ABS(W))**THIRD,W)
F = S+T
ELSE
* Three different real roots: use geometrical method instead.
W = 2D0/SQRT(-3D0*DISCO)
C = 4D0*RP/(DISCO*W*W*W)
S = SQRT(1D0-MIN(C*C,1D0))
T3 = ATAN2(S,C)
* The three solutions.
F1 = W*COS((D2PI-T3)/3D0)
F2 = W*COS((T3)/3D0)
F3 = W*COS((D2PI+T3)/3D0)
* Pick the one that moves [X,Y] least.
W1 = ABS(F1-RP)
W2 = ABS(F2-RP)
W3 = ABS(F3-RP)
IF (W1.LT.W2) THEN
IF (W1.LT.W3) THEN
F = F1
ELSE
F = F3
END IF
ELSE
IF (W2.LT.W3) THEN
F = F2
ELSE
F = F3
END IF
END IF
END IF
* Remove the distortion.
F = F/RP
X = F*X
Y = F*Y
END IF
END
|
