wcslib (8.2.2)
1 /*============================================================================
2 WCSLIB 8.2 - an implementation of the FITS WCS standard.
3 Copyright (C) 1995-2023, Mark Calabretta
4
5 This file is part of WCSLIB.
6
7 WCSLIB is free software: you can redistribute it and/or modify it under the
8 terms of the GNU Lesser General Public License as published by the Free
9 Software Foundation, either version 3 of the License, or (at your option)
10 any later version.
11
12 WCSLIB is distributed in the hope that it will be useful, but WITHOUT ANY
13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15 more details.
16
17 You should have received a copy of the GNU Lesser General Public License
18 along with WCSLIB. If not, see http://www.gnu.org/licenses.
19
20 Author: Mark Calabretta, Australia Telescope National Facility, CSIRO.
21 http://www.atnf.csiro.au/people/Mark.Calabretta
22 $Id: dis.h,v 8.2.1.1 2023/11/16 10:05:57 mcalabre Exp mcalabre $
23 *=============================================================================
24 *
25 * WCSLIB 8.2 - C routines that implement the FITS World Coordinate System
26 * (WCS) standard. Refer to the README file provided with WCSLIB for an
27 * overview of the library.
28 *
29 *
30 * Summary of the dis routines
31 * ---------------------------
32 * Routines in this suite implement extensions to the FITS World Coordinate
33 * System (WCS) standard proposed by
34 *
35 = "Representations of distortions in FITS world coordinate systems",
36 = Calabretta, M.R. et al. (WCS Paper IV, draft dated 2004/04/22),
37 = available from http://www.atnf.csiro.au/people/Mark.Calabretta
38 *
39 * In brief, a distortion function may occupy one of two positions in the WCS
40 * algorithm chain. Prior distortions precede the linear transformation
41 * matrix, whether it be PCi_ja or CDi_ja, and sequent distortions follow it.
42 * WCS Paper IV defines FITS keywords used to specify parameters for predefined
43 * distortion functions. The following are used for prior distortions:
44 *
45 = CPDISja ...(string-valued, identifies the distortion function)
46 = DPja ...(record-valued, parameters)
47 = CPERRja ...(floating-valued, maximum value)
48 *
49 * Their counterparts for sequent distortions are CQDISia, DQia, and CQERRia.
50 * An additional floating-valued keyword, DVERRa, records the maximum value of
51 * the combined distortions.
52 *
53 * DPja and DQia are "record-valued". Syntactically, the keyvalues are
54 * standard FITS strings, but they are to be interpreted in a special way.
55 * The general form is
56 *
57 = DPja = '<field-specifier>: <float>'
58 *
59 * where the field-specifier consists of a sequence of fields separated by
60 * periods, and the ': ' between the field-specifier and the floating-point
61 * value is part of the record syntax. For example:
62 *
63 = DP1 = 'AXIS.1: 1'
64 *
65 * Certain field-specifiers are defined for all distortion functions, while
66 * others are defined only for particular distortions. Refer to WCS Paper IV
67 * for further details. wcspih() parses all distortion keywords and loads them
68 * into a disprm struct for analysis by disset() which knows (or possibly does
69 * not know) how to interpret them. Of the Paper IV distortion functions, only
70 * the general Polynomial distortion is currently implemented here.
71 *
72 * TPV - the TPV "projection":
73 * ---------------------------
74 * The distortion function component of the TPV celestial "projection" is also
75 * supported. The TPV projection, originally proposed in a draft of WCS Paper
76 * II, consists of a TAN projection with sequent polynomial distortion, the
77 * coefficients of which are encoded in PVi_ma keyrecords. Full details may be
78 * found at the registry of FITS conventions:
79 *
80 = http://fits.gsfc.nasa.gov/registry/tpvwcs/tpv.html
81 *
82 * Internally, wcsset() changes TPV to a TAN projection, translates the PVi_ma
83 * keywords to DQia and loads them into a disprm struct. These DQia keyrecords
84 * have the form
85 *
86 = DQia = 'TPV.m: <value>'
87 *
88 * where i, a, m, and the value for each DQia match each PVi_ma. Consequently,
89 * WCSLIB would handle a FITS header containing these keywords, along with
90 * CQDISia = 'TPV' and the required DQia.NAXES and DQia.AXIS.ihat keywords.
91 *
92 * Note that, as defined, TPV assumes that CDi_ja is used to define the linear
93 * transformation. The section on historical idiosyncrasies (below) cautions
94 * about translating CDi_ja to PCi_ja plus CDELTia in this case.
95 *
96 * SIP - Simple Imaging Polynomial:
97 * --------------------------------
98 * These routines also support the Simple Imaging Polynomial (SIP), whose
99 * design was influenced by early drafts of WCS Paper IV. It is described in
100 * detail in
101 *
102 = http://fits.gsfc.nasa.gov/registry/sip.html
103 *
104 * SIP, which is defined only as a prior distortion for 2-D celestial images,
105 * has the interesting feature that it records an approximation to the inverse
106 * polynomial distortion function. This is used by disx2p() to provide an
107 * initial estimate for its more precise iterative inversion. The
108 * special-purpose keywords used by SIP are parsed and translated by wcspih()
109 * as follows:
110 *
111 = A_p_q = <value> -> DP1 = 'SIP.FWD.p_q: <value>'
112 = AP_p_q = <value> -> DP1 = 'SIP.REV.p_q: <value>'
113 = B_p_q = <value> -> DP2 = 'SIP.FWD.p_q: <value>'
114 = BP_p_q = <value> -> DP2 = 'SIP.REV.p_q: <value>'
115 = A_DMAX = <value> -> DPERR1 = <value>
116 = B_DMAX = <value> -> DPERR2 = <value>
117 *
118 * SIP's A_ORDER and B_ORDER keywords are not used. WCSLIB would recognise a
119 * FITS header containing the above keywords, along with CPDISja = 'SIP' and
120 * the required DPja.NAXES keywords.
121 *
122 * DSS - Digitized Sky Survey:
123 * ---------------------------
124 * The Digitized Sky Survey resulted from the production of the Guide Star
125 * Catalogue for the Hubble Space Telescope. Plate solutions based on a
126 * polynomial distortion function were encoded in FITS using non-standard
127 * keywords. Sect. 5.2 of WCS Paper IV describes how DSS coordinates may be
128 * translated to a sequent Polynomial distortion using two auxiliary variables.
129 * That translation is based on optimising the non-distortion component of the
130 * plate solution.
131 *
132 * Following Paper IV, wcspih() translates the non-distortion component of DSS
133 * coordinates to standard WCS keywords (CRPIXja, PCi_ja, CRVALia, etc), and
134 * fills a wcsprm struct with their values. It encodes the DSS polynomial
135 * coefficients as
136 *
137 = AMDXm = <value> -> DQ1 = 'AMD.m: <value>'
138 = AMDYm = <value> -> DQ2 = 'AMD.m: <value>'
139 *
140 * WCSLIB would recognise a FITS header containing the above keywords, along
141 * with CQDISia = 'DSS' and the required DQia.NAXES keywords.
142 *
143 * WAT - the TNX and ZPX "projections":
144 * ------------------------------------
145 * The TNX and ZPX "projections" add a polynomial distortion function to the
146 * standard TAN and ZPN projections respectively. Unusually, the polynomial
147 * may be expressed as the sum of Chebyshev or Legendre polynomials, or as a
148 * simple sum of monomials, as described in
149 *
150 = http://fits.gsfc.nasa.gov/registry/tnx/tnx-doc.html
151 = http://fits.gsfc.nasa.gov/registry/zpxwcs/zpx.html
152 *
153 * The polynomial coefficients are encoded in special-purpose WATi_n keywords
154 * as a set of continued strings, thus providing the name for this distortion
155 * type. WATi_n are parsed and translated by wcspih() into the following set:
156 *
157 = DQi = 'WAT.POLY: <value>'
158 = DQi = 'WAT.XMIN: <value>'
159 = DQi = 'WAT.XMAX: <value>'
160 = DQi = 'WAT.YMIN: <value>'
161 = DQi = 'WAT.YMAX: <value>'
162 = DQi = 'WAT.CHBY.m_n: <value>' or
163 = DQi = 'WAT.LEGR.m_n: <value>' or
164 = DQi = 'WAT.MONO.m_n: <value>'
165 *
166 * along with CQDISia = 'WAT' and the required DPja.NAXES keywords. For ZPX,
167 * the ZPN projection parameters are also encoded in WATi_n, and wcspih()
168 * translates these to standard PVi_ma.
169 *
170 * Note that, as defined, TNX and ZPX assume that CDi_ja is used to define the
171 * linear transformation. The section on historical idiosyncrasies (below)
172 * cautions about translating CDi_ja to PCi_ja plus CDELTia in this case.
173 *
174 * TPD - Template Polynomial Distortion:
175 * -------------------------------------
176 * The "Template Polynomial Distortion" (TPD) is a superset of the TPV, SIP,
177 * DSS, and WAT (TNX & ZPX) polynomial distortions that also supports 1-D usage
178 * and inversions. Like TPV, SIP, and DSS, the form of the polynomial is fixed
179 * (the "template") and only the coefficients for the required terms are set
180 * non-zero. TPD generalizes TPV in going to 9th degree, SIP by accomodating
181 * TPV's linear and radial terms, and DSS in both respects. While in theory
182 * the degree of the WAT polynomial distortion in unconstrained, in practice it
183 * is limited to values that can be handled by TPD.
184 *
185 * Within WCSLIB, TPV, SIP, DSS, and WAT are all implemented as special cases
186 * of TPD. Indeed, TPD was developed precisely for that purpose. WAT
187 * distortions expressed as the sum of Chebyshev or Legendre polynomials are
188 * expanded for TPD as a simple sum of monomials. Moreover, the general
189 * Polynomial distortion is translated and implemented internally as TPD
190 * whenever possible.
191 *
192 * However, WCSLIB also recognizes 'TPD' as a distortion function in its own
193 * right (i.e. a recognized value of CPDISja or CQDISia), for use as both prior
194 * and sequent distortions. Its DPja and DQia keyrecords have the form
195 *
196 = DPja = 'TPD.FWD.m: <value>'
197 = DPja = 'TPD.REV.m: <value>'
198 *
199 * for the forward and reverse distortion functions. Moreover, like the
200 * general Polynomial distortion, TPD supports auxiliary variables, though only
201 * as a linear transformation of pixel coordinates (p1,p2):
202 *
203 = x = a0 + a1*p1 + a2*p2
204 = y = b0 + b1*p1 + b2*p2
205 *
206 * where the coefficients of the auxiliary variables (x,y) are recorded as
207 *
208 = DPja = 'AUX.1.COEFF.0: a0' ...default 0.0
209 = DPja = 'AUX.1.COEFF.1: a1' ...default 1.0
210 = DPja = 'AUX.1.COEFF.2: a2' ...default 0.0
211 = DPja = 'AUX.2.COEFF.0: b0' ...default 0.0
212 = DPja = 'AUX.2.COEFF.1: b1' ...default 0.0
213 = DPja = 'AUX.2.COEFF.2: b2' ...default 1.0
214 *
215 * Though nowhere near as powerful, in typical applications TPD is considerably
216 * faster than the general Polynomial distortion. As TPD has a finite and not
217 * too large number of possible terms (60), the coefficients for each can be
218 * stored (by disset()) in a fixed location in the disprm::dparm[] array. A
219 * large part of the speedup then arises from evaluating the polynomial using
220 * Horner's scheme.
221 *
222 * Separate implementations for polynomials of each degree, and conditionals
223 * for 1-D polynomials and 2-D polynomials with and without the radial
224 * variable, ensure that unused terms mostly do not impose a significant
225 * computational overhead.
226 *
227 * The TPD terms are as follows
228 *
229 = 0: 1 4: xx 12: xxxx 24: xxxxxx 40: xxxxxxxx
230 = 5: xy 13: xxxy 25: xxxxxy 41: xxxxxxxy
231 = 1: x 6: yy 14: xxyy 26: xxxxyy 42: xxxxxxyy
232 = 2: y 15: xyyy 27: xxxyyy 43: xxxxxyyy
233 = 3: r 7: xxx 16: yyyy 28: xxyyyy 44: xxxxyyyy
234 = 8: xxy 29: xyyyyy 45: xxxyyyyy
235 = 9: xyy 17: xxxxx 30: yyyyyy 46: xxyyyyyy
236 = 10: yyy 18: xxxxy 47: xyyyyyyy
237 = 11: rrr 19: xxxyy 31: xxxxxxx 48: yyyyyyyy
238 = 20: xxyyy 32: xxxxxxy
239 = 21: xyyyy 33: xxxxxyy 49: xxxxxxxxx
240 = 22: yyyyy 34: xxxxyyy 50: xxxxxxxxy
241 = 23: rrrrr 35: xxxyyyy 51: xxxxxxxyy
242 = 36: xxyyyyy 52: xxxxxxyyy
243 = 37: xyyyyyy 53: xxxxxyyyy
244 = 38: yyyyyyy 54: xxxxyyyyy
245 = 39: rrrrrrr 55: xxxyyyyyy
246 = 56: xxyyyyyyy
247 = 57: xyyyyyyyy
248 = 58: yyyyyyyyy
249 = 59: rrrrrrrrr
250 *
251 * where r = sqrt(xx + yy). Note that even powers of r are excluded since they
252 * can be accomodated by powers of (xx + yy).
253 *
254 * Note here that "x" refers to the axis to which the distortion function is
255 * attached, with "y" being the complementary axis. So, for example, with
256 * longitude on axis 1 and latitude on axis 2, for TPD attached to axis 1, "x"
257 * refers to axis 1 and "y" to axis 2. For TPD attached to axis 2, "x" refers
258 * to axis 2, and "y" to axis 1.
259 *
260 * TPV uses all terms up to 39. The m in its PVi_ma keywords translates
261 * directly to the TPD coefficient number.
262 *
263 * SIP uses all terms except for 0, 3, 11, 23, 39, and 59, with terms 1 and 2
264 * only used for the inverse. Its A_p_q, etc. keywords must be translated
265 * using a map.
266 *
267 * DSS uses terms 0, 1, 2, 4, 5, 6, 7, 8, 9, 10, 17, 19, and 21. The presence
268 * of a non-zero constant term arises through the use of auxiliary variables
269 * with origin offset from the reference point of the TAN projection. However,
270 * in the translation given by WCS Paper IV, the distortion polynomial is zero,
271 * or very close to zero, at the reference pixel itself. The mapping between
272 * DSS's AMDXm (or AMDYm) keyvalues and TPD coefficients, while still simple,
273 * is not quite as straightforward as for TPV and SIP.
274 *
275 * WAT uses all but the radial terms, namely 3, 11, 23, 39, and 59. While the
276 * mapping between WAT's monomial coefficients and TPD is fairly simple, for
277 * its expression in terms of a sum of Chebyshev or Legendre polynomials it is
278 * much less so.
279 *
280 * Historical idiosyncrasies:
281 * --------------------------
282 * In addition to the above, some historical distortion functions have further
283 * idiosyncrasies that must be taken into account when translating them to TPD.
284 *
285 * WCS Paper IV specifies that a distortion function returns a correction to be
286 * added to pixel coordinates (prior distortion) or intermediate pixel
287 * coordinates (sequent distortion). The correction is meant to be small so
288 * that ignoring the distortion function, i.e. setting the correction to zero,
289 * produces a commensurately small error.
290 *
291 * However, rather than an additive correction, some historical distortion
292 * functions (TPV, DSS) define a polynomial that returns the corrected
293 * coordinates directly.
294 *
295 * The difference between the two approaches is readily accounted for simply by
296 * adding or subtracting 1 from the coefficient of the first degree term of the
297 * polynomial. However, it opens the way for considerable confusion.
298 *
299 * Additional to the formalism of WCS Paper IV, both the Polynomial and TPD
300 * distortion functions recognise a keyword
301 *
302 = DPja = 'DOCORR: 0'
303 *
304 * which is meant to apply generally to indicate that the distortion function
305 * returns the corrected coordinates directly. Any other value for DOCORR (or
306 * its absence) indicates that the distortion function returns an additive
307 * correction.
308 *
309 * WCS Paper IV also specifies that the independent variables of a distortion
310 * function are pixel coordinates (prior distortion) or intermediate pixel
311 * coordinates (sequent distortion).
312 *
313 * On the contrary, the independent variables of the SIP polynomial are pixel
314 * coordinate offsets from the reference pixel. This is readily handled via
315 * the renormalisation parameters
316 *
317 = DPja = 'OFFSET.jhat: <value>'
318 *
319 * where the value corresponds to CRPIXja.
320 *
321 * Likewise, because TPV, TNX, and ZPX are defined in terms of CDi_ja, the
322 * independent variables of the polynomial are intermediate world coordinates
323 * rather than intermediate pixel coordinates. Because sequent distortions
324 * are always applied before CDELTia, if CDi_ja is translated to PCi_ja plus
325 * CDELTia, then either CDELTia must be unity, or the distortion polynomial
326 * coefficients must be adjusted to account for the change of scale.
327 *
328 * Summary of the dis routines:
329 * ----------------------------
330 * These routines apply the distortion functions defined by the extension to
331 * the FITS WCS standard proposed in Paper IV. They are based on the disprm
332 * struct which contains all information needed for the computations. The
333 * struct contains some members that must be set by the user, and others that
334 * are maintained by these routines, somewhat like a C++ class but with no
335 * encapsulation.
336 *
337 * dpfill(), dpkeyi(), and dpkeyd() are provided to manage the dpkey struct.
338 *
339 * disndp(), disini(), disinit(), discpy(), and disfree() are provided to
340 * manage the disprm struct, dissize() computes its total size including
341 * allocated memory, and disprt() prints its contents.
342 *
343 * disperr() prints the error message(s) (if any) stored in a disprm struct.
344 *
345 * wcshdo() normally writes SIP and TPV headers in their native form if at all
346 * possible. However, dishdo() may be used to set a flag that tells it to
347 * write the header in the form of the TPD translation used internally.
348 *
349 * A setup routine, disset(), computes intermediate values in the disprm struct
350 * from parameters in it that were supplied by the user. The struct always
351 * needs to be set up by disset(), though disset() need not be called
352 * explicitly - refer to the explanation of disprm::flag.
353 *
354 * disp2x() and disx2p() implement the WCS distortion functions, disp2x() using
355 * separate functions, such as dispoly() and tpd7(), to do the computation.
356 *
357 * An auxiliary routine, diswarp(), computes various measures of the distortion
358 * over a specified range of coordinates.
359 *
360 * PLEASE NOTE: Distortions are not yet handled by wcsbth(), or wcscompare().
361 *
362 *
363 * disndp() - Memory allocation for DPja and DQia
364 * ----------------------------------------------
365 * disndp() sets or gets the value of NDPMAX (default 256). This global
366 * variable controls the maximum number of dpkey structs, for holding DPja or
367 * DQia keyvalues, that disini() should allocate space for. It is also used by
368 * disinit() as the default value of ndpmax.
369 *
370 * PLEASE NOTE: This function is not thread-safe.
371 *
372 * Given:
373 * n int Value of NDPMAX; ignored if < 0. Use a value less
374 * than zero to get the current value.
375 *
376 * Function return value:
377 * int Current value of NDPMAX.
378 *
379 *
380 * dpfill() - Fill the contents of a dpkey struct
381 * ----------------------------------------------
382 * dpfill() is a utility routine to aid in filling the contents of the dpkey
383 * struct. No checks are done on the validity of the inputs.
384 *
385 * WCS Paper IV specifies the syntax of a record-valued keyword as
386 *
387 = keyword = '<field-specifier>: <float>'
388 *
389 * However, some DPja and DQia record values, such as those of DPja.NAXES and
390 * DPja.AXIS.j, are intrinsically integer-valued. While FITS header parsers
391 * are not expected to know in advance which of DPja and DQia are integral and
392 * which are floating point, if the record's value parses as an integer (i.e.
393 * without decimal point or exponent), then preferably enter it into the dpkey
394 * struct as an integer. Either way, it doesn't matter as disset() accepts
395 * either data type for all record values.
396 *
397 * Given and returned:
398 * dp struct dpkey*
399 * Store for DPja and DQia keyvalues.
400 *
401 * Given:
402 * keyword const char *
403 * field const char *
404 * These arguments are concatenated with an intervening
405 * "." to construct the full record field name, i.e.
406 * including the keyword name, DPja or DQia (but
407 * excluding the colon delimiter which is NOT part of the
408 * name). Either may be given as a NULL pointer. Set
409 * both NULL to omit setting this component of the
410 * struct.
411 *
412 * j int Axis number (1-relative), i.e. the j in DPja or
413 * i in DQia. Can be given as 0, in which case the axis
414 * number will be obtained from the keyword component of
415 * the field name which must either have been given or
416 * preset.
417 *
418 * If j is non-zero, and keyword was given, then the
419 * value of j will be used to fill in the axis number.
420 *
421 * type int Data type of the record's value
422 * 0: Integer,
423 * 1: Floating point.
424 *
425 * i int For type == 0, the integer value of the record.
426 *
427 * f double For type == 1, the floating point value of the record.
428 *
429 * Function return value:
430 * int Status return value:
431 * 0: Success.
432 *
433 *
434 * dpkeyi() - Get the data value in a dpkey struct as int
435 * ------------------------------------------------------
436 * dpkeyi() returns the data value in a dpkey struct as an integer value.
437 *
438 * Given and returned:
439 * dp const struct dpkey *
440 * Parsed contents of a DPja or DQia keyrecord.
441 *
442 * Function return value:
443 * int The record's value as int.
444 *
445 *
446 * dpkeyd() - Get the data value in a dpkey struct as double
447 * ---------------------------------------------------------
448 * dpkeyd() returns the data value in a dpkey struct as a floating point
449 * value.
450 *
451 * Given and returned:
452 * dp const struct dpkey *
453 * Parsed contents of a DPja or DQia keyrecord.
454 *
455 * Function return value:
456 * double The record's value as double.
457 *
458 *
459 * disini() - Default constructor for the disprm struct
460 * ----------------------------------------------------
461 * disini() is a thin wrapper on disinit(). It invokes it with ndpmax set
462 * to -1 which causes it to use the value of the global variable NDPMAX. It
463 * is thereby potentially thread-unsafe if NDPMAX is altered dynamically via
464 * disndp(). Use disinit() for a thread-safe alternative in this case.
465 *
466 *
467 * disinit() - Default constructor for the disprm struct
468 * ----------------------------------------------------
469 * disinit() allocates memory for arrays in a disprm struct and sets all
470 * members of the struct to default values.
471 *
472 * PLEASE NOTE: every disprm struct must be initialized by disinit(), possibly
473 * repeatedly. On the first invokation, and only the first invokation,
474 * disprm::flag must be set to -1 to initialize memory management, regardless
475 * of whether disinit() will actually be used to allocate memory.
476 *
477 * Given:
478 * alloc int If true, allocate memory unconditionally for arrays in
479 * the disprm struct.
480 *
481 * If false, it is assumed that pointers to these arrays
482 * have been set by the user except if they are null
483 * pointers in which case memory will be allocated for
484 * them regardless. (In other words, setting alloc true
485 * saves having to initalize these pointers to zero.)
486 *
487 * naxis int The number of world coordinate axes, used to determine
488 * array sizes.
489 *
490 * Given and returned:
491 * dis struct disprm*
492 * Distortion function parameters. Note that, in order
493 * to initialize memory management disprm::flag must be
494 * set to -1 when dis is initialized for the first time
495 * (memory leaks may result if it had already been
496 * initialized).
497 *
498 * Given:
499 * ndpmax int The number of DPja or DQia keywords to allocate space
500 * for. If set to -1, the value of the global variable
501 * NDPMAX will be used. This is potentially
502 * thread-unsafe if disndp() is being used dynamically to
503 * alter its value.
504 *
505 * Function return value:
506 * int Status return value:
507 * 0: Success.
508 * 1: Null disprm pointer passed.
509 * 2: Memory allocation failed.
510 *
511 * For returns > 1, a detailed error message is set in
512 * disprm::err if enabled, see wcserr_enable().
513 *
514 *
515 * discpy() - Copy routine for the disprm struct
516 * ---------------------------------------------
517 * discpy() does a deep copy of one disprm struct to another, using disinit()
518 * to allocate memory unconditionally for its arrays if required. Only the
519 * "information to be provided" part of the struct is copied; a call to
520 * disset() is required to initialize the remainder.
521 *
522 * Given:
523 * alloc int If true, allocate memory unconditionally for arrays in
524 * the destination. Otherwise, it is assumed that
525 * pointers to these arrays have been set by the user
526 * except if they are null pointers in which case memory
527 * will be allocated for them regardless.
528 *
529 * dissrc const struct disprm*
530 * Struct to copy from.
531 *
532 * Given and returned:
533 * disdst struct disprm*
534 * Struct to copy to. disprm::flag should be set to -1
535 * if disdst was not previously initialized (memory leaks
536 * may result if it was previously initialized).
537 *
538 * Function return value:
539 * int Status return value:
540 * 0: Success.
541 * 1: Null disprm pointer passed.
542 * 2: Memory allocation failed.
543 *
544 * For returns > 1, a detailed error message is set in
545 * disprm::err if enabled, see wcserr_enable().
546 *
547 *
548 * disfree() - Destructor for the disprm struct
549 * --------------------------------------------
550 * disfree() frees memory allocated for the disprm arrays by disinit().
551 * disinit() keeps a record of the memory it allocates and disfree() will only
552 * attempt to free this.
553 *
554 * PLEASE NOTE: disfree() must not be invoked on a disprm struct that was not
555 * initialized by disinit().
556 *
557 * Given:
558 * dis struct disprm*
559 * Distortion function parameters.
560 *
561 * Function return value:
562 * int Status return value:
563 * 0: Success.
564 * 1: Null disprm pointer passed.
565 *
566 *
567 * dissize() - Compute the size of a disprm struct
568 * -----------------------------------------------
569 * dissize() computes the full size of a disprm struct, including allocated
570 * memory.
571 *
572 * Given:
573 * dis const struct disprm*
574 * Distortion function parameters.
575 *
576 * If NULL, the base size of the struct and the allocated
577 * size are both set to zero.
578 *
579 * Returned:
580 * sizes int[2] The first element is the base size of the struct as
581 * returned by sizeof(struct disprm). The second element
582 * is the total allocated size, in bytes, assuming that
583 * the allocation was done by disini(). This figure
584 * includes memory allocated for members of constituent
585 * structs, such as disprm::dp.
586 *
587 * It is not an error for the struct not to have been set
588 * up via tabset(), which normally results in additional
589 * memory allocation.
590 *
591 * Function return value:
592 * int Status return value:
593 * 0: Success.
594 *
595 *
596 * disprt() - Print routine for the disprm struct
597 * ----------------------------------------------
598 * disprt() prints the contents of a disprm struct using wcsprintf(). Mainly
599 * intended for diagnostic purposes.
600 *
601 * Given:
602 * dis const struct disprm*
603 * Distortion function parameters.
604 *
605 * Function return value:
606 * int Status return value:
607 * 0: Success.
608 * 1: Null disprm pointer passed.
609 *
610 *
611 * disperr() - Print error messages from a disprm struct
612 * -----------------------------------------------------
613 * disperr() prints the error message(s) (if any) stored in a disprm struct.
614 * If there are no errors then nothing is printed. It uses wcserr_prt(), q.v.
615 *
616 * Given:
617 * dis const struct disprm*
618 * Distortion function parameters.
619 *
620 * prefix const char *
621 * If non-NULL, each output line will be prefixed with
622 * this string.
623 *
624 * Function return value:
625 * int Status return value:
626 * 0: Success.
627 * 1: Null disprm pointer passed.
628 *
629 *
630 * dishdo() - write FITS headers using TPD
631 * ---------------------------------------
632 * dishdo() sets a flag that tells wcshdo() to write FITS headers in the form
633 * of the TPD translation used internally. Normally SIP and TPV would be
634 * written in their native form if at all possible.
635 *
636 * Given and returned:
637 * dis struct disprm*
638 * Distortion function parameters.
639 *
640 * Function return value:
641 * int Status return value:
642 * 0: Success.
643 * 1: Null disprm pointer passed.
644 * 3: No TPD translation.
645 *
646 *
647 * disset() - Setup routine for the disprm struct
648 * ----------------------------------------------
649 * disset(), sets up the disprm struct according to information supplied within
650 * it - refer to the explanation of disprm::flag.
651 *
652 * Note that this routine need not be called directly; it will be invoked by
653 * disp2x() and disx2p() if the disprm::flag is anything other than a
654 * predefined magic value.
655 *
656 * Given and returned:
657 * dis struct disprm*
658 * Distortion function parameters.
659 *
660 * Function return value:
661 * int Status return value:
662 * 0: Success.
663 * 1: Null disprm pointer passed.
664 * 2: Memory allocation failed.
665 * 3: Invalid parameter.
666 *
667 * For returns > 1, a detailed error message is set in
668 * disprm::err if enabled, see wcserr_enable().
669 *
670 *
671 * disp2x() - Apply distortion function
672 * ------------------------------------
673 * disp2x() applies the distortion functions. By definition, the distortion
674 * is in the pixel-to-world direction.
675 *
676 * Depending on the point in the algorithm chain at which it is invoked,
677 * disp2x() may transform pixel coordinates to corrected pixel coordinates, or
678 * intermediate pixel coordinates to corrected intermediate pixel coordinates,
679 * or image coordinates to corrected image coordinates.
680 *
681 *
682 * Given and returned:
683 * dis struct disprm*
684 * Distortion function parameters.
685 *
686 * Given:
687 * rawcrd const double[naxis]
688 * Array of coordinates.
689 *
690 * Returned:
691 * discrd double[naxis]
692 * Array of coordinates to which the distortion functions
693 * have been applied.
694 *
695 * Function return value:
696 * int Status return value:
697 * 0: Success.
698 * 1: Null disprm pointer passed.
699 * 2: Memory allocation failed.
700 * 3: Invalid parameter.
701 * 4: Distort error.
702 *
703 * For returns > 1, a detailed error message is set in
704 * disprm::err if enabled, see wcserr_enable().
705 *
706 *
707 * disx2p() - Apply de-distortion function
708 * ---------------------------------------
709 * disx2p() applies the inverse of the distortion functions. By definition,
710 * the de-distortion is in the world-to-pixel direction.
711 *
712 * Depending on the point in the algorithm chain at which it is invoked,
713 * disx2p() may transform corrected pixel coordinates to pixel coordinates, or
714 * corrected intermediate pixel coordinates to intermediate pixel coordinates,
715 * or corrected image coordinates to image coordinates.
716 *
717 * disx2p() iteratively solves for the inverse using disp2x(). It assumes
718 * that the distortion is small and the functions are well-behaved, being
719 * continuous and with continuous derivatives. Also that, to first order
720 * in the neighbourhood of the solution, discrd[j] ~= a + b*rawcrd[j], i.e.
721 * independent of rawcrd[i], where i != j. This is effectively equivalent to
722 * assuming that the distortion functions are separable to first order.
723 * Furthermore, a is assumed to be small, and b close to unity.
724 *
725 * If disprm::disx2p() is defined, then disx2p() uses it to provide an initial
726 * estimate for its more precise iterative inversion.
727 *
728 * Given and returned:
729 * dis struct disprm*
730 * Distortion function parameters.
731 *
732 * Given:
733 * discrd const double[naxis]
734 * Array of coordinates.
735 *
736 * Returned:
737 * rawcrd double[naxis]
738 * Array of coordinates to which the inverse distortion
739 * functions have been applied.
740 *
741 * Function return value:
742 * int Status return value:
743 * 0: Success.
744 * 1: Null disprm pointer passed.
745 * 2: Memory allocation failed.
746 * 3: Invalid parameter.
747 * 5: De-distort error.
748 *
749 * For returns > 1, a detailed error message is set in
750 * disprm::err if enabled, see wcserr_enable().
751 *
752 *
753 * diswarp() - Compute measures of distortion
754 * ------------------------------------------
755 * diswarp() computes various measures of the distortion over a specified range
756 * of coordinates.
757 *
758 * For prior distortions, the measures may be interpreted simply as an offset
759 * in pixel coordinates. For sequent distortions, the interpretation depends
760 * on the nature of the linear transformation matrix (PCi_ja or CDi_ja). If
761 * the latter introduces a scaling, then the measures will also be scaled.
762 * Note also that the image domain, which is rectangular in pixel coordinates,
763 * may be rotated, skewed, and/or stretched in intermediate pixel coordinates,
764 * and in general cannot be defined using pixblc[] and pixtrc[].
765 *
766 * PLEASE NOTE: the measures of total distortion may be essentially meaningless
767 * if there are multiple sequent distortions with different scaling.
768 *
769 * See also linwarp().
770 *
771 * Given and returned:
772 * dis struct disprm*
773 * Distortion function parameters.
774 *
775 * Given:
776 * pixblc const double[naxis]
777 * Start of the range of pixel coordinates (for prior
778 * distortions), or intermediate pixel coordinates (for
779 * sequent distortions). May be specified as a NULL
780 * pointer which is interpreted as (1,1,...).
781 *
782 * pixtrc const double[naxis]
783 * End of the range of pixel coordinates (prior) or
784 * intermediate pixel coordinates (sequent).
785 *
786 * pixsamp const double[naxis]
787 * If positive or zero, the increment on the particular
788 * axis, starting at pixblc[]. Zero is interpreted as a
789 * unit increment. pixsamp may also be specified as a
790 * NULL pointer which is interpreted as all zeroes, i.e.
791 * unit increments on all axes.
792 *
793 * If negative, the grid size on the particular axis (the
794 * absolute value being rounded to the nearest integer).
795 * For example, if pixsamp is (-128.0,-128.0,...) then
796 * each axis will be sampled at 128 points between
797 * pixblc[] and pixtrc[] inclusive. Use caution when
798 * using this option on non-square images.
799 *
800 * Returned:
801 * nsamp int* The number of pixel coordinates sampled.
802 *
803 * Can be specified as a NULL pointer if not required.
804 *
805 * maxdis double[naxis]
806 * For each individual distortion function, the
807 * maximum absolute value of the distortion.
808 *
809 * Can be specified as a NULL pointer if not required.
810 *
811 * maxtot double* For the combination of all distortion functions, the
812 * maximum absolute value of the distortion.
813 *
814 * Can be specified as a NULL pointer if not required.
815 *
816 * avgdis double[naxis]
817 * For each individual distortion function, the
818 * mean value of the distortion.
819 *
820 * Can be specified as a NULL pointer if not required.
821 *
822 * avgtot double* For the combination of all distortion functions, the
823 * mean value of the distortion.
824 *
825 * Can be specified as a NULL pointer if not required.
826 *
827 * rmsdis double[naxis]
828 * For each individual distortion function, the
829 * root mean square deviation of the distortion.
830 *
831 * Can be specified as a NULL pointer if not required.
832 *
833 * rmstot double* For the combination of all distortion functions, the
834 * root mean square deviation of the distortion.
835 *
836 * Can be specified as a NULL pointer if not required.
837 *
838 * Function return value:
839 * int Status return value:
840 * 0: Success.
841 * 1: Null disprm pointer passed.
842 * 2: Memory allocation failed.
843 * 3: Invalid parameter.
844 * 4: Distort error.
845 *
846 *
847 * disprm struct - Distortion parameters
848 * -------------------------------------
849 * The disprm struct contains all of the information required to apply a set of
850 * distortion functions. It consists of certain members that must be set by
851 * the user ("given") and others that are set by the WCSLIB routines
852 * ("returned"). While the addresses of the arrays themselves may be set by
853 * disinit() if it (optionally) allocates memory, their contents must be set by
854 * the user.
855 *
856 * int flag
857 * (Given and returned) This flag must be set to zero whenever any of the
858 * following members of the disprm struct are set or modified:
859 *
860 * - disprm::naxis,
861 * - disprm::dtype,
862 * - disprm::ndp,
863 * - disprm::dp.
864 *
865 * This signals the initialization routine, disset(), to recompute the
866 * returned members of the disprm struct. disset() will reset flag to
867 * indicate that this has been done.
868 *
869 * PLEASE NOTE: flag must be set to -1 when disinit() is called for the
870 * first time for a particular disprm struct in order to initialize memory
871 * management. It must ONLY be used on the first initialization otherwise
872 * memory leaks may result.
873 *
874 * int naxis
875 * (Given or returned) Number of pixel and world coordinate elements.
876 *
877 * If disinit() is used to initialize the disprm struct (as would normally
878 * be the case) then it will set naxis from the value passed to it as a
879 * function argument. The user should not subsequently modify it.
880 *
881 * char (*dtype)[72]
882 * (Given) Pointer to the first element of an array of char[72] containing
883 * the name of the distortion function for each axis.
884 *
885 * int ndp
886 * (Given) The number of entries in the disprm::dp[] array.
887 *
888 * int ndpmax
889 * (Given) The length of the disprm::dp[] array.
890 *
891 * ndpmax will be set by disinit() if it allocates memory for disprm::dp[],
892 * otherwise it must be set by the user. See also disndp().
893 *
894 * struct dpkey dp
895 * (Given) Address of the first element of an array of length ndpmax of
896 * dpkey structs.
897 *
898 * As a FITS header parser encounters each DPja or DQia keyword it should
899 * load it into a dpkey struct in the array and increment ndp. However,
900 * note that a single disprm struct must hold only DPja or DQia keyvalues,
901 * not both. disset() interprets them as required by the particular
902 * distortion function.
903 *
904 * double *maxdis
905 * (Given) Pointer to the first element of an array of double specifying
906 * the maximum absolute value of the distortion for each axis computed over
907 * the whole image.
908 *
909 * It is not necessary to reset the disprm struct (via disset()) when
910 * disprm::maxdis is changed.
911 *
912 * double totdis
913 * (Given) The maximum absolute value of the combination of all distortion
914 * functions specified as an offset in pixel coordinates computed over the
915 * whole image.
916 *
917 * It is not necessary to reset the disprm struct (via disset()) when
918 * disprm::totdis is changed.
919 *
920 * int *docorr
921 * (Returned) Pointer to the first element of an array of int containing
922 * flags that indicate the mode of correction for each axis.
923 *
924 * If docorr is zero, the distortion function returns the corrected
925 * coordinates directly. Any other value indicates that the distortion
926 * function computes a correction to be added to pixel coordinates (prior
927 * distortion) or intermediate pixel coordinates (sequent distortion).
928 *
929 * int *Nhat
930 * (Returned) Pointer to the first element of an array of int containing
931 * the number of coordinate axes that form the independent variables of the
932 * distortion function for each axis.
933 *
934 * int **axmap
935 * (Returned) Pointer to the first element of an array of int* containing
936 * pointers to the first elements of the axis mapping arrays for each axis.
937 *
938 * An axis mapping associates the independent variables of a distortion
939 * function with the 0-relative image axis number. For example, consider
940 * an image with a spectrum on the first axis (axis 0), followed by RA
941 * (axis 1), Dec (axis2), and time (axis 3) axes. For a distortion in
942 * (RA,Dec) and no distortion on the spectral or time axes, the axis
943 * mapping arrays, axmap[j][], would be
944 *
945 = j=0: [-1, -1, -1, -1] ...no distortion on spectral axis,
946 = 1: [ 1, 2, -1, -1] ...RA distortion depends on RA and Dec,
947 = 2: [ 2, 1, -1, -1] ...Dec distortion depends on Dec and RA,
948 = 3: [-1, -1, -1, -1] ...no distortion on time axis,
949 *
950 * where -1 indicates that there is no corresponding independent
951 * variable.
952 *
953 * double **offset
954 * (Returned) Pointer to the first element of an array of double*
955 * containing pointers to the first elements of arrays of offsets used to
956 * renormalize the independent variables of the distortion function for
957 * each axis.
958 *
959 * The offsets are subtracted from the independent variables before
960 * scaling.
961 *
962 * double **scale
963 * (Returned) Pointer to the first element of an array of double*
964 * containing pointers to the first elements of arrays of scales used to
965 * renormalize the independent variables of the distortion function for
966 * each axis.
967 *
968 * The scale is applied to the independent variables after the offsets are
969 * subtracted.
970 *
971 * int **iparm
972 * (Returned) Pointer to the first element of an array of int*
973 * containing pointers to the first elements of the arrays of integer
974 * distortion parameters for each axis.
975 *
976 * double **dparm
977 * (Returned) Pointer to the first element of an array of double*
978 * containing pointers to the first elements of the arrays of floating
979 * point distortion parameters for each axis.
980 *
981 * int i_naxis
982 * (Returned) Dimension of the internal arrays (normally equal to naxis).
983 *
984 * int ndis
985 * (Returned) The number of distortion functions.
986 *
987 * struct wcserr *err
988 * (Returned) If enabled, when an error status is returned, this struct
989 * contains detailed information about the error, see wcserr_enable().
990 *
991 * int (**disp2x)(DISP2X_ARGS)
992 * (For internal use only.)
993 * int (**disx2p)(DISX2P_ARGS)
994 * (For internal use only.)
995 * double *dummy
996 * (For internal use only.)
997 * int m_flag
998 * (For internal use only.)
999 * int m_naxis
1000 * (For internal use only.)
1001 * char (*m_dtype)[72]
1002 * (For internal use only.)
1003 * double **m_dp
1004 * (For internal use only.)
1005 * double *m_maxdis
1006 * (For internal use only.)
1007 *
1008 *
1009 * dpkey struct - Store for DPja and DQia keyvalues
1010 * ------------------------------------------------
1011 * The dpkey struct is used to pass the parsed contents of DPja or DQia
1012 * keyrecords to disset() via the disprm struct. A disprm struct must hold
1013 * only DPja or DQia keyvalues, not both.
1014 *
1015 * All members of this struct are to be set by the user.
1016 *
1017 * char field[72]
1018 * (Given) The full field name of the record, including the keyword name.
1019 * Note that the colon delimiter separating the field name and the value in
1020 * record-valued keyvalues is not part of the field name. For example, in
1021 * the following:
1022 *
1023 = DP3A = 'AXIS.1: 2'
1024 *
1025 * the full record field name is "DP3A.AXIS.1", and the record's value
1026 * is 2.
1027 *
1028 * int j
1029 * (Given) Axis number (1-relative), i.e. the j in DPja or i in DQia.
1030 *
1031 * int type
1032 * (Given) The data type of the record's value
1033 * - 0: Integer (stored as an int),
1034 * - 1: Floating point (stored as a double).
1035 *
1036 * union value
1037 * (Given) A union comprised of
1038 * - dpkey::i,
1039 * - dpkey::f,
1040 *
1041 * the record's value.
1042 *
1043 *
1044 * Global variable: const char *dis_errmsg[] - Status return messages
1045 * ------------------------------------------------------------------
1046 * Error messages to match the status value returned from each function.
1047 *
1048 *===========================================================================*/
1049
1050 #ifndef WCSLIB_DIS
1051 #define WCSLIB_DIS
1052
1053 #ifdef __cplusplus
1054 extern "C" {
1055 #endif
1056
1057
1058 extern const char *dis_errmsg[];
1059
1060 enum dis_errmsg_enum {
1061 DISERR_SUCCESS = 0, // Success.
1062 DISERR_NULL_POINTER = 1, // Null disprm pointer passed.
1063 DISERR_MEMORY = 2, // Memory allocation failed.
1064 DISERR_BAD_PARAM = 3, // Invalid parameter value.
1065 DISERR_DISTORT = 4, // Distortion error.
1066 DISERR_DEDISTORT = 5 // De-distortion error.
1067 };
1068
1069 // For use in declaring distortion function prototypes (= DISX2P_ARGS).
1070 #define DISP2X_ARGS int inverse, const int iparm[], const double dparm[], \
1071 int ncrd, const double rawcrd[], double *discrd
1072
1073 // For use in declaring de-distortion function prototypes (= DISP2X_ARGS).
1074 #define DISX2P_ARGS int inverse, const int iparm[], const double dparm[], \
1075 int ncrd, const double discrd[], double *rawcrd
1076
1077
1078 // Struct used for storing DPja and DQia keyvalues.
1079 struct dpkey {
1080 char field[72]; // Full record field name (no colon).
1081 int j; // Axis number, as in DPja (1-relative).
1082 int type; // Data type of value.
1083 union {
1084 int i; // Integer record value.
1085 double f; // Floating point record value.
1086 } value; // Record value.
1087 };
1088
1089 // Size of the dpkey struct in int units, used by the Fortran wrappers.
1090 #define DPLEN (sizeof(struct dpkey)/sizeof(int))
1091
1092
1093 struct disprm {
1094 // Initialization flag (see the prologue above).
1095 //--------------------------------------------------------------------------
1096 int flag; // Set to zero to force initialization.
1097
1098 // Parameters to be provided (see the prologue above).
1099 //--------------------------------------------------------------------------
1100 int naxis; // The number of pixel coordinate elements,
1101 // given by NAXIS.
1102 char (*dtype)[72]; // For each axis, the distortion type.
1103 int ndp; // Number of DPja or DQia keywords, and the
1104 int ndpmax; // number for which space was allocated.
1105 struct dpkey *dp; // DPja or DQia keyvalues (not both).
1106 double totdis; // The maximum combined distortion.
1107 double *maxdis; // For each axis, the maximum distortion.
1108
1109 // Information derived from the parameters supplied.
1110 //--------------------------------------------------------------------------
1111 int *docorr; // For each axis, the mode of correction.
1112 int *Nhat; // For each axis, the number of coordinate
1113 // axes that form the independent variables
1114 // of the distortion function.
1115 int **axmap; // For each axis, the axis mapping array.
1116 double **offset; // For each axis, renormalization offsets.
1117 double **scale; // For each axis, renormalization scales.
1118 int **iparm; // For each axis, the array of integer
1119 // distortion parameters.
1120 double **dparm; // For each axis, the array of floating
1121 // point distortion parameters.
1122 int i_naxis; // Dimension of the internal arrays.
1123 int ndis; // The number of distortion functions.
1124
1125 // Error handling, if enabled.
1126 //--------------------------------------------------------------------------
1127 struct wcserr *err;
1128
1129 // Private - the remainder are for internal use.
1130 //--------------------------------------------------------------------------
1131 int (**disp2x)(DISP2X_ARGS); // For each axis, pointers to the
1132 int (**disx2p)(DISX2P_ARGS); // distortion function and its inverse.
1133
1134 int m_flag, m_naxis; // The remainder are for memory management.
1135 char (*m_dtype)[72];
1136 struct dpkey *m_dp;
1137 double *m_maxdis;
1138 };
1139
1140 // Size of the disprm struct in int units, used by the Fortran wrappers.
1141 #define DISLEN (sizeof(struct disprm)/sizeof(int))
1142
1143
1144 int disndp(int n);
1145
1146 int dpfill(struct dpkey *dp, const char *keyword, const char *field, int j,
1147 int type, int i, double f);
1148
1149 int dpkeyi(const struct dpkey *dp);
1150
1151 double dpkeyd(const struct dpkey *dp);
1152
1153 int disini(int alloc, int naxis, struct disprm *dis);
1154
1155 int disinit(int alloc, int naxis, struct disprm *dis, int ndpmax);
1156
1157 int discpy(int alloc, const struct disprm *dissrc, struct disprm *disdst);
1158
1159 int disfree(struct disprm *dis);
1160
1161 int dissize(const struct disprm *dis, int sizes[2]);
1162
1163 int disprt(const struct disprm *dis);
1164
1165 int disperr(const struct disprm *dis, const char *prefix);
1166
1167 int dishdo(struct disprm *dis);
1168
1169 int disset(struct disprm *dis);
1170
1171 int disp2x(struct disprm *dis, const double rawcrd[], double discrd[]);
1172
1173 int disx2p(struct disprm *dis, const double discrd[], double rawcrd[]);
1174
1175 int diswarp(struct disprm *dis, const double pixblc[], const double pixtrc[],
1176 const double pixsamp[], int *nsamp,
1177 double maxdis[], double *maxtot,
1178 double avgdis[], double *avgtot,
1179 double rmsdis[], double *rmstot);
1180
1181 #ifdef __cplusplus
1182 }
1183 #endif
1184
1185 #endif // WCSLIB_DIS
Browse
Usage
Versions