Skip to content

Commit

Permalink
🐛 Fix rectangle rotation bug using a new conversion algorithm
Browse files Browse the repository at this point in the history
  • Loading branch information
@Xen0Xys
Xen0Xys committed May 17, 2024
1 parent 88db505 commit cfa5513
Showing 1 changed file with 41 additions and 137 deletions.
178 changes: 41 additions & 137 deletions src/ipyaladin/region_converter.py
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,7 @@
import math

import numpy as np

from astropy.coordinates import SkyCoord
from astropy.coordinates import SkyCoord, CartesianRepresentation, Angle
from astropy.coordinates.matrix_utilities import rotation_matrix

try:
from regions import (
Expand All @@ -22,108 +21,10 @@
LineSkyRegion = None
from typing import Union

TWICE_PI = 2 * math.pi


class RefToLocalRotMatrix:
"""A rotation matrix from the reference frame to the local frame.
The reference frame is the ICRS frame. The local frame is
defined by the center of the local reference frame.
Attributes
----------
rotation_matrix : np.ndarray[float]
The rotation matrix.
"""

def __init__(
self,
rotation_matrix: np.ndarray[float],
) -> None:
self.rotation_matrix = rotation_matrix

@classmethod
def from_center(cls: any, lon: float, lat: float) -> "RefToLocalRotMatrix":
"""Create a rotation matrix from the center of the local reference frame.
Parameters
----------
lon : float
The longitude of the center.
lat : float
The latitude of the center.
Returns
-------
RefToLocalRotMatrix
The rotation matrix.
"""
ca, sa = math.cos(lon), math.sin(lon)
cd, sd = math.cos(lat), math.sin(lat)
return cls(
np.array([[ca * cd, sa * cd, sd], [-sa, ca, 0.0], [-ca * sd, -sa * sd, cd]])
)

def to_global_xyz(self, x: float, y: float, z: float) -> tuple[float, float, float]:
"""Convert local cartesian coordinates to global cartesian coordinates.
Parameters
----------
x : float
The x coordinate.
y : float
The y coordinate.
z : float
The z coordinate.
Returns
-------
tuple[float, float, float]
The x, y, and z coordinates.
"""
rotated_matrix = np.sum(
self.rotation_matrix * np.array([[x], [y], [z]]), axis=0
)
return rotated_matrix[0], rotated_matrix[1], rotated_matrix[2]

def to_global_coo(self, x: float, y: float, z: float) -> tuple[float, float]:
"""Convert local cartesian coordinates to global spherical coordinates.
Parameters
----------
x : float
The x coordinate.
y : float
The y coordinate.
z : float
The z coordinate.
Returns
-------
tuple[float, float]
The longitude and latitude.
"""
x, y, z = self.to_global_xyz(x, y, z)
# Convert cartesian coordinates to spherical coordinates.
r2 = x * x + y * y
lat = math.atan2(z, math.sqrt(r2))
lon = math.atan2(y, x)
if lon < 0:
lon += TWICE_PI
return lon, lat


def rectangle_to_polygon_region(region: RectangleSkyRegion) -> PolygonSkyRegion:
"""Convert a RectangleSkyRegion to a PolygonSkyRegion.
Converted code from cdshealpix:
https://github.com/cds-astro/cds-healpix-rust/blob/master/src/nested/mod.rs#L3687-3742
Parameters
----------
region : RectangleSkyRegion
Expand All @@ -135,47 +36,50 @@ def rectangle_to_polygon_region(region: RectangleSkyRegion) -> PolygonSkyRegion:
The converted region.
"""
longitude, latitude = region.center.icrs.ra.rad, region.center.icrs.dec.rad
width = region.width.rad
height = region.height.rad
position_angle = region.angle.rad

if not 0.0 <= longitude < TWICE_PI:
raise ValueError(f"Expected: lon in [0, 2pi[. Actual: {longitude}")
if not -math.pi / 2 <= latitude <= math.pi / 2:
raise ValueError(f"Expected: lat in [-pi/2, pi/2]. Actual: {latitude}")
if not 0.0 < width <= math.pi / 2:
raise ValueError(f"Expected: a in ]0, pi/2]. Actual: {width}")
if not 0.0 < height <= width:
raise ValueError(f"Expected: b in ]0, a]. Actual: {height}")
if not 0.0 <= position_angle < math.pi:
raise ValueError(f"Expected: pa in [0, pi[. Actual: {position_angle}")

frame_rotation = RefToLocalRotMatrix.from_center(longitude, latitude)

sin_lon, cos_lon = math.sin(width), math.cos(width)
latitude = math.atan(cos_lon * math.tan(height))
sin_lat, cos_lat = math.sin(latitude), math.cos(latitude)
sin_pa, cos_pa = math.sin(position_angle), math.cos(position_angle)

x1, y1, z1 = cos_lon * cos_lat, sin_lon * cos_lat, sin_lat
bl_x, bl_y, bl_z = (
region.center.spherical_offsets_by(-region.width / 2, -region.height / 2)
.represent_as("cartesian")
.xyz.value
)
br_x, br_y, br_z = (
region.center.spherical_offsets_by(-region.width / 2, region.height / 2)
.represent_as("cartesian")
.xyz.value
)
tr_x, tr_y, tr_z = (
region.center.spherical_offsets_by(region.width / 2, region.height / 2)
.represent_as("cartesian")
.xyz.value
)
rl_x, tly, tl_z = (
region.center.spherical_offsets_by(region.width / 2, -region.height / 2)
.represent_as("cartesian")
.xyz.value
)
c_x, c_y, c_z = region.center.represent_as("cartesian").xyz.value

pa_matrix = np.array([[y1 * sin_pa, z1 * sin_pa], [y1 * cos_pa, z1 * sin_pa]])
rot_mat = rotation_matrix(Angle(-region.angle.deg, unit="deg"), (c_x, c_y, c_z))

multipliers = [
[[1, -1], [1, 1]],
[[1, 1], [1, -1]],
[[-1, 1], [-1, -1]],
[[-1, -1], [-1, 1]],
corners = np.array(
[
[bl_x, bl_y, bl_z],
[br_x, br_y, br_z],
[tr_x, tr_y, tr_z],
[rl_x, tly, tl_z],
]
)
rotated_corners = np.dot(corners, rot_mat)
vertices = [
SkyCoord(CartesianRepresentation(*corner), frame="icrs")
for corner in rotated_corners
]
vertices = []

for multiplier in multipliers:
y2, z2 = np.sum(np.array(multiplier) * pa_matrix, axis=1)
vertices.append(frame_rotation.to_global_coo(x1, y2, z2))

return PolygonSkyRegion(
vertices=SkyCoord(vertices, unit="rad", frame="icrs"),
vertices=SkyCoord(
[coord.icrs.ra.deg for coord in vertices],
[coord.icrs.dec.deg for coord in vertices],
unit="deg",
),
visual=region.visual,
meta=region.meta,
)
Expand Down

0 comments on commit cfa5513

Please sign in to comment.