WCSLIB 8.2.2
Loading...
Searching...
No Matches
dis.h
Go to the documentation of this file.
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
1054extern "C" {
1055#endif
1056
1057
1058extern const char *dis_errmsg[];
1059
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.
1068
1069// For use in declaring distortion function prototypes (= DISX2P_ARGS).
1070#define DISP2X_ARGS int inverse, const int iparm[], const double dparm[], \
1071int 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[], \
1075int ncrd, const double discrd[], double *rawcrd
1076
1077
1078// Struct used for storing DPja and DQia keyvalues.
1079struct 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
1093struct 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
1144int disndp(int n);
1145
1146int dpfill(struct dpkey *dp, const char *keyword, const char *field, int j,
1147 int type, int i, double f);
1148
1149int dpkeyi(const struct dpkey *dp);
1150
1151double dpkeyd(const struct dpkey *dp);
1152
1153int disini(int alloc, int naxis, struct disprm *dis);
1154
1155int disinit(int alloc, int naxis, struct disprm *dis, int ndpmax);
1156
1157int discpy(int alloc, const struct disprm *dissrc, struct disprm *disdst);
1158
1159int disfree(struct disprm *dis);
1160
1161int dissize(const struct disprm *dis, int sizes[2]);
1162
1163int disprt(const struct disprm *dis);
1164
1165int disperr(const struct disprm *dis, const char *prefix);
1166
1167int dishdo(struct disprm *dis);
1168
1169int disset(struct disprm *dis);
1170
1171int disp2x(struct disprm *dis, const double rawcrd[], double discrd[]);
1172
1173int disx2p(struct disprm *dis, const double discrd[], double rawcrd[]);
1174
1175int 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
int disp2x(struct disprm *dis, const double rawcrd[], double discrd[])
Apply distortion function.
double dpkeyd(const struct dpkey *dp)
Get the data value in a dpkey struct as double.
#define DISP2X_ARGS
Definition dis.h:1070
dis_errmsg_enum
Definition dis.h:1060
@ DISERR_MEMORY
Definition dis.h:1063
@ DISERR_NULL_POINTER
Definition dis.h:1062
@ DISERR_DEDISTORT
Definition dis.h:1066
@ DISERR_BAD_PARAM
Definition dis.h:1064
@ DISERR_DISTORT
Definition dis.h:1065
@ DISERR_SUCCESS
Definition dis.h:1061
int disinit(int alloc, int naxis, struct disprm *dis, int ndpmax)
Default constructor for the disprm struct.
int discpy(int alloc, const struct disprm *dissrc, struct disprm *disdst)
Copy routine for the disprm struct.
int dpfill(struct dpkey *dp, const char *keyword, const char *field, int j, int type, int i, double f)
Fill the contents of a dpkey struct.
int dpkeyi(const struct dpkey *dp)
Get the data value in a dpkey struct as int.
int dissize(const struct disprm *dis, int sizes[2])
Compute the size of a disprm struct.
const char * dis_errmsg[]
Status return messages.
int disfree(struct disprm *dis)
Destructor for the disprm struct.
int disini(int alloc, int naxis, struct disprm *dis)
Default constructor for the disprm struct.
int diswarp(struct disprm *dis, const double pixblc[], const double pixtrc[], const double pixsamp[], int *nsamp, double maxdis[], double *maxtot, double avgdis[], double *avgtot, double rmsdis[], double *rmstot)
Compute measures of distortion.
int disndp(int n)
Memory allocation for DPja and DQia.
int dishdo(struct disprm *dis)
write FITS headers using TPD.
int disx2p(struct disprm *dis, const double discrd[], double rawcrd[])
Apply de-distortion function.
#define DISX2P_ARGS
Definition dis.h:1074
int disperr(const struct disprm *dis, const char *prefix)
Print error messages from a disprm struct.
int disset(struct disprm *dis)
Setup routine for the disprm struct.
int disprt(const struct disprm *dis)
Print routine for the disprm struct.
Distortion parameters.
Definition dis.h:1093
double totdis
Definition dis.h:1106
int ** iparm
Definition dis.h:1118
double ** scale
Definition dis.h:1117
int naxis
Definition dis.h:1100
struct wcserr * err
Definition dis.h:1127
char(* dtype)[72]
Definition dis.h:1102
int flag
Definition dis.h:1096
double * m_maxdis
Definition dis.h:1137
int i_naxis
Definition dis.h:1122
int(** disp2x)(DISP2X_ARGS)
Definition dis.h:1131
int ndpmax
Definition dis.h:1104
int m_naxis
Definition dis.h:1134
struct dpkey * dp
Definition dis.h:1105
int ndis
Definition dis.h:1123
int * docorr
Definition dis.h:1111
double ** offset
Definition dis.h:1116
char(* m_dtype)[72]
Definition dis.h:1135
double ** dparm
Definition dis.h:1120
int ** axmap
Definition dis.h:1115
struct dpkey * m_dp
Definition dis.h:1136
double * maxdis
Definition dis.h:1107
int ndp
Definition dis.h:1103
int(** disx2p)(DISX2P_ARGS)
Definition dis.h:1132
int * Nhat
Definition dis.h:1112
int m_flag
Definition dis.h:1134
Store for DPja and DQia keyvalues.
Definition dis.h:1079
union dpkey::@0 value
int i
Definition dis.h:1084
char field[72]
Definition dis.h:1080
int type
Definition dis.h:1082
int j
Definition dis.h:1081
double f
Definition dis.h:1085
Error message handling.
Definition wcserr.h:243