From 183f7b446e311ae2247d2d73720539eede9ba37d Mon Sep 17 00:00:00 2001 From: cgarling Date: Mon, 7 Sep 2026 07:50:50 -0400 Subject: [PATCH] Fix TPV/TPD implementation --- ext/FITSWCSFITSFilesExt.jl | 3 +- ext/FITSWCSFITSIOExt.jl | 3 +- src/api.jl | 4 +- src/auxiliary_data.jl | 26 +++- src/distortion.jl | 276 ++++++++++++++++++++++++++++++++++--- src/linear.jl | 56 ++++++-- src/parsing.jl | 81 ++++++++--- src/projections.jl | 28 ++-- src/transforms.jl | 4 +- test/runtests.jl | 37 ++--- test/slicing_tests.jl | 2 +- test/tpv_tpd_tests.jl | 191 +++++++++++++++++++++++++ 12 files changed, 621 insertions(+), 90 deletions(-) create mode 100644 test/tpv_tpd_tests.jl diff --git a/ext/FITSWCSFITSFilesExt.jl b/ext/FITSWCSFITSFilesExt.jl index dfd21fa..900c0d1 100644 --- a/ext/FITSWCSFITSFilesExt.jl +++ b/ext/FITSWCSFITSFilesExt.jl @@ -35,7 +35,8 @@ function _set_fitsfiles_header_card!(dict::Dict{String, Any}, key::String, value isnothing(parsed) && return dict name, parsed_value = parsed - if name == "EXTVER" || name == "NAXES" || startswith(name, "AXIS.") + if name in ("EXTVER", "NAXES", "DOCORR") || + any(prefix -> startswith(name, prefix), ("AXIS.", "OFFSET.", "SCALE.", "TPD.", "AUX.")) dict["$(key).$(name)"] = parsed_value end return dict diff --git a/ext/FITSWCSFITSIOExt.jl b/ext/FITSWCSFITSIOExt.jl index 0c62a69..1d51f2e 100644 --- a/ext/FITSWCSFITSIOExt.jl +++ b/ext/FITSWCSFITSIOExt.jl @@ -35,7 +35,8 @@ function _set_fitsio_header_card!(dict::Dict{String, Any}, key::String, value) isnothing(parsed) && return dict name, parsed_value = parsed - if name == "EXTVER" || name == "NAXES" || startswith(name, "AXIS.") + if name in ("EXTVER", "NAXES", "DOCORR") || + any(prefix -> startswith(name, prefix), ("AXIS.", "OFFSET.", "SCALE.", "TPD.", "AUX.")) dict["$(key).$(name)"] = parsed_value end return dict diff --git a/src/api.jl b/src/api.jl index 3d3f172..c3c5a16 100644 --- a/src/api.jl +++ b/src/api.jl @@ -33,9 +33,7 @@ matching the APE-14 convention. For `SlicedWCSTransform`: the parent's matrix indexed by the kept axes. """ function axis_correlation_matrix(wcs::WCSTransform{N}) where {N} - # If any pre-linear distortion is present (SIP, D2IM, CPDIS), we must be - # conservative: distortions can introduce arbitrary cross-axis coupling - # that isn't captured by the CD matrix alone. + # Any distortion stage may introduce cross-axis coupling not captured by CD. if has_distortion(wcs.pipeline) return fill(true, SMatrix{N, N, Bool}) end diff --git a/src/auxiliary_data.jl b/src/auxiliary_data.jl index 862325e..1ace3eb 100644 --- a/src/auxiliary_data.jl +++ b/src/auxiliary_data.jl @@ -67,7 +67,9 @@ function _paper_iv_lookup_spec( ) dist_key = "$(dist_prefix)$(axis)$(alt_str)" haskey(header, dist_key) || return nothing - uppercase(String(header[dist_key])) == "LOOKUP" || + dtype = uppercase(strip(String(header[dist_key]))) + dtype == "TPD" && return nothing + dtype == "LOOKUP" || throw(ArgumentError("unsupported Paper IV distortion type $(header[dist_key]) in $dist_key")) # Skip tables whose declared error is below the requested threshold. @@ -154,10 +156,26 @@ function _header_references_tabular_axis(header::AbstractDict, alt_str::Abstract end function _header_references_external_wcs_data(header::AbstractDict, alt_str::AbstractString) - # Paper IV distortion keywords and TAB axes are the current external-data signals. - for key in keys(header) + suffix = uppercase(String(alt_str)) + + # Only lookup distortions require external arrays; TPD is header-resident. + for (key, value) in header key isa AbstractString || continue - is_lookup_distortion_keyword(key, alt_str) && return true + ukey = uppercase(String(key)) + occursin(Regex("^D2IMDIS[1-9][0-9]*$(suffix)\$"), ukey) && return true + occursin(Regex("^D2IMERR[1-9][0-9]*$(suffix)\$"), ukey) && return true + ukey == "AXISCORR$(suffix)" && return true + if occursin(Regex("^CPDIS[1-9][0-9]*$(suffix)\$"), ukey) + uppercase(strip(String(value))) == "LOOKUP" && return true + end + + # Undeclared DP/DQ records retain the legacy external-data guard, while + # records belonging to an explicit TPD function remain header-resident. + parameter = match(Regex("^D([PQ])([1-9][0-9]*)$(suffix)\\."), ukey) + isnothing(parameter) && continue + dist_prefix = parameter.captures[1] == "P" ? "CPDIS" : "CQDIS" + dist_key = "$(dist_prefix)$(parameter.captures[2])$(suffix)" + haskey(header, dist_key) && uppercase(strip(String(header[dist_key]))) == "TPD" || return true end return _header_references_tabular_axis(header, alt_str) diff --git a/src/distortion.jl b/src/distortion.jl index 519da04..a616f91 100644 --- a/src/distortion.jl +++ b/src/distortion.jl @@ -23,46 +23,66 @@ struct SIPDistortion bp::Union{Nothing, Matrix{Float64}} end +""" + TPDDistortion + +One axis of a Paper IV Template Polynomial Distortion. `axmap`, `offset`, and +`scale` define the independent variables. Forward and optional reverse +coefficients use the 60-term TPD basis. By default the polynomial is an +additive correction; `docorr=false` makes it return the corrected coordinate +directly. +""" +struct TPDDistortion{M, A <: Union{Nothing, SMatrix{2, 3, Float64, 6}}} + axmap::SVector{M, Int} + offset::SVector{M, Float64} + scale::SVector{M, Float64} + aux::A + forward::Vector{Float64} + reverse::Union{Nothing, Vector{Float64}} + docorr::Bool +end + # ────────────────────────────────────────────────────────────────────────────── -"""Abstract supertype for pre-linear pixel/focal-plane distortion pipelines.""" +"""Abstract supertype for distortion pipelines surrounding the linear transform.""" abstract type AbstractDistortionPipeline end -"""Identity distortion pipeline for WCS transforms with no pre-linear distortion.""" +"""Identity distortion pipeline for WCS transforms with no distortion stages.""" struct NoDistortionPipeline <: AbstractDistortionPipeline end """ DistortionPipeline -Pre-linear distortion pipeline. Detector-to-image lookup tables are applied -first; SIP and CPDIS offsets are then evaluated at the detector-corrected -coordinate before the linear WCS matrix. +Distortion pipeline surrounding the linear WCS matrix. Detector-to-image +lookup tables, SIP, CPDIS lookup offsets, and prior TPD are applied before PC. +Sequent TPD is applied after PC and before CDELT. """ -struct DistortionPipeline{S <: Union{Nothing, SIPDistortion}, D <: Tuple, C <: Tuple} <: AbstractDistortionPipeline +struct DistortionPipeline{S <: Union{Nothing, SIPDistortion}, D <: Tuple, C <: Tuple, P <: Tuple, Q <: Tuple} <: AbstractDistortionPipeline det2im::D sip::S cpdis::C + tpd_pre::P + tpd_seq::Q end DistortionPipeline(sip::SIPDistortion) = - DistortionPipeline{typeof(sip), Tuple{Nothing, Nothing}, Tuple{Nothing, Nothing}}((nothing, nothing), sip, (nothing, nothing)) - -function distortion_pipeline(sip::Union{Nothing, SIPDistortion}, ::NoAuxiliaryWCSData) - # Header-only WCS uses the existing no-op or SIP-only pipeline. - return distortion_pipeline(sip) -end + DistortionPipeline((nothing, nothing), sip, (nothing, nothing), (), ()) -function distortion_pipeline(sip::Union{Nothing, SIPDistortion}, aux::AuxiliaryWCSData) - det2im = aux.det2im isa Tuple && length(aux.det2im) == 2 ? aux.det2im : (nothing, nothing) - cpdis = aux.cpdis isa Tuple && length(aux.cpdis) == 2 ? aux.cpdis : (nothing, nothing) +function distortion_pipeline(sip::Union{Nothing, SIPDistortion}, aux::AbstractAuxiliaryWCSData, tpd_pre::Tuple, tpd_seq::Tuple) + # External lookup tables are absent from header-only auxiliary payloads. + det2im = aux isa AuxiliaryWCSData && aux.det2im isa Tuple && length(aux.det2im) == 2 ? aux.det2im : (nothing, nothing) + cpdis = aux isa AuxiliaryWCSData && aux.cpdis isa Tuple && length(aux.cpdis) == 2 ? aux.cpdis : (nothing, nothing) - # Avoid allocating a pipeline for auxiliary payloads that contain only TAB data. - if isnothing(sip) && all(isnothing, det2im) && all(isnothing, cpdis) + # Avoid allocating a pipeline when every distortion stage is empty. + if isnothing(sip) && all(isnothing, det2im) && all(isnothing, cpdis) && + all(isnothing, tpd_pre) && all(isnothing, tpd_seq) return NoDistortionPipeline() end - # Preserve exact tuple types so no-lookup and partial-lookup stages stay concrete. - return DistortionPipeline{typeof(sip), typeof(det2im), typeof(cpdis)}(det2im, sip, cpdis) + # Preserve empty TPD stages explicitly for cheap common-case dispatch. + prior = all(isnothing, tpd_pre) ? () : tpd_pre + sequent = all(isnothing, tpd_seq) ? () : tpd_seq + return DistortionPipeline(det2im, sip, cpdis, prior, sequent) end function has_sip_keywords(header::AbstractDict, alt_str::AbstractString) @@ -128,6 +148,107 @@ function parse_sip_distortion(header::AbstractDict, crpix::Vector{Float64}, naxi return SIPDistortion(SVector{2, Float64}(crpix[1:2]), a, b, ap, bp) end +function _collect_tpd_coefficients(header::AbstractDict, prefix::AbstractString, direction::AbstractString) + coeff = Float64[] + + # Preserve coefficient numbering while zero-filling omitted terms. + for m in 0:59 + key = "$(prefix).TPD.$(direction).$(m)" + haskey(header, key) || continue + while length(coeff) < m + push!(coeff, 0.0) + end + push!(coeff, Float64(header[key])) + end + return coeff +end + +function _tpd_auxiliary_matrix(header::AbstractDict, prefix::AbstractString) + has_aux = any(haskey(header, "$(prefix).AUX.$(k).COEFF.$(m)") for k in 1:2, m in 0:2) + has_aux || return nothing + + # WCSLIB defaults the auxiliary transform to the two-dimensional identity. + aux = zeros(Float64, 2, 3) + aux[1, 2] = 1.0 + aux[2, 3] = 1.0 + for k in 1:2, m in 0:2 + key = "$(prefix).AUX.$(k).COEFF.$(m)" + haskey(header, key) && (aux[k, m + 1] = Float64(header[key])) + end + return SMatrix{2, 3, Float64, 6}(aux) +end + +function _valid_tpd_parameter_field(field::AbstractString) + field in ("DOCORR", "NAXES") && return true + occursin(r"^(AXIS|OFFSET|SCALE)\.[1-9][0-9]*$", field) && return true + occursin(r"^TPD\.(FWD|REV)\.([0-9]|[1-5][0-9])$", field) && return true + return occursin(r"^AUX\.[12]\.COEFF\.[0-2]$", field) +end + +function _build_tpd_distortion(header::AbstractDict, prefix::AbstractString, naxis::Int, ::Val{M}) where {M} + # Resolve the independent-axis map and its normalization parameters. + axmap = SVector{M, Int}(ntuple(k -> Int(get(header, "$(prefix).AXIS.$(k)", k)), M)) + all(i -> 1 <= i <= naxis, axmap) || throw(ArgumentError("$prefix axis map contains an axis outside 1:$naxis")) + length(unique(axmap)) == M || throw(ArgumentError("$prefix axis map contains duplicate axes")) + offset = SVector{M, Float64}(ntuple(k -> Float64(get(header, "$(prefix).OFFSET.$(k)", 0.0)), M)) + scale = SVector{M, Float64}(ntuple(k -> Float64(get(header, "$(prefix).SCALE.$(k)", 1.0)), M)) + + # TPD forward coefficients are required structurally but may all be zero. + forward = _collect_tpd_coefficients(header, prefix, "FWD") + reverse_raw = _collect_tpd_coefficients(header, prefix, "REV") + reverse = isempty(reverse_raw) ? nothing : reverse_raw + aux = _tpd_auxiliary_matrix(header, prefix) + docorr = Int(get(header, "$(prefix).DOCORR", 1)) != 0 + return TPDDistortion(axmap, offset, scale, aux, forward, reverse, docorr) +end + +function _parse_tpd_axis(header::AbstractDict, axis::Int, naxis::Int, alt_str::AbstractString, + dist_prefix::AbstractString, param_prefix::AbstractString, + error_prefix::AbstractString, minerr::Real) + dist_key = "$(dist_prefix)$(axis)$(alt_str)" + haskey(header, dist_key) || return nothing + dtype = uppercase(strip(String(header[dist_key]))) + if dtype == "LOOKUP" + dist_prefix == "CPDIS" && return nothing + throw(ArgumentError("sequent LOOKUP distortion $dist_key is not supported")) + end + dtype == "TPD" || throw(ArgumentError("unsupported Paper IV distortion type $(header[dist_key]) in $dist_key")) + + # Error thresholds suppress the complete distortion on this axis. + error_key = "$(error_prefix)$(axis)$(alt_str)" + Float64(get(header, error_key, 0.0)) < Float64(minerr) && return nothing + + prefix = "$(param_prefix)$(axis)$(alt_str)" + naxes_key = "$(prefix).NAXES" + haskey(header, naxes_key) || throw(ArgumentError("TPD distortion $dist_key requires $naxes_key")) + nhat = Int(header[naxes_key]) + nhat in (1, 2) || throw(ArgumentError("$naxes_key must be 1 or 2, got $nhat")) + + # Reject misspelled or unsupported TPD parameter records. + dotted_prefix = "$(prefix)." + for key in keys(header) + key isa AbstractString || continue + startswith(key, dotted_prefix) || continue + field = key[length(dotted_prefix) + 1:end] + _valid_tpd_parameter_field(field) || throw(ArgumentError("unrecognized TPD parameter $key")) + mapped = match(r"^(AXIS|OFFSET|SCALE)\.([1-9][0-9]*)$", field) + mapped !== nothing && parse(Int, mapped.captures[2]) > nhat && + throw(ArgumentError("$key exceeds the $nhat independent axes declared by $naxes_key")) + end + + return nhat == 1 ? _build_tpd_distortion(header, prefix, naxis, Val(1)) : + _build_tpd_distortion(header, prefix, naxis, Val(2)) +end + +function parse_tpd_distortions(header::AbstractDict, naxis::Int, alt::Char, minerr::Real) + alt_str = alt == ' ' ? "" : string(alt) + + # CPDIS/DP is prior to PC; CQDIS/DQ is sequent to PC and prior to CDELT. + prior = ntuple(axis -> _parse_tpd_axis(header, axis, naxis, alt_str, "CPDIS", "DP", "CPERR", minerr), naxis) + sequent = ntuple(axis -> _parse_tpd_axis(header, axis, naxis, alt_str, "CQDIS", "DQ", "CQERR", minerr), naxis) + return prior, sequent +end + distortion_pipeline(::Nothing) = NoDistortionPipeline() distortion_pipeline(sip::SIPDistortion) = DistortionPipeline(sip) @@ -164,6 +285,107 @@ function evaluate_sip_polynomial(coeff::AbstractMatrix, u::Real, v::Real) return value end +function _evaluate_tpd_coefficients(coeff::AbstractVector, u::T, v::T, ::Val{M}) where {T, M} + result = zero(T) + + # One-dimensional TPD uses only pure powers of its attached coordinate. + @inbounds for i in eachindex(coeff) + c = coeff[i] + iszero(c) && continue + term = _TPD_TERMS[i] + if M == 1 + term[1] === :mono && term[3] == 0 || continue + end + result += T(c) * _tpd_term_value(i - 1, u, v) + end + return result +end + +function _evaluate_tpd(model::TPDDistortion{M}, raw::StaticVector, coeff::AbstractVector) where {M} + T = _coordinate_float_type(raw) + + # Normalize and reorder the independent coordinates through the axis map. + vars = SVector{M, T}(ntuple(k -> (T(raw[model.axmap[k]]) - T(model.offset[k])) * T(model.scale[k]), M)) + u = vars[1] + v = M == 1 ? zero(T) : vars[2] + + # Optional auxiliary variables are a two-dimensional affine transform. + if !isnothing(model.aux) + aux = SMatrix{2, 3, T, 6}(model.aux) + transformed = aux * SVector{3, T}(one(T), u, v) + u, v = transformed + end + return _evaluate_tpd_coefficients(coeff, u, v, Val(M)) +end + +apply_tpd_stage(::Tuple{}, raw::StaticVector) = raw + +function apply_tpd_stage(models::Tuple, raw::StaticVector{N, T}) where {N, T} + length(models) == N || throw(DimensionMismatch("TPD stage has $(length(models)) axes, expected $N")) + + # Every output is evaluated from the same undistorted input coordinate. + return SVector{N, T}(ntuple(j -> begin + model = models[j] + if isnothing(model) + raw[j] + else + value = _evaluate_tpd(model, raw, model.forward) + model.docorr ? raw[j] + value : value + end + end, N)) +end + +function _tpd_reverse_guess(models::Tuple, target::StaticVector{N, T}) where {N, T} + # Reverse polynomials provide only the starting point for forward iteration. + return SVector{N, T}(ntuple(j -> begin + model = models[j] + if isnothing(model) || isnothing(model.reverse) + target[j] + else + value = _evaluate_tpd(model, target, model.reverse) + model.docorr ? target[j] + value : value + end + end, N)) +end + +invert_tpd_stage(::Tuple{}, target::StaticVector) = target + +function invert_tpd_stage(models::Tuple, target::StaticVector{N, T}) where {N, T} + raw = _tpd_reverse_guess(models, target) + tol = T(_convergence_tol(T)) + + # Refine the reverse-polynomial estimate using the forward model and a numerical Jacobian. + for _ in 1:30 + distorted = apply_tpd_stage(models, raw) + residual = distorted - target + all(j -> abs(residual[j]) <= tol * max(one(T), abs(target[j])), 1:N) && return raw + + jacobian = MMatrix{N, N, T}(undef) + for column in 1:N + step = clamp(abs(residual[column]) / 2, T(1e-6), one(T)) + trial = Base.setindex(raw, raw[column] + step, column) + shifted = apply_tpd_stage(models, trial) + for row in 1:N + jacobian[row, column] = (shifted[row] - distorted[row]) / step + end + end + raw -= SMatrix{N, N, T}(jacobian) \ residual + end + + residual = apply_tpd_stage(models, raw) - target + @warn "TPD inverse failed to converge after 30 iterations (residual $(sqrt(sum(abs2, residual))) > tolerance $tol); returning best estimate" + return raw +end + +has_sequent_distortion(::NoDistortionPipeline) = false +has_sequent_distortion(pipeline::DistortionPipeline) = any(!isnothing, pipeline.tpd_seq) + +apply_sequent_distortion(::NoDistortionPipeline, coord::StaticVector) = coord +apply_sequent_distortion(pipeline::DistortionPipeline, coord::StaticVector) = apply_tpd_stage(pipeline.tpd_seq, coord) + +invert_sequent_distortion(::NoDistortionPipeline, coord::StaticVector) = coord +invert_sequent_distortion(pipeline::DistortionPipeline, coord::StaticVector) = invert_tpd_stage(pipeline.tpd_seq, coord) + sip_pixel_to_focal(::Nothing, pixel::AbstractVector) = SVector{2, _coordinate_float_type(pixel)}(pixel[1], pixel[2]) sip_pixel_to_focal(::Nothing, pixel::StaticVector) = pixel function sip_pixel_to_focal(sip::SIPDistortion, pixel::AbstractVector) @@ -266,7 +488,12 @@ function pixel_to_focal(pipeline::DistortionPipeline, pixel::StaticVector{N}, :: i == 2 ? T(fy) : detector[i], N)) end - return coord + _lookup_stage_offset(pipeline.cpdis, detector) + lookup_coord = coord + _lookup_stage_offset(pipeline.cpdis, detector) + isempty(pipeline.tpd_pre) && return lookup_coord + + # CPDIS functions on different output axes are evaluated from the same input. + tpd_coord = apply_tpd_stage(pipeline.tpd_pre, coord) + return SVector{N, T}(ntuple(i -> isnothing(pipeline.tpd_pre[i]) ? lookup_coord[i] : tpd_coord[i], N)) end function focal_to_pixel(::NoDistortionPipeline, focal::AbstractVector, ::Val{N}) where {N} @@ -294,7 +521,7 @@ end # No Paper IV lookup stage is present, so we invert only SIP stage. # This dispatch is needed to avoid allocations in the common case of SIP-only distortion. -function focal_to_pixel(pipeline::DistortionPipeline{S, Tuple{Nothing, Nothing}, Tuple{Nothing, Nothing}}, focal::StaticVector{N}, ::Val{N}) where {S, N} +function focal_to_pixel(pipeline::DistortionPipeline{S, Tuple{Nothing, Nothing}, Tuple{Nothing, Nothing}, Tuple{}, Q}, focal::StaticVector{N}, ::Val{N}) where {S, Q, N} T = _coordinate_float_type(focal) # Preserve identity behavior for SIP-free pipeline variants. @@ -310,6 +537,13 @@ function focal_to_pixel(pipeline::DistortionPipeline{S, Tuple{Nothing, Nothing}, T(focal[i]), N)) end +# A TPD-only prior stage can use its reverse polynomial and Newton refinement directly. +function focal_to_pixel(pipeline::DistortionPipeline{Nothing, Tuple{Nothing, Nothing}, Tuple{Nothing, Nothing}, P, Q}, focal::StaticVector{N}, ::Val{N}) where {P <: Tuple{Any, Vararg{Any}}, Q, N} + T = _coordinate_float_type(focal) + target = SVector{N, T}(focal) + return invert_tpd_stage(pipeline.tpd_pre, target) +end + # Has Paper IV lookup stage, so we must iterate to invert the full pipeline. function focal_to_pixel(pipeline::DistortionPipeline, focal::StaticVector{N}, ::Val{N}) where {N} T = _coordinate_float_type(focal) diff --git a/src/linear.jl b/src/linear.jl index ab80405..005bd2f 100644 --- a/src/linear.jl +++ b/src/linear.jl @@ -3,7 +3,8 @@ Linear pixel-to-intermediate and intermediate-to-pixel transforms. Paper I (Greisen & Calabretta 2002), Section 2. -The intermediate world coordinate xᵢ is related to pixel coordinate pⱼ by +Without a sequent distortion, the intermediate world coordinate xᵢ is related +to pixel coordinate pⱼ by xᵢ = Σⱼ CDᵢⱼ (pⱼ − CRPIXⱼ) @@ -13,6 +14,12 @@ units of their corresponding `CRVALᵢ` (typically degrees for celestial axes). The inverse is pⱼ = CRPIXⱼ + Σᵢ (CD⁻¹)ⱼᵢ xᵢ + +With TPV or a Paper IV sequent distortion D, WCSLIB instead evaluates + + q = PC (p − CRPIX), x = CDELT ⊙ D(q) + +so the PC/CDELT decomposition must be preserved. """ """ @@ -31,9 +38,27 @@ function pixel_to_intermediate(wcs::WCSTransform{N}, pixel::AbstractVector) wher T = _coordinate_float_type(pixel) focal = pixel_to_focal(wcs.pipeline, pixel, Val(N)) delta = SVector{N,T}(ntuple(i -> T(focal[i]) - T(wcs.crpix[i]), N)) - # Convert CD to T so the matrix-vector product preserves input precision. - cd_T = SMatrix{N, N, T}(wcs.cd) - return cd_T * delta + + # The ordinary affine path can use the combined CD matrix directly. + if !(wcs.projection isa TPV) && !has_sequent_distortion(wcs.pipeline) + cd_T = SMatrix{N, N, T}(wcs.cd) + return cd_T * delta + end + + # Sequent distortions operate after PC and before the per-axis CDELT scale. + pc_T = SMatrix{N, N, T}(wcs.pc) + coord = pc_T * delta + if wcs.projection isa TPV + lon = wcs.lon_axis + lat = wcs.lat_axis + x = _evaluate_tpv_polynomial(wcs.projection.xcoeff, coord[lon], coord[lat]) + y = _evaluate_tpv_polynomial(wcs.projection.ycoeff, coord[lat], coord[lon]) + coord = SVector{N, T}(ntuple(i -> i == lon ? T(x) : i == lat ? T(y) : coord[i], N)) + end + coord = apply_sequent_distortion(wcs.pipeline, coord) + + # Apply the final scale in the input coordinate's floating-point type. + return coord .* SVector{N, T}(wcs.cdelt) end @@ -47,11 +72,26 @@ if the matrix is singular. """ function intermediate_to_pixel(wcs::WCSTransform{N}, intermediate::StaticVector{N}) where {N} T = _coordinate_float_type(intermediate) - # Undo the linear matrix to recover focal/image-plane pixel coordinates. - # Convert CD and CRPIX to T so that the output preserves the input precision. - cd_T = SMatrix{N, N, T}(wcs.cd) cpix_T = SVector{N, T}(wcs.crpix) - focal = cpix_T .+ (cd_T \ intermediate) + + # The ordinary affine path can invert the combined CD matrix directly. + if !(wcs.projection isa TPV) && !has_sequent_distortion(wcs.pipeline) + cd_T = SMatrix{N, N, T}(wcs.cd) + focal = cpix_T .+ (cd_T \ intermediate) + return focal_to_pixel(wcs.pipeline, focal, Val(N)) + end + + # Undo CDELT, sequent TPD, and TPV before solving the PC matrix. + coord = intermediate ./ SVector{N, T}(wcs.cdelt) + coord = invert_sequent_distortion(wcs.pipeline, coord) + if wcs.projection isa TPV + lon = wcs.lon_axis + lat = wcs.lat_axis + x, y = _tpv_inverse(wcs.projection.xcoeff, wcs.projection.ycoeff, coord[lon], coord[lat]) + coord = SVector{N, T}(ntuple(i -> i == lon ? T(x) : i == lat ? T(y) : coord[i], N)) + end + pc_T = SMatrix{N, N, T}(wcs.pc) + focal = cpix_T .+ (pc_T \ coord) return focal_to_pixel(wcs.pipeline, focal, Val(N)) end diff --git a/src/parsing.jl b/src/parsing.jl index 8f8b85c..48abf4d 100644 --- a/src/parsing.jl +++ b/src/parsing.jl @@ -36,6 +36,9 @@ array indices the values are numerically identical; no offset is required. - `cd` – combined CD matrix (naxis × naxis), where `cd[i,j] = CDELT_i * PC_i_j` (or the explicit `CD_i_j` value). Units match `cunit`. +- `pc` – linear matrix applied before a sequent distortion. +- `cdelt` – per-axis scale applied after a sequent distortion. Explicit + CD headers use `pc = cd` and unit `cdelt`. - `ctype` – FITS `CTYPEi` strings, length `naxis`. - `cunit` – FITS `CUNITi` strings (empty string means degrees for celestial axes), length `naxis`. @@ -49,7 +52,7 @@ array indices the values are numerically identical; no offset is required. Paper II. Precomputed during construction. - `projection` – spherical projection for the celestial axes, or `nothing` for purely linear WCS. -- `pipeline` – pre-linear pixel/focal-plane distortion pipeline. +- `pipeline` – distortion stages surrounding the linear transform. - `aux` – resolved auxiliary WCS data, or `NoAuxiliaryWCSData`. - `lon_axis` – 1-based index of the longitude axis; 0 if no celestial axes. - `lat_axis` – 1-based index of the latitude axis; 0 if no celestial axes. @@ -65,6 +68,8 @@ struct WCSTransform{N, L, P <: Union{Nothing, AbstractProjection}, D <: Abstract crpix::SVector{N, Float64} crval::SVector{N, Float64} cd::SMatrix{N, N, Float64, L} # naxis x naxis + pc::SMatrix{N, N, Float64, L} + cdelt::SVector{N, Float64} ctype::Vector{String} cunit::Vector{String} lonpole::Float64 # degrees; phi_p (Paper II) @@ -87,6 +92,19 @@ struct WCSTransform{N, L, P <: Union{Nothing, AbstractProjection}, D <: Abstract unit_scaling::SVector{N, Float64} end +function WCSTransform(naxis::Int, crpix::SVector{N, Float64}, crval::SVector{N, Float64}, + cd::SMatrix{N, N, Float64, L}, ctype::Vector{String}, cunit::Vector{String}, + lonpole::Float64, latpole::Float64, alpha_p::Float64, delta_p::Float64, + projection::P, pipeline::D, aux::A, lon_axis::Int, lat_axis::Int, obs::O, + radesys::String, equinox::Float64, wcsname::String, preserve_units::Bool, + unit_scaling::SVector{N, Float64}) where {N, L, P, D, A, O} + # Preserve the former positional constructor by treating its CD as an explicit matrix. + return WCSTransform(naxis, crpix, crval, cd, cd, ones(SVector{N, Float64}), + ctype, cunit, lonpole, latpole, alpha_p, delta_p, + projection, pipeline, aux, lon_axis, lat_axis, obs, + radesys, equinox, wcsname, preserve_units, unit_scaling) +end + # ────────────────────────────────────────────────────────────────────────────── # CTYPE parsing helpers # ────────────────────────────────────────────────────────────────────────────── @@ -182,7 +200,6 @@ function projection_from_code(code::AbstractString) c == "HPX" && return HPX() c == "XPH" && return XPH() c == "TPV" && return TPV() - c == "TPD" && return TPV() return UnknownProjection(c) end @@ -190,6 +207,9 @@ function projection_from_header(code::AbstractString, header::AbstractDict, lon_axis::Int, lat_axis::Int, alt::Char) c = uppercase(strip(code)) + # TPD is declared through CPDIS/CQDIS rather than as a CTYPE projection. + c == "TPD" && throw(ArgumentError("TPD is a Paper IV distortion function, not a celestial projection code")) + if c == "AZP" alt_str = alt == ' ' ? "" : string(alt) mu = Float64(get(header, "PV$(lat_axis)_1$(alt_str)", 0.0)) @@ -286,16 +306,12 @@ function projection_from_header(code::AbstractString, header::AbstractDict, return HPX(H, K) end - # TPV/TPD: TAN + sequent polynomial distortion encoded in PVi_m keywords. + # TPV: TAN + sequent polynomial distortion encoded in PVi_m keywords. # Coefficients are collected from PV keywords on both celestial axes. - if c == "TPV" || c == "TPD" + if c == "TPV" alt_str = alt == ' ' ? "" : string(alt) - xcoeff_raw = _collect_tpv_coeffs(header, lon_axis, alt_str) - ycoeff_raw = _collect_tpv_coeffs(header, lat_axis, alt_str) - # If both empty, return identity TPV (= plain TAN). - xcoeff = isempty(xcoeff_raw) ? Float64[0.0, 1.0] : xcoeff_raw - ycoeff = isempty(ycoeff_raw) ? Float64[0.0, 0.0, 1.0] : ycoeff_raw - return TPV(xcoeff, ycoeff) + return TPV(_collect_tpv_coeffs(header, lon_axis, alt_str), + _collect_tpv_coeffs(header, lat_axis, alt_str)) end # Other supported projections currently have no parsed PV parameters. @@ -305,14 +321,14 @@ end """ _collect_tpv_coeffs(header, axis, alt_str) -> Vector{Float64} -Collect TPV/TPD polynomial coefficients from `PV{axis}_m` keywords for -``m = 0..59`` (the TPD coefficient range). Gaps are zero-filled so the +Collect TPV polynomial coefficients from `PV{axis}_m` keywords for +``m = 0..39``. Gaps are zero-filled so the returned vector is indexed directly by ``m`` (i.e. `result[m+1] = PVm`). Returns an empty vector if no PV keywords are found on this axis. """ function _collect_tpv_coeffs(header::AbstractDict, axis::Int, alt_str::AbstractString) coeff = Float64[] - for m in 0:59 + for m in 0:39 key = "PV$(axis)_$(m)$(alt_str)" if haskey(header, key) while length(coeff) < m @@ -321,6 +337,11 @@ function _collect_tpv_coeffs(header::AbstractDict, axis::Int, alt_str::AbstractS push!(coeff, Float64(header[key])) end end + + # TPV is limited to seventh degree even though the TPD basis has 60 terms. + invalid = Regex("^PV$(axis)_([4-9][0-9]|[1-9][0-9]{2,})$(alt_str)\$") + any(key -> key isa AbstractString && occursin(invalid, uppercase(String(key))), keys(header)) && + throw(ArgumentError("TPV coefficients must use parameter numbers 0 through 39")) return coeff end @@ -838,10 +859,8 @@ The primary WCS (no suffix) uses `alt=' '`. Alternate versions use `alt='A'`, Each call to `WCS` returns exactly one `WCSTransform` for the selected alternate. To extract all alternates, call `WCS` once per character. -`minerr` is reserved for Paper IV distortion support: when auxiliary lookup -tables are implemented, distortion components whose declared error estimate is -below `minerr` may be skipped. It is currently passed through to auxiliary-data -resolver methods but is not used by the core header-only parser. +`minerr` skips Paper IV lookup or TPD distortion components whose declared +error estimate is below the threshold. `preserve_units` controls whether `pixel_to_world` returns world coordinates in canonical units (`false`, default) or in the original header CUNIT (`true`). @@ -956,13 +975,16 @@ function WCS(header::AbstractDict; fobj = nothing, alt::Char = ' ', minerr::Real projection = isempty(proj_code) ? nothing : projection_from_header(proj_code, header, lon_axis, lat_axis, alt) - # ── Parse optional SIP distortion before building transform output ──────── - # TPV/TPD takes precedence over SIP (SCAMP convention). - if proj_code in ("TPV", "TPD") + # ── Parse prior and sequent distortions before building output ────────────── + # TPV takes precedence over SIP in legacy SCAMP headers. + if proj_code == "TPV" _remove_sip_keywords!(header, alt_str) end sip = parse_sip_distortion(header, crpix, naxis, alt) - pipeline = distortion_pipeline(sip, aux) + tpd_pre, tpd_seq = parse_tpd_distortions(header, naxis, alt, minerr) + proj_code == "TPV" && any(!isnothing, tpd_seq) && + throw(ArgumentError("TPV cannot be combined with a CQDIS sequent distortion")) + pipeline = distortion_pipeline(sip, aux, tpd_pre, tpd_seq) # ── Normalize crval to canonical units before building spectral specs ───── # TAB axes are skipped: their crval values are lookup-table indices, not @@ -1014,6 +1036,21 @@ function WCS(header::AbstractDict; fobj = nothing, alt::Char = ' ', minerr::Real end end + # Preserve PC and CDELT separately only where a sequent distortion needs them. + has_cd_keywords = has_indexed_wcs_keyword(header, "CD", naxis, alt_str) + needs_linear_split = projection isa TPV || has_sequent_distortion(pipeline) + linear_cdelt = has_cd_keywords || !needs_linear_split ? ones(Float64, naxis) : cdelt .* unit_scaling_vec + pc = if has_cd_keywords || !needs_linear_split + copy(cd) + else + any(iszero, linear_cdelt) && throw(ArgumentError("CDELT values must be nonzero with a sequent distortion")) + pre = copy(cd) + for i in 1:naxis + pre[i, :] ./= linear_cdelt[i] + end + pre + end + # ── Celestial pole parameters ───────────────────────────────────────────── lonpole_raw = get(header, "LONPOLE$(alt_str)", get(header, "LONPOLE", nothing)) latpole_raw = get(header, "LATPOLE$(alt_str)", get(header, "LATPOLE", nothing)) @@ -1080,6 +1117,8 @@ function WCS(header::AbstractDict; fobj = nothing, alt::Char = ' ', minerr::Real SVector{naxis, Float64}(crpix), SVector{naxis, Float64}(crval), SMatrix{naxis, naxis, Float64, naxis * naxis}(cd), + SMatrix{naxis, naxis, Float64, naxis * naxis}(pc), + SVector{naxis, Float64}(linear_cdelt), ctype, cunit, lonpole, latpole, alpha_p, delta_p, projection, pipeline, aux, lon_axis, lat_axis, obs, diff --git a/src/projections.jl b/src/projections.jl index 37df172..dc45034 100644 --- a/src/projections.jl +++ b/src/projections.jl @@ -326,15 +326,16 @@ TAN projection with sequent polynomial distortion (SCAMP convention). TPV is not a fundamental projection; it is a TAN (gnomonic) projection combined with a sequent polynomial distortion whose coefficients are stored in FITS -`PVi_m` keywords. The polynomial operates on intermediate world coordinates -(x, y) in degrees **after** the CD matrix, using the TPD coefficient indexing -convention. +`PVi_m` keywords. WCSLIB applies the polynomial after PC and before CDELT, +using the TPD coefficient indexing convention. -Coefficients are stored in **direct** form (wcslib convention): the polynomial -returns the corrected coordinate directly. The identity polynomial is -``xcoeff = [0, 1]`` (i.e. ``x' = x``) and ``ycoeff = [0, 0, 1]`` (``y' = y``). +Coefficients are stored in **direct** form (WCSLIB convention): the polynomial +returns the corrected coordinate directly. Since TPD variables are local to +the attached axis, identity is encoded by coefficient 1 on term 1 for both +axes. -The CTYPE projection codes `TPV` and `TPD` both map to this type. +Only the CTYPE projection code `TPV` maps to this type. TPD is a Paper IV +distortion function rather than a celestial projection code. # Fields - `xcoeff` – ``PV_{lon\\_axis,m}`` coefficients (direct form) @@ -345,8 +346,13 @@ struct TPV <: AbstractProjection ycoeff::Vector{Float64} end -# Default is identity polynomial (= plain TAN). x' = 1·x, y' = 1·y. -TPV() = TPV(Float64[0.0, 1.0], Float64[0.0, 0.0, 1.0]) +TPV(xcoeff::AbstractVector, ycoeff::AbstractVector) = TPV(Float64.(xcoeff), Float64.(ycoeff)) + +# Missing PVi_m coefficients default to zero in WCSLIB. +TPV() = TPV(Float64[], Float64[]) + +_spherical_projection(proj::AbstractProjection) = proj +_spherical_projection(::TPV) = TAN() # ── Unknown / deferred ──────────────────────────────────────────────────────── @@ -479,7 +485,7 @@ function _tpv_inverse(xcoeff::AbstractVector, ycoeff::AbstractVector, @inbounds for k in 1:max_iter x_corr = _evaluate_tpv_polynomial(xcoeff, u, v) - y_corr = _evaluate_tpv_polynomial(ycoeff, u, v) + y_corr = _evaluate_tpv_polynomial(ycoeff, v, u) du = x_corr - xt dv = y_corr - yt u -= du @@ -773,7 +779,7 @@ delegate to the TAN (gnomonic) inverse. """ function intermediate_to_native(proj::TPV, x::Real, y::Real) x_corr = _evaluate_tpv_polynomial(proj.xcoeff, x, y) - y_corr = _evaluate_tpv_polynomial(proj.ycoeff, x, y) + y_corr = _evaluate_tpv_polynomial(proj.ycoeff, y, x) return intermediate_to_native(TAN(), x_corr, y_corr) end diff --git a/src/transforms.jl b/src/transforms.jl index 0359e3d..4bf07cd 100644 --- a/src/transforms.jl +++ b/src/transforms.jl @@ -228,7 +228,7 @@ function pixel_to_world(wcs::WCSTransform, pixel::AbstractVector) x_lat = x[lat_idx] # degrees # Deproject: (x, y) -> (phi, theta) in radians - phi, theta = intermediate_to_native(wcs.projection, x_lon, x_lat) + phi, theta = intermediate_to_native(_spherical_projection(wcs.projection), x_lon, x_lat) # Spherical rotation: (phi, theta) -> (alpha, delta) in radians # alpha_p, delta_p are the celestial coords of the native north pole. @@ -315,7 +315,7 @@ function world_to_pixel(wcs::WCSTransform{N}, world::AbstractVector) where {N} phi, theta = celestial_to_native(alpha, delta, alpha_p, delta_p, phi_p) # Re-project: (phi, theta) -> (x, y) in degrees - x_lon, x_lat = native_to_intermediate(wcs.projection, phi, theta) + x_lon, x_lat = native_to_intermediate(_spherical_projection(wcs.projection), phi, theta) # Build intermediate coordinate vector # Non-celestial axes: x_i = world_i - crval_i (trivial linear axes) diff --git a/test/runtests.jl b/test/runtests.jl index 3025d87..a7a3ce8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2917,7 +2917,7 @@ end @testset "TPV inverse" begin xcoeff = Float64[0.0, 1.0] - ycoeff = Float64[0.0, 0.0, 1.0] + ycoeff = Float64[0.0, 1.0] # Identity: inverse recovers target. u, v = FITSWCS._tpv_inverse(xcoeff, ycoeff, 5.0, -3.0) @test u ≈ 5.0 @@ -2927,7 +2927,7 @@ end xc = zeros(Float64, 8) xc[2] = 1.0 # m=1: x xc[8] = 1e-5 # m=7: x³ - yc = Float64[0.0, 0.0, 1.0] + yc = Float64[0.0, 1.0] target = 2.0 + 1e-5 * 8.0 # 2 + 1e-5*8 u, v = FITSWCS._tpv_inverse(xc, yc, target, 0.0) @test u ≈ 2.0 atol=1e-10 @@ -2942,7 +2942,7 @@ end @testset "TPV projection round-trip" begin # Identity TPV matches TAN. - tpv0 = TPV() + tpv0 = TPV(Float64[0.0, 1.0], Float64[0.0, 1.0]) @test FITSWCS.native_theta0(tpv0) == 90.0 @test FITSWCS.native_phi0(tpv0) == 0.0 @@ -2954,7 +2954,7 @@ end # Non-trivial TPV: x² term on x, y² term on y. xc = Float64[0.0, 1.0, 0.0, 0.0, 5e-6] # x' = x + 5e-6·x² - yc = Float64[0.0, 0.0, 1.0, 0.0, 0.0, 0.0, -3e-6] # y' = y - 3e-6·y² + yc = Float64[0.0, 1.0, 0.0, 0.0, -3e-6] # y' = y - 3e-6·y² tpv1 = TPV(xc, yc) for _ in 1:5 phi = (rand() - 0.5) * 0.1 @@ -2967,7 +2967,7 @@ end end @testset "TPV header parsing" begin - # Minimal TPV header (no PV keywords → identity). + # Missing TPV coefficients default to zero, as in WCSLIB. hdr = Dict{String,Any}( "NAXIS" => 2, "CTYPE1" => "RA---TPV", "CTYPE2" => "DEC--TPV", @@ -2977,11 +2977,11 @@ end ) wcs = WCS(hdr) @test wcs.projection isa TPV - # Should behave like TAN with no extra distortion. pix = [256.0, 256.0] world = pixel_to_world(wcs, pix) @test world[1] ≈ 83.8221 atol=1e-4 @test world[2] ≈ -5.3911 atol=1e-4 + @test pixel_to_world(wcs, [300.0, 200.0]) ≈ world # TPV with PV coefficients. hdr2 = Dict{String,Any}( @@ -2992,8 +2992,8 @@ end "CDELT1" => -0.05, "CDELT2" => 0.05, "PV1_1" => 1.0, # identity x term "PV1_4" => 1e-6, # x² term on axis 1 - "PV2_2" => 1.0, # identity y term - "PV2_6" => -2e-6, # y² term on axis 2 + "PV2_1" => 1.0, # identity y term in axis-local coordinates + "PV2_4" => -2e-6, # y² term on axis 2 ) wcs2 = WCS(hdr2) t = wcs2.projection @@ -3001,8 +3001,8 @@ end @test t.xcoeff[1] == 0.0 # m=0 (constant), zero-filled gap @test t.xcoeff[2] == 1.0 # m=1 (x) @test t.xcoeff[5] == 1e-6 # m=4 (x²) - @test t.ycoeff[3] == 1.0 # m=2 (y) - @test t.ycoeff[7] == -2e-6 # m=6 (y²) + @test t.ycoeff[2] == 1.0 # m=1 (attached-axis coordinate) + @test t.ycoeff[5] == -2e-6 # m=4 (attached-axis coordinate squared) # Round-trip at reference pixel. pix2 = [128.0, 96.0] @@ -3010,7 +3010,7 @@ end pix2b = world_to_pixel(wcs2, world2) @test pix2b ≈ pix2 atol=1e-10 - # TPD CTYPE also works. + # TPD is a Paper IV distortion function, not a CTYPE projection. hdr3 = Dict{String,Any}( "NAXIS" => 2, "CTYPE1" => "RA---TPD", "CTYPE2" => "DEC--TPD", @@ -3018,8 +3018,7 @@ end "CRVAL1" => 45.0, "CRVAL2" => 30.0, "CDELT1" => -0.01, "CDELT2" => 0.01, ) - wcs3 = WCS(hdr3) - @test wcs3.projection isa TPV + @test_throws "TPD is a Paper IV distortion function" WCS(hdr3) end @testset "SCAMP compatibility" begin @@ -3032,7 +3031,7 @@ end "CDELT1" => -2.7778e-4, "CDELT2" => 2.7778e-4, "PV1_1" => 1.0, # TPV coefficients hiding as TAN "PV1_5" => 1e-6, # j=5 signals SCAMP - "PV2_2" => 1.0, + "PV2_1" => 1.0, ) wcs_pre = WCS(hdr_pre) @test wcs_pre.projection isa TPV @@ -3070,7 +3069,7 @@ end @testset "TPV Float32 type stability" begin # Polynomial eval in Float32. xc32 = Float32[0.0, 1.0] - yc32 = Float32[0.0, 0.0, 1.0] + yc32 = Float32[0.0, 1.0] @test @inferred(FITSWCS._evaluate_tpv_polynomial(xc32, 1.0f0, 2.0f0)) isa Float32 # Inverse in Float32. @@ -3095,6 +3094,7 @@ end "CRPIX1" => 256.0, "CRPIX2" => 256.0, "CRVAL1" => 83.8221, "CRVAL2" => -5.3911, "CDELT1" => -2.7778e-4, "CDELT2" => 2.7778e-4, + "PV1_1" => 1.0, "PV2_1" => 1.0, ) wcs32 = WCS(hdr32) pix32 = Float32[256.0, 256.0] @@ -3118,8 +3118,8 @@ end "CDELT1" => -2.7778e-4, "CDELT2" => 2.7778e-4, "PV1_1" => 1.0, # identity x "PV1_4" => 1e-7, # small x² correction - "PV2_2" => 1.0, # identity y - "PV2_6" => -1e-7, # small y² correction + "PV2_1" => 1.0, # identity y in axis-local coordinates + "PV2_4" => -1e-7, # small y² correction ) wcs = WCS(hdr) @test wcs.projection isa TPV @@ -3210,3 +3210,6 @@ include("regression_astropy_values.jl") # WCS slicing tests. include("slicing_tests.jl") + +# TPV ordering and Paper IV TPD regression tests. +include("tpv_tpd_tests.jl") diff --git a/test/slicing_tests.jl b/test/slicing_tests.jl index a786ca6..787a7ac 100644 --- a/test/slicing_tests.jl +++ b/test/slicing_tests.jl @@ -426,7 +426,7 @@ end "CRVAL1" => 83.8221, "CRVAL2" => -5.3911, "CDELT1" => -2.7778e-4, "CDELT2" => 2.7778e-4, "PV1_0" => 0.0, "PV1_1" => 1.0, - "PV2_0" => 0.0, "PV2_1" => 0.0, "PV2_2" => 1.0, + "PV2_0" => 0.0, "PV2_1" => 1.0, ) wcs = WCS(hdr) diff --git a/test/tpv_tpd_tests.jl b/test/tpv_tpd_tests.jl new file mode 100644 index 0000000..3d0ddb6 --- /dev/null +++ b/test/tpv_tpd_tests.jl @@ -0,0 +1,191 @@ +@testset "TPV sequent distortion ordering" begin + # A nontrivial PC/CDELT split must match WCSLIB's PC -> TPV -> CDELT pipeline. + header = Dict{String,Any}( + "NAXIS" => 2, + "CTYPE1" => "RA---TPV", "CTYPE2" => "DEC--TPV", + "CRPIX1" => 100.0, "CRPIX2" => 200.0, + "CRVAL1" => 123.4, "CRVAL2" => -22.5, + "CDELT1" => -0.02, "CDELT2" => 0.035, + "PC1_1" => 0.8, "PC1_2" => 0.3, + "PC2_1" => -0.2, "PC2_2" => 1.1, + "PV1_0" => 0.01, "PV1_1" => 1.02, "PV1_2" => -0.03, + "PV1_4" => 0.0007, "PV1_5" => -0.0002, "PV1_6" => 0.0005, + "PV2_0" => -0.02, "PV2_1" => 0.5, "PV2_2" => 0.5, + "PV2_4" => 0.0004, "PV2_5" => -0.0003, "PV2_6" => 0.0004, + ) + wcs = WCS(header) + pixel = [112.3, 176.4] + expected_intermediate = [-0.08205336719999999, -0.43751199240000005] + expected_world = [123.31090434225322, -22.937478626425648] + @test pixel_to_intermediate(wcs, pixel) ≈ expected_intermediate atol=2e-15 + @test pixel_to_world(wcs, pixel) ≈ expected_world atol=2e-14 + @test world_to_pixel(wcs, expected_world) ≈ pixel atol=1e-11 + + # Explicit CD is interpreted as PC with unit CDELT, as WCSLIB requires for TPV. + cd_header = copy(header) + for key in ("CDELT1", "CDELT2", "PC1_1", "PC1_2", "PC2_1", "PC2_2") + delete!(cd_header, key) + end + cd_header["CD1_1"] = -0.016 + cd_header["CD1_2"] = -0.006 + cd_header["CD2_1"] = -0.007 + cd_header["CD2_2"] = 0.0385 + cd_wcs = WCS(cd_header) + @test pixel_to_intermediate(cd_wcs, pixel) ≈ [-0.01597713451499998, -0.5445694821799998] atol=2e-15 + @test pixel_to_world(cd_wcs, pixel) ≈ [123.38263812267824, -23.0445521374377] atol=2e-14 +end + +@testset "TPV axis-local coefficients" begin + # Standard identity coefficients use term 1 on both attached axes. + header = Dict{String,Any}( + "NAXIS" => 2, + "CTYPE1" => "RA---TPV", "CTYPE2" => "DEC--TPV", + "CRPIX1" => 0.0, "CRPIX2" => 0.0, + "CRVAL1" => 0.0, "CRVAL2" => 0.0, + "CDELT1" => 1.0, "CDELT2" => 1.0, + "PV1_1" => 1.0, "PV2_1" => 1.0, + ) + wcs = WCS(header) + @test pixel_to_intermediate(wcs, [2.0, 3.0]) == [2.0, 3.0] + @test pixel_to_world(wcs, [2.0, 3.0]) ≈ [1.9991882802315373, 2.9954418975929276] atol=2e-14 + @test world_to_pixel(wcs, pixel_to_world(wcs, [2.0, 3.0])) ≈ [2.0, 3.0] atol=2e-14 + + # TPV has only the first 40 TPD terms. + header["PV1_40"] = 1.0 + @test_throws "parameter numbers 0 through 39" WCS(header) +end + +@testset "Paper IV prior TPD" begin + # Prior additive TPD must run before the linear WCS matrix and invert cleanly. + header = Dict{String,Any}( + "NAXIS" => 2, "CTYPE1" => "LINEAR", "CTYPE2" => "LINEAR", + "CRPIX1" => 0.0, "CRPIX2" => 0.0, + "CRVAL1" => 0.0, "CRVAL2" => 0.0, + "CDELT1" => 1.0, "CDELT2" => 1.0, + "CPDIS1" => "TPD", "CPDIS2" => "TPD", + "DP1.NAXES" => 2, "DP1.AXIS.1" => 1, "DP1.AXIS.2" => 2, + "DP1.TPD.FWD.4" => 0.01, + "DP2.NAXES" => 2, "DP2.AXIS.1" => 2, "DP2.AXIS.2" => 1, + "DP2.TPD.FWD.4" => -0.02, + ) + wcs = WCS(header) + @test pixel_to_world(wcs, [2.0, 3.0]) ≈ [2.04, 2.82] atol=2e-15 + @test world_to_pixel(wcs, [2.04, 2.82]) ≈ [2.0, 3.0] atol=1e-13 + + # This nonzero-CRPIX case verifies WCSLIB's prior TPD -> CRPIX -> PC -> CDELT order. + rich = Dict{String,Any}( + "NAXIS" => 2, "CTYPE1" => "LINEAR", "CTYPE2" => "LINEAR", + "CRPIX1" => 10.0, "CRPIX2" => 20.0, + "CRVAL1" => 5.0, "CRVAL2" => -7.0, + "CDELT1" => -0.02, "CDELT2" => 0.035, + "PC1_1" => 0.8, "PC1_2" => 0.3, + "PC2_1" => -0.2, "PC2_2" => 1.1, + "CPDIS1" => "TPD", "CPDIS2" => "TPD", + "DP1.NAXES" => 2, "DP1.AXIS.1" => 1, "DP1.AXIS.2" => 2, + "DP1.OFFSET.1" => 1.0, "DP1.OFFSET.2" => -2.0, + "DP1.SCALE.1" => 0.5, "DP1.SCALE.2" => 2.0, + "DP1.TPD.FWD.0" => 0.01, "DP1.TPD.FWD.1" => 1.02, "DP1.TPD.FWD.2" => -0.03, + "DP1.TPD.FWD.4" => 0.0007, "DP1.TPD.FWD.5" => -0.0002, "DP1.TPD.FWD.6" => 0.0005, + "DP2.NAXES" => 2, "DP2.AXIS.1" => 2, "DP2.AXIS.2" => 1, + "DP2.TPD.FWD.4" => -0.002, "DP2.TPD.FWD.5" => 0.0003, + ) + rich_wcs = WCS(rich) + pixel = [12.3, -3.4] + expected_world = [5.00963112, -7.95957814925] + @test pixel_to_intermediate(rich_wcs, pixel) ≈ [0.009631119999999993, -0.9595781492500003] atol=2e-15 + @test pixel_to_world(rich_wcs, pixel) ≈ expected_world atol=2e-14 + @test world_to_pixel(rich_wcs, expected_world) ≈ pixel atol=2e-12 +end + +@testset "Paper IV sequent TPD" begin + # Sequent TPD exercises axis maps, normalization, direct output, auxiliary variables, and reverse seeds. + header = Dict{String,Any}( + "NAXIS" => 2, "CTYPE1" => "LINEAR", "CTYPE2" => "LINEAR", + "CRPIX1" => 10.0, "CRPIX2" => 20.0, + "CRVAL1" => 5.0, "CRVAL2" => -7.0, + "CDELT1" => -0.02, "CDELT2" => 0.035, + "PC1_1" => 0.8, "PC1_2" => 0.3, + "PC2_1" => -0.2, "PC2_2" => 1.1, + "CQDIS1" => "TPD", "CQDIS2" => "TPD", + "DQ1.NAXES" => 2, "DQ1.AXIS.1" => 1, "DQ1.AXIS.2" => 2, + "DQ1.OFFSET.1" => 1.0, "DQ1.OFFSET.2" => -2.0, + "DQ1.SCALE.1" => 0.5, "DQ1.SCALE.2" => 2.0, "DQ1.DOCORR" => 0, + "DQ1.AUX.1.COEFF.0" => 0.1, "DQ1.AUX.1.COEFF.1" => 1.1, "DQ1.AUX.1.COEFF.2" => -0.2, + "DQ1.AUX.2.COEFF.0" => -0.3, "DQ1.AUX.2.COEFF.1" => 0.4, "DQ1.AUX.2.COEFF.2" => 0.9, + "DQ1.TPD.FWD.0" => 0.01, "DQ1.TPD.FWD.1" => 1.02, "DQ1.TPD.FWD.2" => -0.03, + "DQ1.TPD.FWD.4" => 0.0007, "DQ1.TPD.FWD.5" => -0.0002, "DQ1.TPD.FWD.6" => 0.0005, + "DQ1.TPD.REV.1" => 1.0, + "DQ2.NAXES" => 2, "DQ2.AXIS.1" => 2, "DQ2.AXIS.2" => 1, + "DQ2.TPD.FWD.4" => -0.002, "DQ2.TPD.FWD.5" => 0.0003, + ) + wcs = WCS(header) + pixel = [12.3, -3.4] + expected_world = [4.820512437282, -7.963625782] + @test pixel_to_intermediate(wcs, pixel) ≈ [-0.17948756271800004, -0.9636257820000002] atol=2e-15 + @test pixel_to_world(wcs, pixel) ≈ expected_world atol=2e-14 + @test world_to_pixel(wcs, expected_world) ≈ pixel atol=2e-12 + + # minerr suppresses each complete TPD axis before constructing the pipeline. + header["CQERR1"] = 0.01 + header["CQERR2"] = 0.01 + skipped = WCS(header; minerr = 0.02) + @test pixel_to_world(skipped, pixel) ≈ [5.1036, -7.917] atol=2e-15 +end + +@testset "One-dimensional ninth-degree TPD" begin + # One-dimensional TPD keeps pure attached-axis powers and ignores radial terms. + header = Dict{String,Any}( + "NAXIS" => 2, "CTYPE1" => "LINEAR", "CTYPE2" => "LINEAR", + "CRPIX1" => 0.0, "CRPIX2" => 0.0, + "CRVAL1" => 0.0, "CRVAL2" => 0.0, + "CDELT1" => 1.0, "CDELT2" => 1.0, + "CPDIS1" => "TPD", "CPDIS2" => "TPD", + "DP1.NAXES" => 1, "DP1.AXIS.1" => 1, + "DP1.TPD.FWD.49" => 1e-9, "DP1.TPD.FWD.59" => 10.0, + "DP2.NAXES" => 1, "DP2.AXIS.1" => 2, + ) + @test pixel_to_world(WCS(header), [2.0, 3.0]) ≈ [2.000000512, 3.0] atol=2e-15 + + # Two-dimensional radial terms use both mapped coordinates. + radial = copy(header) + radial["DP1.NAXES"] = 2 + radial["DP1.AXIS.2"] = 2 + delete!(radial, "DP1.TPD.FWD.49") + radial["DP1.TPD.FWD.59"] = 1e-9 + @test pixel_to_world(WCS(radial), [3.0, 4.0]) ≈ [3.001953125, 4.0] atol=2e-15 +end + +@testset "FITSFiles repeated TPD parameter cards" begin + # Repeated DP cards must survive backend conversion as distinct flattened parameters. + cards = FITSFiles.Card[ + FITSFiles.Card("NAXIS", 2), + FITSFiles.Card("CTYPE1", "LINEAR"), FITSFiles.Card("CTYPE2", "LINEAR"), + FITSFiles.Card("CRPIX1", 0.0), FITSFiles.Card("CRPIX2", 0.0), + FITSFiles.Card("CRVAL1", 0.0), FITSFiles.Card("CRVAL2", 0.0), + FITSFiles.Card("CDELT1", 1.0), FITSFiles.Card("CDELT2", 1.0), + FITSFiles.Card("CPDIS1", "TPD"), FITSFiles.Card("CPDIS2", "TPD"), + FITSFiles.Card("DP1", "NAXES: 2"), FITSFiles.Card("DP1", "AXIS.1: 1"), + FITSFiles.Card("DP1", "AXIS.2: 2"), FITSFiles.Card("DP1", "TPD.FWD.4: 0.01"), + FITSFiles.Card("DP2", "NAXES: 2"), FITSFiles.Card("DP2", "AXIS.1: 2"), + FITSFiles.Card("DP2", "AXIS.2: 1"), FITSFiles.Card("DP2", "TPD.FWD.4: -0.02"), + ] + @test pixel_to_world(WCS(cards), [2.0, 3.0]) ≈ [2.04, 2.82] atol=2e-15 +end + +@testset "FITSIO repeated TPD parameter cards" begin + # FITSIO's repeated string-card representation must retain every TPD field. + extension = Base.get_extension(FITSWCS, :FITSWCSFITSIOExt) + header = Dict{String,Any}() + for value in ("NAXES: 2", "AXIS.1: 1", "AXIS.2: 2", "DOCORR: 0", + "OFFSET.1: 3", "SCALE.2: 4", "AUX.1.COEFF.0: 5", + "TPD.FWD.4: 0.01", "TPD.REV.4: -0.01") + extension._set_fitsio_header_card!(header, "DP1", value) + end + @test header["DP1.NAXES"] == 2.0 + @test header["DP1.DOCORR"] == 0.0 + @test header["DP1.OFFSET.1"] == 3.0 + @test header["DP1.SCALE.2"] == 4.0 + @test header["DP1.AUX.1.COEFF.0"] == 5.0 + @test header["DP1.TPD.FWD.4"] == 0.01 + @test header["DP1.TPD.REV.4"] == -0.01 +end