| 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 |
|
|---|
| 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 | }
|
|---|
