feat: add (ps *PlateSolver) Solve to solver module in @observerly/skysolve

cmd/main.go

@@ -0,0 +1,101 @@
+/*****************************************************************************************************************/
+
+//	@author		Michael Roberts <[email protected]>
+//	@package	@observerly/skysolve
+//	@license	Copyright © 2021-2025 observerly
+
+/*****************************************************************************************************************/
+
+package main
+
+import (
+	"fmt"
+	"os"
+	"time"
+
+	"github.com/observerly/iris/pkg/fits"
+	"github.com/observerly/skysolve/pkg/geometry"
+	"github.com/observerly/skysolve/pkg/solver"
+)
+
+/*****************************************************************************************************************/
+
+func main() {
+	// Attempt to open the file from the given filepath:
+	file, err := os.Open("./samples/noise16.fits")
+
+	if err != nil {
+		fmt.Printf("failed to open file: %v", err)
+		return
+	}
+
+	// Defer closing the file:
+	defer file.Close()
+
+	// Assume an image of 2x2 pixels with 16-bit depth, and no offset:
+	fit := fits.NewFITSImage(2, 0, 0, 65535)
+
+	// Read in our exposure data into the image:
+	err = fit.Read(file)
+
+	if err != nil {
+		fmt.Printf("failed to read file: %v", err)
+		return
+	}
+
+	// Attempt to get the RA header from the FITS file:
+	ra, exists := fit.Header.Floats["RA"]
+	if !exists {
+		fmt.Printf("ra header not found")
+		return
+	}
+
+	// Attempt to get the Dec header from the FITS file:
+	dec, exists := fit.Header.Floats["DEC"]
+	if !exists {
+		fmt.Printf("dec header not found")
+		return
+	}
+
+	// Attempt to create a new PlateSolver:
+	solver, err := solver.NewPlateSolver(solver.GAIA, fit, solver.Params{
+		RA:                  float64(ra.Value),  // The appoximate RA of the center of the image
+		Dec:                 float64(dec.Value), // The appoximate Dec of the center of the image
+		PixelScale:          2.061 / 3600.0,     // 2.061 arcseconds per pixel (0.0005725 degrees)
+		ExtractionThreshold: 80,                 // Extract a minimum of 80 of the brightest stars
+		Radius:              16,                 // 16 pixels radius for the star extraction
+		Sigma:               8,                  // 8 pixels sigma for the Gaussian kernel
+	})
+
+	if err != nil {
+		fmt.Printf("error: %v", err)
+		return
+	}
+
+	// Define the tolerances for the solver, we can adjust these as needed:
+	tolerances := geometry.InvariantFeatureTolerance{
+		LengthRatio: 0.025, // 5% tolerance in side length ratios
+		Angle:       0.5,   // 1 degree tolerance in angles
+	}
+
+	// Record the start time
+	startTime := time.Now()
+	wcs, err := solver.Solve(tolerances)
+
+	if err != nil {
+		fmt.Printf("an error occured while plate solving: %v", err)
+		return
+	}
+
+	// Calculate the elapsed time
+	elapsedTime := time.Since(startTime)
+
+	fmt.Println(elapsedTime)
+
+	// Calculate the reference equatorial coordinate:
+	eq := wcs.PixelToEquatorialCoordinate(578.231147766, 485.620500565)
+
+	fmt.Println(eq)
+}
+
+/*****************************************************************************************************************/

pkg/solver/solver.go

@@ -12,19 +12,24 @@ package solver
 
 import (
 	"errors"
+	"fmt"
 	"math"
 	"sort"
 	"sync"
 
 	"github.com/observerly/iris/pkg/fits"
 	"github.com/observerly/iris/pkg/photometry"
 	stats "github.com/observerly/iris/pkg/statistics"
+	iutils "github.com/observerly/iris/pkg/utils"
 
 	"github.com/observerly/skysolve/pkg/astrometry"
 	"github.com/observerly/skysolve/pkg/catalog"
 	"github.com/observerly/skysolve/pkg/geometry"
+	"github.com/observerly/skysolve/pkg/matrix"
 	"github.com/observerly/skysolve/pkg/projection"
+	"github.com/observerly/skysolve/pkg/transform"
 	"github.com/observerly/skysolve/pkg/utils"
+	"github.com/observerly/skysolve/pkg/wcs"
 )
 
 /*****************************************************************************************************************/
@@ -417,3 +422,109 @@ func (ps *PlateSolver) FindSourceMatches(tolerance geometry.InvariantFeatureTole
 }
 
 /*****************************************************************************************************************/
+
+func (ps *PlateSolver) Solve(tolerance geometry.InvariantFeatureTolerance) (*wcs.WCS, error) {
+	matches, err := ps.FindSourceMatches(tolerance)
+
+	if err != nil {
+		return nil, err
+	}
+
+	if len(matches) < 3 {
+		return nil, errors.New("insufficient matches to perform plate solving")
+	}
+
+	n := 2 * len(matches)
+	A := make([][]float64, n)
+	B := make([]float64, n)
+
+	// Calculate the affine transformation matrix from the matches:
+	for i, match := range matches {
+		x := float64(match.Star.X)
+		y := float64(match.Star.Y)
+
+		ra := match.Source.RA
+		dec := match.Source.Dec
+
+		// Create the matrix A and vector B for the least squares method:
+		A[2*i] = []float64{x, y, 1, 0, 0, 0}
+		B[2*i] = ra
+		A[2*i+1] = []float64{0, 0, 0, x, y, 1}
+		B[2*i+1] = dec
+	}
+
+	// Convert A to a matrix:
+	a, err := matrix.NewFromSlice(iutils.Flatten2DFloat64Array(A), n, 6)
+	if err != nil {
+		return nil, err
+	}
+
+	// Convert B to a matrix:
+	b, err := matrix.NewFromSlice(B, n, 1)
+	if err != nil {
+		return nil, err
+	}
+
+	// Compute A^T
+	aT, err := a.Transpose()
+
+	if err != nil {
+		return nil, fmt.Errorf("failed to compute A^T: %v", err)
+	}
+
+	// Compute A^T * A
+	aTa, err := aT.Multiply(a)
+	if err != nil {
+		return nil, fmt.Errorf("failed to compute A^T * A: %v", err)
+	}
+
+	// Compute A^T * B
+	aTb, err := aT.Multiply(b)
+	if err != nil {
+		return nil, fmt.Errorf("failed to compute A^T * B: %v", err)
+	}
+
+	// Compute the inverse of A^T * A
+	aTaInv, err := aTa.Invert()
+	if err != nil {
+		return nil, fmt.Errorf("failed to invert A^T * A: %v", err)
+	}
+
+	// Compute the affine transformation matrix parameters using the least squares method:
+	params := make([]float64, 6)
+	if aTaInv == nil || aTb == nil {
+		return nil, errors.New("failed to compute affine transformation matrix parameters")
+	}
+
+	// Calculate the affine transformation matrix parameters:
+	for i := 0; i < 6; i++ {
+		for j := 0; j < 6; j++ {
+			params[i] += aTaInv.Value[i*6+j] * aTb.Value[j]
+		}
+	}
+
+	// Calculate the x-coordinate of the center of the image:
+	x := float64(ps.Width) / 2
+
+	// Calculate the y-coordinate of the center of the image:
+	y := float64(ps.Height) / 2
+
+	// Now that we have the affine parameters, we can calculate the actual RA and dec coordinate
+	// for the center of the image:
+	t := wcs.NewWorldCoordinateSystem(
+		x,
+		y,
+		transform.Affine2DParameters{
+			A: params[0],
+			B: params[1],
+			C: params[3],
+			D: params[4],
+			E: params[2],
+			F: params[5],
+		},
+	)
+
+	return &t, nil
+}
+
+/*****************************************************************************************************************/

pkg/solver/solver_test.go

@@ -0,0 +1,112 @@
+/*****************************************************************************************************************/
+
+//	@author		Michael Roberts <[email protected]>
+//	@package	@observerly/skysolve
+//	@license	Copyright © 2021-2025 observerly
+
+/*****************************************************************************************************************/
+
+package solver
+
+/*****************************************************************************************************************/
+
+import (
+	"math"
+	"os"
+	"testing"
+	"time"
+
+	"github.com/observerly/iris/pkg/fits"
+	"github.com/observerly/skysolve/pkg/geometry"
+)
+
+/*****************************************************************************************************************/
+
+func TestSolverOnMatches(t *testing.T) {
+	// Attempt to open the file from the given filepath:
+	file, err := os.Open("../../samples/noise16.fits")
+
+	if err != nil {
+		t.Errorf("NewFITSImageFromReader() os.Open(): %v", err)
+	}
+
+	// Defer closing the file:
+	defer file.Close()
+
+	// Assume an image of 2x2 pixels with 16-bit depth, and no offset:
+	fit := fits.NewFITSImage(2, 0, 0, 65535)
+
+	// Read in our exposure data into the image:
+	err = fit.Read(file)
+
+	if err != nil {
+		t.Errorf("Read() error: %v", err)
+	}
+
+	// Attempt to get the RA header from the FITS file:
+	ra, exists := fit.Header.Floats["RA"]
+	if !exists {
+		t.Errorf("ra header not found")
+	}
+
+	// Attempt to get the Dec header from the FITS file:
+	dec, exists := fit.Header.Floats["DEC"]
+	if !exists {
+		t.Errorf("dec header not found")
+	}
+
+	// Attempt to create a new PlateSolver:
+	solver, err := NewPlateSolver(GAIA, fit, Params{
+		RA:                  float64(ra.Value),  // The appoximate RA of the center of the image
+		Dec:                 float64(dec.Value), // The appoximate Dec of the center of the image
+		PixelScale:          2.061 / 3600.0,     // 2.061 arcseconds per pixel (0.0005725 degrees)
+		ExtractionThreshold: 80,                 // Extract a minimum of 80 of the brightest stars
+		Radius:              16,                 // 16 pixels radius for the star extraction
+		Sigma:               8,                  // 8 pixels sigma for the Gaussian kernel
+	})
+
+	if err != nil {
+		t.Errorf("error: %v", err)
+	}
+
+	// Define the tolerances for the solver, we can adjust these as needed:
+	tolerances := geometry.InvariantFeatureTolerance{
+		LengthRatio: 0.025, // 5% tolerance in side length ratios
+		Angle:       0.5,   // 1 degree tolerance in angles
+	}
+
+	// Record the start time
+	startTime := time.Now()
+	wcs, err := solver.Solve(tolerances)
+
+	if err != nil {
+		t.Errorf("error: %v", err)
+		return
+	}
+
+	// Calculate the elapsed time
+	elapsedTime := time.Since(startTime)
+
+	// Calculate the reference equatorial coordinate:
+	eq := wcs.PixelToEquatorialCoordinate(578.231147766, 485.620500565)
+
+	// We cross-reference here with calibration data from the astrometry.net API:
+	if math.Abs(eq.RA-98.64649870834154) > 0.0001 {
+		t.Errorf("RA not set correctly")
+	}
+
+	// We cross-reference here with calibration data from the astrometry.net API:
+	if math.Abs(eq.Dec-2.537618674632346) > 0.0001 {
+		t.Errorf("Dec not set correctly")
+	}
+
+	if elapsedTime.Seconds() > 0.4 {
+		t.Errorf("plate solver took too long to execute")
+	}
+
+	t.Logf("solver.Solve(tolerances) completed in %v", elapsedTime)
+
+	t.Logf("RA: %v, Dec: %v", eq.RA, eq.Dec)
+}
+
+/*****************************************************************************************************************/