From 5f6c7d54efa2cc2d7e8d3fb45335134b2b9cc7df Mon Sep 17 00:00:00 2001
From: michealroberts <michael@observerly.com>
Date: Tue, 12 Nov 2024 15:19:04 +0000
Subject: [PATCH] feat: add convertPixelToWorldCoordinateSystem() to wcs module
 in @observerly/astrometry

feat: add convertPixelToWorldCoordinateSystem() to wcs module in @observerly/astrometry
---
 src/wcs.ts        | 75 +++++++++++++++++++++++++++++++++++++++++++++++
 tests/wcs.spec.ts | 72 +++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 147 insertions(+)

diff --git a/src/wcs.ts b/src/wcs.ts
index 1fd07de..090464a 100644
--- a/src/wcs.ts
+++ b/src/wcs.ts
@@ -109,5 +109,80 @@ export type WCS = {
   SIP?: SIP2DParameters
 }
 
+/*****************************************************************************************************************/
+
+/**
+ *
+ * convertPixelToWorldCoordinateSystem
+ *
+ * Converts a given image pixel to an equatorial coordinate base on the prescribed WCS
+ *
+ * @param wcs the World Coordinate System (WCS) parameters to use for the conversion
+ * @param pixel the pixel to convert to an equatorial coordinate, e.g., { x: 0, y: 0 }
+ * @returns the equatorial coordinate of a given pixel in the prescribed WCS
+ */
+export const convertPixelToWorldCoordinateSystem = (
+  wcs: WCS,
+  pixel: CartesianCoordinate
+): EquatorialCoordinate => {
+  // Calculate the pixel delta in the X dimension relative to central reference pixel:
+  let deltaX = pixel.x - wcs.crpix.x
+
+  // Calculate the pixel delta in the Y dimension relative to central reference pixel:
+  let deltaY = pixel.y - wcs.crpix.y
+
+  // Initialize SIP correction terms for A and B:
+  let A = 0
+  let B = 0
+
+  // Apply SIP polynomial corrections if SIP parameters are defined
+  if (wcs.SIP) {
+    // Apply A polynomial corrections:
+    for (const term in wcs.SIP.APower) {
+      const coeff = wcs.SIP.APower[term]
+      const indices = parseSIPTerm(term, 'A')
+      if (indices) {
+        const [i, j] = indices
+        A += coeff * deltaX ** i * deltaY ** j
+      }
+    }
+
+    // Apply B polynomial corrections:
+    for (const term in wcs.SIP.BPower) {
+      const coeff = wcs.SIP.BPower[term]
+      const indices = parseSIPTerm(term, 'B')
+      if (indices) {
+        const [i, j] = indices
+        B += coeff * deltaX ** i * deltaY ** j
+      }
+    }
+  }
+
+  // Apply forward SIP transformation to correct for non-linear distortions:
+  deltaX += A
+  deltaY += B
+
+  // Calculate the reference equatorial coordinate for the right ascension:
+  let ra = wcs.cd11 * deltaX + wcs.cd12 * deltaY + wcs.E
+
+  // Correct for large values of RA:
+  ra = ra % 360
+
+  // Correct for negative values of RA:
+  if (ra < 0) {
+    ra += 360
+  }
+
+  // Calculate the reference equatorial coordinate for the declination:
+  let dec = wcs.cd21 * deltaX + wcs.cd22 * deltaY + wcs.F
+
+  // Correct for large values of DEC:
+  dec = dec % 90
+
+  return {
+    ra,
+    dec
+  }
+}
 
 /*****************************************************************************************************************/
diff --git a/tests/wcs.spec.ts b/tests/wcs.spec.ts
index e69de29..5407f30 100644
--- a/tests/wcs.spec.ts
+++ b/tests/wcs.spec.ts
@@ -0,0 +1,72 @@
+/*****************************************************************************************************************/
+
+// @author         Michael Roberts <michael@observerly.com>
+// @package        @observerly/astrometry/wcs
+// @license        Copyright © 2021-2023 observerly
+
+/*****************************************************************************************************************/
+
+import { describe, expect, it } from 'vitest'
+
+import {
+  type SIP2DParameters,
+  WCS,
+  convertPixelToWorldCoordinateSystem
+} from '../src'
+
+/*****************************************************************************************************************/
+
+describe('convertPixelToWorldCoordinateSystem', () => {
+  it('should be defined', () => {
+    expect(convertPixelToWorldCoordinateSystem).toBeDefined()
+  })
+
+  it('should convert pixel coordinates to equatorial coordinates correctly without SIP transformations', () => {
+    const wcs = {
+      crpix: { x: 200, y: 200 },
+      crval: { ra: 0, dec: 0 },
+      cd11: 0.2,
+      cd12: 30,
+      cd21: 0.2,
+      cd22: 0.2,
+      E: 0,
+      F: 0
+    } satisfies WCS
+
+    const pixel = { x: 100, y: 100 }
+
+    const equatorial = convertPixelToWorldCoordinateSystem(wcs, pixel)
+
+    expect(equatorial).toEqual({ ra: 220, dec: -40 })
+  })
+
+  it('should convert pixel coordinates to equatorial coordinates correctly with SIP transformations', () => {
+    const sip = {
+      AOrder: 2,
+      BOrder: 2,
+      APower: { A_1_0: 0, A_0_1: 0, A_1_1: 2.47110429124e-07, A_2_0: 6.71868472834e-07, A_0_2: 2.8223469986e-07 },
+      BPower: { B_1_0: 0, B_0_1: 0, B_1_1: 2.47110429124e-07, B_2_0: 6.71868472834e-07, B_0_2: 2.8223469986e-07 }
+    } satisfies SIP2DParameters
+
+    const wcs = {
+      crpix: { x: 200, y: 200 },
+      crval: { ra: 0, dec: 0 },
+      cd11: 0.2,
+      cd12: 30,
+      cd21: 0.2,
+      cd22: 0.2,
+      E: 0,
+      F: 0,
+      SIP: sip
+    } satisfies WCS
+
+    const pixel = { x: 100, y: 100 }
+
+    const equatorial = convertPixelToWorldCoordinateSystem(wcs, pixel)
+
+    expect(equatorial.ra).toBeCloseTo(220.36, 1)
+    expect(equatorial.dec).toBeCloseTo(-39.99, 1)
+  })
+})
+
+/*****************************************************************************************************************/
\ No newline at end of file