1 /*
2 *+
3 * Name:
4 * palUnpcd
5
6 * Purpose:
7 * Remove pincushion/barrel distortion
8
9 * Language:
10 * Starlink ANSI C
11
12 * Type of Module:
13 * Library routine
14
15 * Invocation:
16 * palUnpcd( double disco, double * x, double * y );
17
18 * Arguments:
19 * disco = double (Given)
20 * Pincushion/barrel distortion coefficient.
21 * x = double * (Given & Returned)
22 * On input the distorted X coordinate, on output
23 * the tangent-plane X coordinate.
24 * y = double * (Given & Returned)
25 * On input the distorted Y coordinate, on output
26 * the tangent-plane Y coordinate.
27
28 * Description:
29 * Remove pincushion/barrel distortion from a distorted [x,y] to give
30 * tangent-plane [x,y].
31
32 * Authors:
33 * PTW: Pat Wallace (RAL)
34 * TIMJ: Tim Jenness
35 * {enter_new_authors_here}
36
37 * Notes:
38 * - The distortion is of the form RP = R*(1+C*R^2), where R is
39 * the radial distance from the tangent point, C is the DISCO
40 * argument, and RP is the radial distance in the presence of
41 * the distortion.
42 *
43 * - For pincushion distortion, C is +ve; for barrel distortion,
44 * C is -ve.
45 *
46 * - For X,Y in "radians" - units of one projection radius,
47 * which in the case of a photograph is the focal length of
48 * the camera - the following DISCO values apply:
49 *
50 * Geometry DISCO
51 *
52 * astrograph 0.0
53 * Schmidt -0.3333
54 * AAT PF doublet +147.069
55 * AAT PF triplet +178.585
56 * AAT f/8 +21.20
57 * JKT f/8 +13.32
58 *
59 * - The present routine is a rigorous inverse of the companion
60 * routine palPcd. The expression for RP in Note 1 is rewritten
61 * in the form x^3+a*x+b=0 and solved by standard techniques.
62 *
63 * - Cases where the cubic has multiple real roots can sometimes
64 * occur, corresponding to extreme instances of barrel distortion
65 * where up to three different undistorted [X,Y]s all produce the
66 * same distorted [X,Y]. However, only one solution is returned,
67 * the one that produces the smallest change in [X,Y].
68
69 * See Also:
70 * palPcd
71
72 * History:
73 * 2000-09-03 (PTW):
74 * SLALIB implementation.
75 * 2015-01-01 (TIMJ):
76 * Initial version
77 * {enter_further_changes_here}
78
79 * Copyright:
80 * Copyright (C) 2000 Rutherford Appleton Laboratory.
81 * Copyright (C) 2015 Tim Jenness
82 * All Rights Reserved.
83
84 * Licence:
85 * This program is free software; you can redistribute it and/or
86 * modify it under the terms of the GNU General Public License as
87 * published by the Free Software Foundation; either version 3 of
88 * the License, or (at your option) any later version.
89 *
90 * This program is distributed in the hope that it will be
91 * useful, but WITHOUT ANY WARRANTY; without even the implied
92 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
93 * PURPOSE. See the GNU General Public License for more details.
94 *
95 * You should have received a copy of the GNU General Public License
96 * along with this program. If not, see <http://www.gnu.org/licenses/>.
97
98 * Bugs:
99 * {note_any_bugs_here}
100 *-
101 */
102
103 #if HAVE_CONFIG_H
104 #include <config.h>
105 #endif
106
107 #include <math.h>
108
109 #include "pal.h"
110 #include "palmac.h"
111
112 /* copysign is C99 */
113 #if HAVE_COPYSIGN
114 # define COPYSIGN copysign
115 #else
116 # define COPYSIGN(a,b) DSIGN(a,b)
117 #endif
118
palUnpcd(double disco,double * x,double * y)119 void palUnpcd( double disco, double * x, double *y ) {
120
121 const double THIRD = 1.0/3.0;
122
123 double rp,q,r,d,w,s,t,f,c,t3,f1,f2,f3,w1,w2,w3;
124 double c2;
125
126 /* Distance of the point from the origin. */
127 rp = sqrt( (*x)*(*x)+(*y)*(*y));
128
129 /* If zero, or if no distortion, no action is necessary. */
130 if (rp != 0.0 && disco != 0.0) {
131
132 /* Begin algebraic solution. */
133 q = 1.0/(3.0*disco);
134 r = rp/(2.0*disco);
135 w = q*q*q+r*r;
136
137 /* Continue if one real root, or three of which only one is positive. */
138 if (w > 0.0) {
139
140 d = sqrt(w);
141 w = r+d;
142 s = COPYSIGN(pow(fabs(w),THIRD),w);
143 w = r-d;
144 t = COPYSIGN(pow(fabs(w),THIRD),w);
145 f = s+t;
146
147 } else {
148 /* Three different real roots: use geometrical method instead. */
149 w = 2.0/sqrt(-3.0*disco);
150 c = 4.0*rp/(disco*w*w*w);
151 c2 = c*c;
152 s = sqrt(1.0-DMIN(c2,1.0));
153 t3 = atan2(s,c);
154
155 /* The three solutions. */
156 f1 = w*cos((PAL__D2PI-t3)/3.0);
157 f2 = w*cos((t3)/3.0);
158 f3 = w*cos((PAL__D2PI+t3)/3.0);
159
160 /* Pick the one that moves [X,Y] least. */
161 w1 = fabs(f1-rp);
162 w2 = fabs(f2-rp);
163 w3 = fabs(f3-rp);
164 if (w1 < w2) {
165 f = ( w1 < w3 ? f1 : f3 );
166 } else {
167 f = ( w2 < w3 ? f2 : f3 );
168 }
169 }
170
171 /* Remove the distortion. */
172 f = f/rp;
173 *x *= f;
174 *y *= f;
175 }
176 }
177