Skip to content

Commit

Permalink
Merge pull request #80 from observerly/refactor/transform/SIP2DForwar…
Browse files Browse the repository at this point in the history 
…dParameters

refactor!: amend to SIP2DForwardParameters usage across @observerly/skysolve
  • Loading branch information
@michealroberts
michealroberts authored Nov 19, 2024
2 parents 7d31d0f + ae1877f commit 276c9c1
Show file tree
Hide file tree
Showing 5 changed files with 35 additions and 29 deletions.
19 changes: 10 additions & 9 deletions pkg/solver/solver.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -517,19 +517,20 @@ func (ps *PlateSolver) solveForAffineParameters(

/*****************************************************************************************************************/

// solveForSIPParameters fits higher-order SIP polynomials to the non-linear residuals after the affine transformation.
// solveForForwardSIPParameters fits higher-order SIP polynomials to the non-linear residuals after the
// affine transformation.
//
// SIP’s Purpose for Non-linear Distortions: SIP is specifically designed to correct non-linear distortions.
// Terms where p + q <= 1 represent linear transformations, which are unnecessary in SIP since they’re covered
// by the affine transformations.
func (ps *PlateSolver) solveForSIPParameters(
func (ps *PlateSolver) solveForForwardSIPParameters(
aRA [][]float64,
aDec [][]float64,
bRA []float64,
bDec []float64,
n int,
sipOrder int,
) (*transform.SIP2DParameters, error) {
) (*transform.SIP2DForwardParameters, error) {
// Calculate the number of terms in the SIP polynomial:
numTerms := (sipOrder + 1) * (sipOrder + 2) / 2

Expand Down Expand Up @@ -640,7 +641,7 @@ func (ps *PlateSolver) solveForSIPParameters(
bPowerMap[term] = sipParamsDec[idx]
}

sipParams := transform.SIP2DParameters{
sipParams := transform.SIP2DForwardParameters{
AOrder: sipOrder,
BOrder: sipOrder,
APower: aPowerMap,
Expand Down Expand Up @@ -727,14 +728,14 @@ func (ps *PlateSolver) Solve(tolerance geometry.InvariantFeatureTolerance, sipOr
}

// Create the SIP parameters by solving the SIP polynomials from the residuals:
sipParams, err := ps.solveForSIPParameters(A_SIP_RA, A_SIP_Dec, B_SIP_RA, B_SIP_Dec, n, sipOrder)
fsipParams, err := ps.solveForForwardSIPParameters(A_SIP_RA, A_SIP_Dec, B_SIP_RA, B_SIP_Dec, n, sipOrder)

if err != nil {
return nil, err
}

// Check if the SIP parameters are nil:
if sipParams == nil {
if fsipParams == nil {
return nil, errors.New("failed to compute SIP transformation matrix parameters")
}

Expand All @@ -750,9 +751,9 @@ func (ps *PlateSolver) Solve(tolerance geometry.InvariantFeatureTolerance, sipOr
x,
y,
wcs.WCSParams{
Projection: wcs.RADEC_TANSIP,
AffineParams: *affineParams,
SIPParams: *sipParams,
Projection: wcs.RADEC_TANSIP,
AffineParams: *affineParams,
SIPForwardParams: *fsipParams,
},
)

Expand Down
7 changes: 6 additions & 1 deletion pkg/transform/sip.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,11 @@ package transform
// Coordinate System (WCS) by introducing higher-order polynomial terms that account for non-linear
// optical distortions, such as those introduced by telescope optics or atmospheric effects.
// @see https://fits.gsfc.nasa.gov/registry/sip/SIP_distortion_v1_0.pdf
type SIP2DParameters struct {

/*****************************************************************************************************************/

// The forward parameters are polynomial coefficients used to map from pixel coordinates to world coordinates.
type SIP2DForwardParameters struct {
AOrder int
APower map[string]float64
BOrder int
Expand All @@ -24,6 +28,7 @@ type SIP2DParameters struct {

/*****************************************************************************************************************/

// The inverse paramaters are polynomial coefficients used to map from world coordinates to pixel coordinates.
type SIP2DInverseParameters struct {
APOrder int
APPower map[string]float64
Expand Down
4 changes: 2 additions & 2 deletions pkg/transform/sip_test.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,8 +14,8 @@ import "testing"

/*****************************************************************************************************************/

func TestSIP2DParameters(t *testing.T) {
sip := SIP2DParameters{
func TestSIP2DForwardParameters(t *testing.T) {
sip := SIP2DForwardParameters{
AOrder: 1,
BOrder: 1,
APower: map[string]float64{
Expand Down
10 changes: 5 additions & 5 deletions pkg/wcs/wcs.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ func (c CoordinateProjectionType) ToCTypes() CTypeP {
type WCSParams struct {
Projection CoordinateProjectionType // Projection type e.g., "TAN", or "TAN-SIP"
AffineParams transform.Affine2DParameters // Affine transformation parameters
SIPParams transform.SIP2DParameters // SIP transformation (distortion) coefficients, x, y to RA, Dec
SIPForwardParams transform.SIP2DForwardParameters // SIP forward transformation (distortion) coefficients, x, y to RA, Dec
SIPInverseParams transform.SIP2DInverseParameters // SIP inverse transformation (distortion) coefficients RA, Dec to x, y
}

Expand All @@ -86,7 +86,7 @@ type WCS struct {
CD2_2 float64 `hdu:"CD2_2"` // Affine transform parameter D (no default)
E float64 `hdu:"E"` // Affine translation parameter e (optional, no default)
F float64 `hdu:"F"` // Affine translation parameter f (optional, no default)
SIP transform.SIP2DParameters `` // SIP transformation (distortion) coefficients
FSIP transform.SIP2DForwardParameters `` // SIP forward transformation (distortion) coefficients
ISIP transform.SIP2DInverseParameters `` // SIP inverse transformation (distortion) coefficients
}

Expand All @@ -113,7 +113,7 @@ func NewWorldCoordinateSystem(xc float64, yc float64, params WCSParams) WCS {
CD2_2: params.AffineParams.D,
E: params.AffineParams.E,
F: params.AffineParams.F,
SIP: params.SIPParams,
FSIP: params.SIPForwardParams,
ISIP: params.SIPInverseParams,
}

Expand Down Expand Up @@ -168,7 +168,7 @@ func (wcs *WCS) PixelToEquatorialCoordinate(
B := 0.0

// Apply A polynomial corrections:
for term, coeff := range wcs.SIP.APower {
for term, coeff := range wcs.FSIP.APower {
i, j, err := parseSIPTerm(term, "A")
if err != nil {
continue
Expand All @@ -177,7 +177,7 @@ func (wcs *WCS) PixelToEquatorialCoordinate(
}

// Apply B polynomial corrections:
for term, coeff := range wcs.SIP.BPower {
for term, coeff := range wcs.FSIP.BPower {
i, j, err := parseSIPTerm(term, "B")
if err != nil {
continue
Expand Down
24 changes: 12 additions & 12 deletions pkg/wcs/wcs_test.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ func TestNewWCS(t *testing.T) {
F: 0,
}

sip := transform.SIP2DParameters{
fsip := transform.SIP2DForwardParameters{
AOrder: 1,
BOrder: 1,
APower: map[string]float64{
Expand All @@ -45,9 +45,9 @@ func TestNewWCS(t *testing.T) {
}

wcs := NewWorldCoordinateSystem(1000, 1000, WCSParams{
AffineParams: affine,
Projection: RADEC_TAN,
SIPParams: sip,
AffineParams: affine,
Projection: RADEC_TAN,
SIPForwardParams: fsip,
})

if wcs.CRPIX1 != 1000 {
Expand Down Expand Up @@ -160,7 +160,7 @@ func TestPixelToEquatorialCoordinate(t *testing.T) {
/*****************************************************************************************************************/

func TestPixelToEquatorialCoordinateWithSIPDistortionAtImageCenter(t *testing.T) {
sip := transform.SIP2DParameters{
fsip := transform.SIP2DForwardParameters{
AOrder: 3,
BOrder: 3,
APower: map[string]float64{
Expand Down Expand Up @@ -192,7 +192,7 @@ func TestPixelToEquatorialCoordinateWithSIPDistortionAtImageCenter(t *testing.T)
CD1_2: 0.0, // No rotation
CD2_1: 0.0, // No rotation
CD2_2: 0.0002777778, // 1/3600 deg/pixel
SIP: sip,
FSIP: fsip,
}

coordinate := wcs.PixelToEquatorialCoordinate(1024.0, 1024.0)
Expand All @@ -209,7 +209,7 @@ func TestPixelToEquatorialCoordinateWithSIPDistortionAtImageCenter(t *testing.T)
/*****************************************************************************************************************/

func TestPixelToEquatorialCoordinateWithSIPDistortion(t *testing.T) {
sip := transform.SIP2DParameters{
fsip := transform.SIP2DForwardParameters{
AOrder: 3,
BOrder: 3,
APower: map[string]float64{
Expand Down Expand Up @@ -241,7 +241,7 @@ func TestPixelToEquatorialCoordinateWithSIPDistortion(t *testing.T) {
CD1_2: 0.0, // No rotation
CD2_1: 0.0, // No rotation
CD2_2: 0.0002777778, // 1/3600 deg/pixel
SIP: sip,
FSIP: fsip,
}

coordinate := wcs.PixelToEquatorialCoordinate(1000.0, 1000.0)
Expand Down Expand Up @@ -314,7 +314,7 @@ func TestEquatorialCoordinateToPixel(t *testing.T) {
/*****************************************************************************************************************/

func TestEquatorialCoordinateToPixelWithSIPDistortionAtImageCenter(t *testing.T) {
sip := transform.SIP2DParameters{
fsip := transform.SIP2DForwardParameters{
AOrder: 3,
BOrder: 3,
APower: map[string]float64{
Expand Down Expand Up @@ -372,7 +372,7 @@ func TestEquatorialCoordinateToPixelWithSIPDistortionAtImageCenter(t *testing.T)
CD1_2: 0.0, // No rotation
CD2_1: 0.0, // No rotation
CD2_2: 0.0002777778, // 1/3600 deg/pixel
SIP: sip,
FSIP: fsip,
ISIP: isip,
}

Expand All @@ -394,7 +394,7 @@ func TestEquatorialCoordinateToPixelWithSIPDistortionAtImageCenter(t *testing.T)
/*****************************************************************************************************************/

func TestEquatorialCoordinateToPixelWithSIPDistortion(t *testing.T) {
sip := transform.SIP2DParameters{
fsip := transform.SIP2DForwardParameters{
AOrder: 3,
BOrder: 3,
APower: map[string]float64{
Expand Down Expand Up @@ -452,7 +452,7 @@ func TestEquatorialCoordinateToPixelWithSIPDistortion(t *testing.T) {
CD1_2: 0.0, // No rotation
CD2_1: 0.0, // No rotation
CD2_2: 0.0002777778, // 1/3600 deg/pixel
SIP: sip,
FSIP: fsip,
ISIP: isip,
}

Expand Down

0 comments on commit 276c9c1

Please sign in to comment.