From 86a4ae3052cfe48273da3e17330167577d2b8f79 Mon Sep 17 00:00:00 2001 From: Vladislav Tsendrovskii Date: Sun, 13 Aug 2023 18:48:46 +0500 Subject: [PATCH 1/7] Structurizing code - place projection math code to separate files Adding c movement methods fix usage test Optimization --- .vscode/settings.json | 5 +- setup.py | 29 +- .../library/image_process/distorsion.py | 1 - .../library/movement/basic_movement.py | 5 +- src/vstarstack/library/movement/find_shift.py | 4 +- src/vstarstack/library/movement/flat.py | 32 +- src/vstarstack/library/movement/move_image.py | 61 ++-- .../library/movement/movements/lib/flat.c | 17 + .../library/movement/movements/lib/flat.h | 17 + .../library/movement/movements/lib/sphere.c | 125 ++++++++ .../library/movement/movements/lib/sphere.h | 47 +++ .../library/movement/movements/module.c | 256 +++++++++++++++ src/vstarstack/library/movement/sphere.py | 41 +-- src/vstarstack/library/projection/__init__.py | 0 .../library/projection/equirectangular.c | 117 ------- .../library/projection/orthographic.c | 144 --------- .../library/projection/perspective.c | 128 -------- .../projections/lib/equirectangular.c | 47 +++ .../projections/lib/equirectangular.h | 33 ++ .../projection/projections/lib/orthographic.c | 82 +++++ .../projection/projections/lib/orthographic.h | 40 +++ .../projection/projections/lib/perspective.c | 64 ++++ .../projection/projections/lib/perspective.h | 40 +++ .../projection/projections/lib/projection.h | 23 ++ .../projections/lib/projection_common.h | 25 ++ .../library/projection/projections/module.c | 291 ++++++++++++++++++ .../projections/projection_module.h | 39 +++ src/vstarstack/library/projection/tools.py | 24 +- tests/test_equirectangular.py | 26 +- tests/test_movement_flat.py | 231 +++++++------- tests/test_movement_sphere.py | 149 ++++----- tests/test_projection_orthographic.py | 28 +- tests/test_projection_perspective.py | 27 +- tests/test_stars_pipeline.py | 5 +- tests/test_usage.py | 2 +- 35 files changed, 1484 insertions(+), 721 deletions(-) create mode 100644 src/vstarstack/library/movement/movements/lib/flat.c create mode 100644 src/vstarstack/library/movement/movements/lib/flat.h create mode 100644 src/vstarstack/library/movement/movements/lib/sphere.c create mode 100644 src/vstarstack/library/movement/movements/lib/sphere.h create mode 100644 src/vstarstack/library/movement/movements/module.c delete mode 100644 src/vstarstack/library/projection/__init__.py delete mode 100644 src/vstarstack/library/projection/equirectangular.c delete mode 100644 src/vstarstack/library/projection/orthographic.c delete mode 100644 src/vstarstack/library/projection/perspective.c create mode 100644 src/vstarstack/library/projection/projections/lib/equirectangular.c create mode 100644 src/vstarstack/library/projection/projections/lib/equirectangular.h create mode 100644 src/vstarstack/library/projection/projections/lib/orthographic.c create mode 100644 src/vstarstack/library/projection/projections/lib/orthographic.h create mode 100644 src/vstarstack/library/projection/projections/lib/perspective.c create mode 100644 src/vstarstack/library/projection/projections/lib/perspective.h create mode 100644 src/vstarstack/library/projection/projections/lib/projection.h create mode 100644 src/vstarstack/library/projection/projections/lib/projection_common.h create mode 100644 src/vstarstack/library/projection/projections/module.c create mode 100644 src/vstarstack/library/projection/projections/projection_module.h diff --git a/.vscode/settings.json b/.vscode/settings.json index 9642cacb..cb4692e5 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -1,7 +1,10 @@ { "files.associations": { "ndarraytypes.h": "c", - "random": "c" + "random": "c", + "math.h": "c", + "unistd.h": "c", + "sphere.h": "c" }, "editor.tabSize": 4, "editor.insertSpaces": true, diff --git a/setup.py b/setup.py index 2891dbda..9eb17bdf 100644 --- a/setup.py +++ b/setup.py @@ -16,14 +16,24 @@ import numpy as np from setuptools import setup, Extension -perspective = Extension(name="vstarstack.library.projection.perspective", - sources=["src/vstarstack/library/projection/perspective.c"]) +projection = Extension( name="vstarstack.library.projection.projections", + sources=[ + "src/vstarstack/library/projection/projections/module.c", + "src/vstarstack/library/projection/projections/lib/perspective.c", + "src/vstarstack/library/projection/projections/lib/orthographic.c", + "src/vstarstack/library/projection/projections/lib/equirectangular.c", + ]) -equirectangular = Extension(name="vstarstack.library.projection.equirectangular", - sources=["src/vstarstack/library/projection/equirectangular.c"]) - -orthographic = Extension(name="vstarstack.library.projection.orthographic", - sources=["src/vstarstack/library/projection/orthographic.c"]) +movements = Extension( name="vstarstack.library.movement.movements", + sources=[ + "src/vstarstack/library/movement/movements/module.c", + "src/vstarstack/library/movement/movements/lib/sphere.c", + "src/vstarstack/library/movement/movements/lib/flat.c", + ], + include_dirs=[ + "src/vstarstack/library/projection/projections", + np.get_include(), + ]) image_wave = Extension(name="vstarstack.library.fine_shift.image_wave", sources=["src/vstarstack/library/fine_shift/image_wave.c", @@ -51,9 +61,8 @@ scripts = ['bin/vstarstack'], package_dir = {'': 'src'}, packages=packages, - ext_modules = [perspective, - equirectangular, - orthographic, + ext_modules = [projection, + movements, image_wave], install_requires = [ 'numpy', diff --git a/src/vstarstack/library/image_process/distorsion.py b/src/vstarstack/library/image_process/distorsion.py index d8041744..31a1506b 100644 --- a/src/vstarstack/library/image_process/distorsion.py +++ b/src/vstarstack/library/image_process/distorsion.py @@ -17,7 +17,6 @@ import vstarstack.library.data import vstarstack.library.common -import vstarstack.library.projection.perspective import vstarstack.library.projection.tools class Distorsion: diff --git a/src/vstarstack/library/movement/basic_movement.py b/src/vstarstack/library/movement/basic_movement.py index 83252a1d..ce9d422c 100644 --- a/src/vstarstack/library/movement/basic_movement.py +++ b/src/vstarstack/library/movement/basic_movement.py @@ -13,6 +13,7 @@ # along with this program. If not, see . # +import numpy as np from abc import ABC, abstractmethod class MovementException(Exception): @@ -26,12 +27,12 @@ class Movement(ABC): """Interface of movements""" @abstractmethod - def apply(self, positions : list, proj) -> list: + def apply(self, positions : np.ndarray) -> np.ndarray: """Apply movement to positions""" return [] @abstractmethod - def reverse(self, positions : list, proj) -> list: + def reverse(self, positions : np.ndarray) -> np.ndarray: """Apply reverse movement to positions""" @abstractmethod diff --git a/src/vstarstack/library/movement/find_shift.py b/src/vstarstack/library/movement/find_shift.py index 34e7185a..27cacd09 100644 --- a/src/vstarstack/library/movement/find_shift.py +++ b/src/vstarstack/library/movement/find_shift.py @@ -50,8 +50,8 @@ def build_movements(movement : Movement, clusters : list): if name2 not in cluster: continue - star_to = (cluster[name1]["lat"], cluster[name1]["lon"]) - star_from = (cluster[name2]["lat"], cluster[name2]["lon"]) + star_to = (cluster[name1]["lon"], cluster[name1]["lat"]) + star_from = (cluster[name2]["lon"], cluster[name2]["lat"]) stars.append((star_from, star_to)) if len(stars) >= 2: diff --git a/src/vstarstack/library/movement/flat.py b/src/vstarstack/library/movement/flat.py index dfb3066c..87d69ddf 100644 --- a/src/vstarstack/library/movement/flat.py +++ b/src/vstarstack/library/movement/flat.py @@ -25,33 +25,33 @@ class Movement(vstarstack.library.movement.basic_movement.Movement): """Flat movements of plane""" - def apply(self, positions : list, proj) -> List: + def apply(self, positions : np.ndarray) -> List: """Apply movement""" npositions = [] - for y, x in positions: + for x, y in positions: nx = x * math.cos(self.a) - y * math.sin(self.a) + self.dx ny = y * math.cos(self.a) + x * math.sin(self.a) + self.dy - npositions.append((ny, nx)) + npositions.append((nx, ny)) - return npositions + return np.array(npositions) - def reverse(self, positions : list, proj) -> List: + def reverse(self, positions : np.ndarray) -> List: """Apply reverse movement""" npositions = [] - for y, x in positions: + for x, y in positions: nx = (x-self.dx) * math.cos(self.a) + \ (y-self.dy) * math.sin(self.a) ny = (y-self.dy) * math.cos(self.a) - \ (x-self.dx) * math.sin(self.a) - npositions.append((ny, nx)) + npositions.append((nx, ny)) - return npositions + return np.array(npositions) def magnitude(self) -> float: """Magnitude of movement""" return self.dx**2 + self.dy**2 + self.a**2 - def __init__(self, angle, dy, dx): + def __init__(self, angle, dx, dy): self.dx = dx self.dy = dy self.a = angle @@ -77,7 +77,7 @@ def build(point1_from, point2_from, point1_to, point2_to, debug=False): dir_to = norm((point2_to[0] - point1_to[0], point2_to[1] - point1_to[1])) cosa = dir_from[0]*dir_to[0] + dir_from[1]*dir_to[1] - sina = dir_from[1]*dir_to[0] - dir_from[0]*dir_to[1] + sina = dir_from[0]*dir_to[1] - dir_from[1]*dir_to[0] cosa = np.clip(cosa, -1, 1) sina = np.clip(sina, -1, 1) @@ -87,10 +87,10 @@ def build(point1_from, point2_from, point1_to, point2_to, debug=False): angle = math.pi - angle transformation = Movement(angle, 0, 0) - point1_rotated = transformation.apply([point1_from], None)[0] - dy = point1_to[0] - point1_rotated[0] - dx = point1_to[1] - point1_rotated[1] - return Movement(angle, dy, dx) + point1_rotated = transformation.apply([point1_from])[0] + dx = point1_to[0] - point1_rotated[0] + dy = point1_to[1] - point1_rotated[1] + return Movement(angle, dx, dy) @staticmethod def average(transformations : list): @@ -121,9 +121,9 @@ def __mul__(self, other): angle = angle1 + angle2 dx = dx2 * math.cos(angle1) - dy2 * math.sin(angle1) + dx1 dy = dx2 * math.sin(angle1) + dy2 * math.cos(angle1) + dy1 - return Movement(angle, dy, dx) + return Movement(angle, dx, dy) def inverse(self): idx = -self.dx * math.cos(self.a) - self.dy * math.sin(self.a) idy = -self.dy * math.cos(self.a) + self.dx * math.sin(self.a) - return Movement(-self.a, idy, idx) + return Movement(-self.a, idx, idy) diff --git a/src/vstarstack/library/movement/move_image.py b/src/vstarstack/library/movement/move_image.py index 6ecde706..2d40d3d8 100644 --- a/src/vstarstack/library/movement/move_image.py +++ b/src/vstarstack/library/movement/move_image.py @@ -14,6 +14,8 @@ # import numpy as np +import scipy.ndimage + import vstarstack.library.common import vstarstack.library.data from vstarstack.library.data import DataFrame @@ -21,22 +23,22 @@ import vstarstack.library.projection.tools import vstarstack.library.movement.basic_movement as basic_movement -def _generate_points(height, width, len0): +from vstarstack.library.common import getpixel + +def _generate_points(height, width): """Generate grid of pixel coordinates""" - points = [] + points = np.zeros((height*width, 2), dtype='int') for y in range(height): for x in range(width): - points.append((y, x)) - if len(points) >= len0: - yield points - points = [] - if len(points) > 0: - yield points - -def move_image(image : np.ndarray, - transformation : basic_movement.Movement, - proj, - image_weight : float, + points[y * width + x, 0] = x + points[y * width + x, 1] = y + return points + + +def move_image(image: np.ndarray, + transformation: basic_movement.Movement, + input_proj, output_proj, + image_weight: float, image_weight_layer=None): """Apply movement to image""" shape = image.shape @@ -49,23 +51,26 @@ def move_image(image : np.ndarray, if image_weight_layer is None: image_weight_layer = np.ones(shape)*image_weight - for positions in _generate_points(h, w, w*4): - original_positions = transformation.reverse(positions, proj) - for position, original_position in zip(positions, original_positions): - y, x = position - orig_y, orig_x = original_position - _, pixel = vstarstack.library.common.getpixel(image, orig_y, orig_x, False) - _, pixel_weight = vstarstack.library.common.getpixel( - image_weight_layer, orig_y, orig_x, False) - - shifted[y][x] = pixel - shifted_weight_layer[y][x] = pixel_weight + positions = _generate_points(h, w) + original_positions = transformation.reverse(positions.astype('double'), input_proj, output_proj) + num = positions.shape[0] + transform_array = np.zeros([h, w, 2], dtype='double') + for index in range(num): + position = positions[index] + original_position = original_positions[index] + x, y = position[0], position[1] + orig_x, orig_y = original_position[0], original_position[1] + transform_array[y, x, 0] = orig_y + transform_array[y, x, 1] = orig_x + crdtf = lambda pos : tuple(transform_array[pos[0], pos[1], :]) + shifted = scipy.ndimage.geometric_transform(image, crdtf, order=3) + shifted_weight_layer = scipy.ndimage.geometric_transform(image_weight_layer, crdtf, order=3) return shifted, shifted_weight_layer -def move_dataframe(dataframe : DataFrame, - transformation : basic_movement.Movement, - proj = None): +def move_dataframe(dataframe: DataFrame, + transformation: basic_movement.Movement, + proj=None): """Apply movement to dataframe""" if proj is None: proj = vstarstack.library.projection.tools.get_projection(dataframe) @@ -86,7 +91,7 @@ def move_dataframe(dataframe : DataFrame, else: weight = np.ones(image.shape)*1 - shifted, shifted_weight = move_image(image, transformation, proj, weight) + shifted, shifted_weight = move_image(image, transformation, proj, proj, weight) dataframe.add_channel(shifted, channel, **opts) dataframe.add_channel(shifted_weight, weight_channel, weight=True) dataframe.add_channel_link(channel, weight_channel, "weight") diff --git a/src/vstarstack/library/movement/movements/lib/flat.c b/src/vstarstack/library/movement/movements/lib/flat.c new file mode 100644 index 00000000..f7abb7e8 --- /dev/null +++ b/src/vstarstack/library/movement/movements/lib/flat.c @@ -0,0 +1,17 @@ +/* + * Copyright (c) 2023 Vladislav Tsendrovskii + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, version 3 of the License. + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU General Public License for more details. + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include "flat.h" + diff --git a/src/vstarstack/library/movement/movements/lib/flat.h b/src/vstarstack/library/movement/movements/lib/flat.h new file mode 100644 index 00000000..8bc191ea --- /dev/null +++ b/src/vstarstack/library/movement/movements/lib/flat.h @@ -0,0 +1,17 @@ +/* + * Copyright (c) 2023 Vladislav Tsendrovskii + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, version 3 of the License. + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU General Public License for more details. + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#pragma once + +#include \ No newline at end of file diff --git a/src/vstarstack/library/movement/movements/lib/sphere.c b/src/vstarstack/library/movement/movements/lib/sphere.c new file mode 100644 index 00000000..f9b63336 --- /dev/null +++ b/src/vstarstack/library/movement/movements/lib/sphere.c @@ -0,0 +1,125 @@ +/* + * Copyright (c) 2023 Vladislav Tsendrovskii + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, version 3 of the License. + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU General Public License for more details. + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include +#include "sphere.h" + +static struct quat quat_inv(const struct quat q) +{ + double len2 = q.w*q.w+q.x*q.x+q.y*q.y+q.z*q.z; + struct quat rev = { + .w = q.w/len2, + .x = -q.x/len2, + .y = -q.y/len2, + .z = -q.z/len2, + }; + return rev; +} + +static struct quat quat_mul(const struct quat a, const struct quat b) +{ + struct quat q = { + .w = a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z, + .x = a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y, + .y = a.w * b.y - a.x * b.z + a.y * b.w + a.z * b.x, + .z = a.w * b.z + a.x * b.y - a.y * b.x + a.z * b.w, + }; + return q; +} + +static struct quat quat_from_vec(double x, double y, double z) +{ + struct quat q = { + .w = 0, + .x = x, + .y = y, + .z = z, + }; + return q; +} + +void sphere_movement_init(struct SphereMovement *mov, struct quat q) +{ + mov->forward = q; + mov->reverse = quat_inv(q); +} + +void sphere_movement_apply_forward(struct SphereMovement *mov, + const double *posi, double *poso, size_t num, + const struct ProjectionDef *in_proj, + const struct ProjectionDef *out_proj) +{ + unsigned i; + for (i = 0; i < num; i++) + { + double lon, lat; + double x, y, z; + double xi = posi[i*2]; + double yi = posi[i*2+1]; + in_proj->forward(in_proj->projection, yi, xi, &lat, &lon); + x = cos(lon)*cos(lat); + y = sin(lon)*cos(lat); + z = sin(lat); + struct quat v = quat_from_vec(x, y, z); + struct quat vf = quat_mul(mov->forward, quat_mul(v, mov->reverse)); + x = vf.x; + y = vf.y; + z = vf.z; + if (z > 1) + z = 1; + if (z < -1) + z = -1; + lon = atan2(y, x); + lat = asin(z); + double yo, xo; + out_proj->reverse(out_proj->projection, lat, lon, &yo, &xo); + poso[2*i] = xo; + poso[2*i+1] = yo; + } +} + +void sphere_movement_apply_reverse(struct SphereMovement *mov, + const double *posi, double *poso, size_t num, + const struct ProjectionDef *in_proj, + const struct ProjectionDef *out_proj) +{ + unsigned i; + for (i = 0; i < num; i++) + { + double lon, lat; + double x, y, z; + double xi = posi[i*2]; + double yi = posi[i*2+1]; + out_proj->forward(out_proj->projection, yi, xi, &lat, &lon); + x = cos(lon)*cos(lat); + y = sin(lon)*cos(lat); + z = sin(lat); + struct quat v = quat_from_vec(x, y, z); + struct quat vf = quat_mul(mov->reverse, quat_mul(v, mov->forward)); + x = vf.x; + y = vf.y; + z = vf.z; + if (z > 1) + z = 1; + if (z < -1) + z = -1; + lon = atan2(y, x); + lat = asin(z); + double yo, xo; + in_proj->reverse(in_proj->projection, lat, lon, &yo, &xo); + poso[2*i] = xo; + poso[2*i+1] = yo; + } +} diff --git a/src/vstarstack/library/movement/movements/lib/sphere.h b/src/vstarstack/library/movement/movements/lib/sphere.h new file mode 100644 index 00000000..f8d6d51d --- /dev/null +++ b/src/vstarstack/library/movement/movements/lib/sphere.h @@ -0,0 +1,47 @@ +/* + * Copyright (c) 2023 Vladislav Tsendrovskii + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, version 3 of the License. + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU General Public License for more details. + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#pragma once + +#include +#include "lib/projection.h" + +struct quat { + double w, x, y, z; +}; + +struct SphereMovement +{ + struct quat forward; // forward quaternion + struct quat reverse; // reverse quaternion +}; + +struct ProjectionDef +{ + void *projection; + forward_project_f forward; + reverse_project_f reverse; +}; + +void sphere_movement_init(struct SphereMovement *mov, struct quat q); + +void sphere_movement_apply_forward(struct SphereMovement *mov, + const double *posi, double *poso, size_t num, + const struct ProjectionDef *in_proj, + const struct ProjectionDef *out_proj); + +void sphere_movement_apply_reverse(struct SphereMovement *mov, + const double *posi, double *poso, size_t num, + const struct ProjectionDef *in_proj, + const struct ProjectionDef *out_proj); diff --git a/src/vstarstack/library/movement/movements/module.c b/src/vstarstack/library/movement/movements/module.c new file mode 100644 index 00000000..e061395a --- /dev/null +++ b/src/vstarstack/library/movement/movements/module.c @@ -0,0 +1,256 @@ +/* + * Copyright (c) 2023 Vladislav Tsendrovskii + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, version 3 of the License. + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU General Public License for more details. + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + + +#define PY_SSIZE_T_CLEAN +#include +#include +#include +#include +#include + +#include "projection_module.h" + +#include "lib/sphere.h" +#include "lib/flat.h" + +#define BASENAME "vstarstack.library.movement.movements" + +struct SphereMovementObject +{ + PyObject_HEAD + struct SphereMovement mov; +}; + +static bool generic_forward_project(void *self, + double y, double x, + double *lat, double *lon) +{ + PyObject *projection = (PyObject *)self; + PyObject *res = PyObject_CallMethod(projection, "project", "(dd)", x, y); + PyObject *_lat = PyTuple_GetItem(res, 1); + PyObject *_lon = PyTuple_GetItem(res, 0); + if (_lat == NULL || _lon == NULL) + return false; + *lat = PyFloat_AsDouble(_lat); + *lon = PyFloat_AsDouble(_lon); + //printf("xy %lf:%lf -> lonlat %lf:%lf\n", x, y, *lon, *lat); + return true; +} + +static bool generic_reverse_project(void *self, + double lat, double lon, + double *y, double *x) +{ + PyObject *projection = (PyObject *)self; + PyObject *res = PyObject_CallMethod(projection, "reverse", "(dd)", lon, lat); + PyObject *_y = PyTuple_GetItem(res, 1); + PyObject *_x = PyTuple_GetItem(res, 0); + if (_y == NULL || _x == NULL) + return false; + *y = PyFloat_AsDouble(_y); + *x = PyFloat_AsDouble(_x); + //printf("lonlat %lf:%lf -> xy %lf:%lf\n", lon, lat, *x, *y); + return true; +} + +static int SphereMovements_init(PyObject *_self, PyObject *args, PyObject *kwds) +{ + double w, x, y, z; + struct SphereMovementObject *self = (struct SphereMovementObject *)_self; + static char *kwlist[] = {"w", "x", "y", "z", NULL}; + if (!PyArg_ParseTupleAndKeywords(args, kwds, "dddd", kwlist, + &w, &x, &y, &z)) + { + return -1; + } + + struct quat q = { + .w = w, + .x = x, + .y = y, + .z = z, + }; + + sphere_movement_init(&self->mov, q); + return 0; +} + +bool is_perspective(PyObject *obj) +{ + return false; +} + +typedef void (*action_f)(struct SphereMovement *mov, + const double *posi, double *poso, size_t num, + const struct ProjectionDef *in_proj, + const struct ProjectionDef *out_proj); + +static PyObject *apply_action(PyObject *_self, + PyObject *args, + PyObject *kwds, + action_f fun) +{ + PyObject *input_proj, *output_proj; + PyArrayObject *points; + struct SphereMovementObject *self = (struct SphereMovementObject *)_self; + + static char *kwlist[] = {"points", "input_projection", "output_projection", NULL}; + if (!PyArg_ParseTupleAndKeywords(args, kwds, "OOO", kwlist, + &points, &input_proj, &output_proj)) + { + PyErr_SetString(PyExc_ValueError, "invalid function arguments"); + Py_INCREF(Py_None); + return Py_None; + } + + if (PyArray_TYPE(points) != NPY_DOUBLE) + { + PyErr_SetString(PyExc_ValueError, "invalid function arguments - should be dtype == double"); + Py_INCREF(Py_None); + return Py_None; + } + + if (PyArray_NDIM(points) != 2) + { + PyErr_SetString(PyExc_ValueError, "invalid function arguments - should be len(shape) == 2"); + Py_INCREF(Py_None); + return Py_None; + } + + npy_intp *dims = PyArray_SHAPE(points); + + if (dims[1] != 2) + { + PyErr_SetString(PyExc_ValueError, "invalid function arguments - should be shape[1] == 2"); + Py_INCREF(Py_None); + return Py_None; + } + + size_t num = dims[0]; + + const double *posi = PyArray_DATA(points); + + // We will use C functions directly for known projection types + // Otherwise we will use python call + + struct ProjectionDef in_proj = + { + .projection = input_proj, + .forward = generic_forward_project, + .reverse = generic_reverse_project, + }; + + struct ProjectionDef out_proj = + { + .projection = output_proj, + .forward = generic_forward_project, + .reverse = generic_reverse_project, + }; +/* + if (is_perspective(input_proj)) + { + in_proj.projection = &(((struct PerspectiveProjectionObject *)input_proj)->proj); + in_proj.forward = perspective_projection_project; + in_proj.reverse = perspective_projection_reverse; + } + + if (is_perspective(output_proj)) + { + out_proj.projection = &(((struct PerspectiveProjectionObject *)output_proj)->proj); + out_proj.forward = perspective_projection_project; + out_proj.reverse = perspective_projection_reverse; + } + */ + + PyArrayObject *output_points = (PyArrayObject *)PyArray_ZEROS(2, dims, NPY_DOUBLE, 0); + if (output_points == NULL) + { + PyErr_SetString(PyExc_ValueError, "can not allocate memory"); + Py_INCREF(Py_None); + return Py_None; + } + double *poso = PyArray_DATA(output_points); + + fun(&self->mov, posi, poso, num, &in_proj, &out_proj); + return (PyObject *)output_points; +} + +static PyObject *Sphere_forward(PyObject *_self, + PyObject *args, + PyObject *kwds) +{ + return apply_action(_self, args, kwds, sphere_movement_apply_forward); +} + +static PyObject *Sphere_reverse(PyObject *_self, + PyObject *args, + PyObject *kwds) +{ + return apply_action(_self, args, kwds, sphere_movement_apply_reverse); +} + +static PyMethodDef _SphereMovements_methods[] = { + {"apply_forward", (PyCFunction)Sphere_forward, METH_VARARGS | METH_KEYWORDS, + "Apply forward rotation"}, + {"apply_reverse", (PyCFunction)Sphere_reverse, METH_VARARGS | METH_KEYWORDS, + "Apply reverse rotation"}, + {NULL} /* Sentinel */ +}; + +static PyTypeObject _SphereMovement = { + PyVarObject_HEAD_INIT(NULL, 0) + .tp_name = BASENAME ".SphereMovement", + .tp_doc = PyDoc_STR("Sphere movement object"), + .tp_basicsize = sizeof(struct SphereMovementObject), + .tp_itemsize = 0, + .tp_flags = Py_TPFLAGS_DEFAULT, + .tp_new = PyType_GenericNew, + .tp_init = SphereMovements_init, + .tp_methods = _SphereMovements_methods, +}; + +static PyModuleDef movementsModule = { + PyModuleDef_HEAD_INIT, + .m_name = BASENAME, + .m_doc = "Movements module", + .m_size = -1, +}; + + +PyMODINIT_FUNC +PyInit_movements(void) +{ + import_array(); + PyObject *m = PyModule_Create(&movementsModule); + if (m == NULL) + return NULL; + + if (PyType_Ready(&_SphereMovement) < 0) + { + fprintf(stderr, "Bad type _SphereMovement\n"); + return NULL; + } + + Py_INCREF(&_SphereMovement); + if (PyModule_AddObject(m, "SphereMovement", (PyObject *)&_SphereMovement) < 0) + { + Py_DECREF(&_SphereMovement); + Py_DECREF(m); + fprintf(stderr, "Can not add SphereMovement\n"); + return NULL; + } + + return m; +} diff --git a/src/vstarstack/library/movement/sphere.py b/src/vstarstack/library/movement/sphere.py index 6554d41a..31e16d9b 100644 --- a/src/vstarstack/library/movement/sphere.py +++ b/src/vstarstack/library/movement/sphere.py @@ -20,11 +20,12 @@ from scipy.spatial.transform import Rotation import vstarstack.library.movement.basic_movement +from vstarstack.library.movement.movements import SphereMovement def p2vec(pos): """Lon,lat -> x,y,z""" - lat = pos[0] - lon = pos[1] + lon = pos[0] + lat = pos[1] return np.array([math.cos(lon)*math.cos(lat), math.sin(lon)*math.cos(lat), math.sin(lat)]) @@ -48,41 +49,19 @@ def vecangle(vec1, vec2): class Movement(vstarstack.library.movement.basic_movement.Movement): """Class of spherical movements""" - def apply(self, positions, proj): + def apply(self, positions : np.ndarray, input_proj, output_proj): """Apply movement""" - points = [] - for y, x in positions: - lat, lon = proj.project(y, x) - point = p2vec((lat, lon)) - points.append(point) - points = self.rot.apply(points) - new_positions = [] - for point in points: - lon = math.atan2(point[1], point[0]) - lat = math.asin(point[2]) - y, x = proj.reverse(lat, lon) - new_positions.append((y, x)) - return new_positions - - def reverse(self, positions, proj): + return self.mov.apply_forward(positions, input_proj, output_proj) + + def reverse(self, positions : np.ndarray, input_proj, output_proj): """Apply reverse movement""" - points = [] - for y, x in positions: - lat, lon = proj.project(y, x) - point = p2vec((lat, lon)) - points.append(point) - points = self.rev.apply(points) - new_positions = [] - for point in points: - lon = math.atan2(point[1], point[0]) - lat = math.asin(point[2]) - y, x = proj.reverse(lat, lon) - new_positions.append((y, x)) - return new_positions + return self.mov.apply_reverse(positions, input_proj, output_proj) def __init__(self, rot): self.rot = rot self.rev = rot.inv() + q = self.rot.as_quat() + self.mov = SphereMovement(q[3], q[0], q[1], q[2]) def magnitude(self): """Magnitude of movement""" diff --git a/src/vstarstack/library/projection/__init__.py b/src/vstarstack/library/projection/__init__.py deleted file mode 100644 index e69de29b..00000000 diff --git a/src/vstarstack/library/projection/equirectangular.c b/src/vstarstack/library/projection/equirectangular.c deleted file mode 100644 index 7ea8ce14..00000000 --- a/src/vstarstack/library/projection/equirectangular.c +++ /dev/null @@ -1,117 +0,0 @@ -/* - * Copyright (c) 2022 Vladislav Tsendrovskii - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, version 3 of the License. - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. - * See the GNU General Public License for more details. - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ - -#define PY_SSIZE_T_CLEAN -#include -#include -#include -#include - -struct ProjectionObject -{ - PyObject_HEAD - int h; // height of image in pixels - int w; // width of image in pixels -}; - -// arguments: w, h -static int Projection_init(PyObject *_self, PyObject *args, PyObject *kwds) -{ - struct ProjectionObject *self = (struct ProjectionObject *)_self; - static char *kwlist[] = {"w", "h", NULL}; - if (!PyArg_ParseTupleAndKeywords(args, kwds, "ii", kwlist, - &self->w, &self->h)) - return -1; - if (self->h <= 0 || self->w <= 0) - return -1; - - return 0; -} - -static PyObject *Projection_project(PyObject *_self, PyObject *args, PyObject *kwds) -{ - struct ProjectionObject *self = (struct ProjectionObject *)_self; - double x, y; - static char *kwlist[] = {"y", "x", NULL}; - if (!PyArg_ParseTupleAndKeywords(args, kwds, "dd", kwlist, - &y, &x)) - return Py_None; - double lon = (1 - 2*x/self->w) * M_PI; - double lat = (1 - 2*y/self->h) * M_PI_2; - - return Py_BuildValue("(dd)", lat, lon); -} - -static PyObject *Projection_reverse(PyObject *_self, PyObject *args, PyObject *kwds) -{ - struct ProjectionObject *self = (struct ProjectionObject *)_self; - double lat, lon; - static char *kwlist[] = {"lat", "lon", NULL}; - if (!PyArg_ParseTupleAndKeywords(args, kwds, "dd", kwlist, - &lat, &lon)) - return Py_None; - - double x = (1 - lon / M_PI)/2*self->w; - double y = (1 - lat / M_PI_2)/2*self->h; - return Py_BuildValue("(dd)", y, x); -} - -static PyMethodDef Projection_methods[] = { - {"project", (PyCFunction)Projection_project, METH_VARARGS | METH_KEYWORDS, - "Project y,x to lat,lon"}, - {"reverse", (PyCFunction)Projection_reverse, METH_VARARGS | METH_KEYWORDS, - "Project lat,lon to y,x"}, - {NULL} /* Sentinel */ -}; - -static PyTypeObject Projection = { - PyVarObject_HEAD_INIT(NULL, 0) - .tp_name = "vstarstack.library.projection.equirectangular.Projection", - .tp_doc = PyDoc_STR("Equirectangular projection object"), - .tp_basicsize = sizeof(struct ProjectionObject), - .tp_itemsize = 0, - .tp_flags = Py_TPFLAGS_DEFAULT, - .tp_new = PyType_GenericNew, - .tp_init = Projection_init, - .tp_methods = Projection_methods, -}; - -static PyModuleDef equirectangularModule = { - PyModuleDef_HEAD_INIT, - .m_name = "vstarstack.library.projection.equirectangular", - .m_doc = "Equirectangular projection module", - .m_size = -1, -}; - -PyMODINIT_FUNC -PyInit_equirectangular(void) -{ - PyObject *m; - if (PyType_Ready(&Projection) < 0) - return NULL; - - m = PyModule_Create(&equirectangularModule); - if (m == NULL) - return NULL; - - Py_INCREF(&Projection); - if (PyModule_AddObject(m, "Projection", (PyObject *)&Projection) < 0) - { - Py_DECREF(&Projection); - Py_DECREF(m); - return NULL; - } - - return m; -} diff --git a/src/vstarstack/library/projection/orthographic.c b/src/vstarstack/library/projection/orthographic.c deleted file mode 100644 index 8e99f77d..00000000 --- a/src/vstarstack/library/projection/orthographic.c +++ /dev/null @@ -1,144 +0,0 @@ -/* - * Copyright (c) 2022 Vladislav Tsendrovskii - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, version 3 of the License. - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. - * See the GNU General Public License for more details. - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ - -#define PY_SSIZE_T_CLEAN -#include -#include -#include -#include - -struct ProjectionObject -{ - PyObject_HEAD - int w; // image width - int h; // image height - double a; // planet ellipse major axis - double b; // planet ellipse minor axis - double angle; // planet ellipse slope - double rot; // planet rotation angle -}; - -// arguments: W, H, F, w, h -static int Projection_init(PyObject *_self, PyObject *args, PyObject *kwds) -{ - struct ProjectionObject *self = (struct ProjectionObject *)_self; - static char *kwlist[] = {"w", "h", "a", "b", "angle", "rot", NULL}; - if (!PyArg_ParseTupleAndKeywords(args, kwds, "iidddd", kwlist, - &self->w, &self->h, &self->a, &self->b, - &self->angle, &self->rot)) - return -1; - if (self->h <= 0 || self->w <= 0 || self->a <= 0 || self->b <= 0) - return -1; - - return 0; -} - -static PyObject *Projection_project(PyObject *_self, PyObject *args, PyObject *kwds) -{ - struct ProjectionObject *self = (struct ProjectionObject *)_self; - double x, y; - static char *kwlist[] = {"y", "x", NULL}; - if (!PyArg_ParseTupleAndKeywords(args, kwds, "dd", kwlist, - &y, &x)) - return Py_None; - - x -= self->w/2; - y -= self->h/2; - double X = x*cos(self->angle) - y*sin(self->angle); - double Y = x*sin(self->angle) + y*cos(self->angle); - - double sin_lat = Y / (self->b/2); - if (fabs(sin_lat) > 1) - { - Py_INCREF(Py_None); - Py_INCREF(Py_None); - return Py_BuildValue("(OO)", Py_None, Py_None); - } - double lat = -asin(sin_lat); - double sin_lon = X / (self->a/2) / cos(lat); - if (fabs(sin_lon) > 1) - { - Py_INCREF(Py_None); - Py_INCREF(Py_None); - return Py_BuildValue("(OO)", Py_None, Py_None); - } - - double lon = asin(sin_lon) + self->rot; - return Py_BuildValue("(dd)", lat, lon); -} - -static PyObject *Projection_reverse(PyObject *_self, PyObject *args, PyObject *kwds) -{ - struct ProjectionObject *self = (struct ProjectionObject *)_self; - double lat, lon; - static char *kwlist[] = {"lat", "lon", NULL}; - if (!PyArg_ParseTupleAndKeywords(args, kwds, "dd", kwlist, - &lat, &lon)) - return Py_None; - - double x = self->a/2 * sin(lon - self->rot) * cos(lat); - double z = -self->b/2 * sin(lat); - double X = x*cos(self->angle) + z*sin(self->angle) + self->w/2; - double Y = -x*sin(self->angle) + z*cos(self->angle) + self->h/2; - return Py_BuildValue("(dd)", Y, X); -} - -static PyMethodDef Projection_methods[] = { - {"project", (PyCFunction)Projection_project, METH_VARARGS | METH_KEYWORDS, - "Project y,x to lat,lon"}, - {"reverse", (PyCFunction)Projection_reverse, METH_VARARGS | METH_KEYWORDS, - "Project lat,lon to y,x"}, - {NULL} /* Sentinel */ -}; - -static PyTypeObject Projection = { - PyVarObject_HEAD_INIT(NULL, 0) - .tp_name = "vstarstack.library.projection.orthographic.Projection", - .tp_doc = PyDoc_STR("Orthographic projection object"), - .tp_basicsize = sizeof(struct ProjectionObject), - .tp_itemsize = 0, - .tp_flags = Py_TPFLAGS_DEFAULT, - .tp_new = PyType_GenericNew, - .tp_init = Projection_init, - .tp_methods = Projection_methods, -}; - -static PyModuleDef orthographicModule = { - PyModuleDef_HEAD_INIT, - .m_name = "vstarstack.library.projection.orthographic", - .m_doc = "Orthographic projection module", - .m_size = -1, -}; - -PyMODINIT_FUNC -PyInit_orthographic(void) -{ - PyObject *m; - if (PyType_Ready(&Projection) < 0) - return NULL; - - m = PyModule_Create(&orthographicModule); - if (m == NULL) - return NULL; - - Py_INCREF(&Projection); - if (PyModule_AddObject(m, "Projection", (PyObject *)&Projection) < 0) - { - Py_DECREF(&Projection); - Py_DECREF(m); - return NULL; - } - - return m; -} diff --git a/src/vstarstack/library/projection/perspective.c b/src/vstarstack/library/projection/perspective.c deleted file mode 100644 index fed4c4a9..00000000 --- a/src/vstarstack/library/projection/perspective.c +++ /dev/null @@ -1,128 +0,0 @@ -/* - * Copyright (c) 2022 Vladislav Tsendrovskii - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, version 3 of the License. - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. - * See the GNU General Public License for more details. - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ - -#define PY_SSIZE_T_CLEAN -#include -#include -#include -#include - -struct ProjectionObject -{ - PyObject_HEAD - double H; // height of image in mm - double W; // width of image in mm - double F; // focal length in mm - int h; // height of image in pixels - int w; // width of image in pixels - double kx; // transformation from pixels to mm coefficient - double ky; // transformation from pixels to mm coefficient -}; - -// arguments: W, H, F, w, h -static int Projection_init(PyObject *_self, PyObject *args, PyObject *kwds) -{ - struct ProjectionObject *self = (struct ProjectionObject *)_self; - static char *kwlist[] = {"W", "H", "F", "w", "h", NULL}; - if (!PyArg_ParseTupleAndKeywords(args, kwds, "dddii", kwlist, - &self->W, &self->H, &self->F, &self->w, &self->h)) - return -1; - if (self->h <= 0 || self->w <= 0 || self->W <= 0 || self->H <= 0 || self->F <= 0) - return -1; - - self->kx = self->W / self->w; - self->ky = self->H / self->h; - return 0; -} - -static PyObject *Projection_project(PyObject *_self, PyObject *args, PyObject *kwds) -{ - struct ProjectionObject *self = (struct ProjectionObject *)_self; - double x, y; - static char *kwlist[] = {"y", "x", NULL}; - if (!PyArg_ParseTupleAndKeywords(args, kwds, "dd", kwlist, - &y, &x)) - return Py_None; - double X = (self->w / 2 - x) * self->kx; - double Y = (self->h / 2 - y) * self->ky; - double lon = atan(X / self->F); - double lat = atan(Y * cos(lon) / self->F); - - return Py_BuildValue("(dd)", lat, lon); -} - -static PyObject *Projection_reverse(PyObject *_self, PyObject *args, PyObject *kwds) -{ - struct ProjectionObject *self = (struct ProjectionObject *)_self; - double lat, lon; - static char *kwlist[] = {"lat", "lon", NULL}; - if (!PyArg_ParseTupleAndKeywords(args, kwds, "dd", kwlist, - &lat, &lon)) - return Py_None; - - double X = self->F * tan(lon); - double Y = self->F * tan(lat) / cos(lon); - double x = (self->w / 2 - X / self->kx); - double y = (self->h / 2 - Y / self->ky); - return Py_BuildValue("(dd)", y, x); -} - -static PyMethodDef Projection_methods[] = { - {"project", (PyCFunction)Projection_project, METH_VARARGS | METH_KEYWORDS, - "Project y,x to lat,lon"}, - {"reverse", (PyCFunction)Projection_reverse, METH_VARARGS | METH_KEYWORDS, - "Project lat,lon to y,x"}, - {NULL} /* Sentinel */ -}; - -static PyTypeObject Projection = { - PyVarObject_HEAD_INIT(NULL, 0) - .tp_name = "vstarstack.library.projection.perspective.Projection", - .tp_doc = PyDoc_STR("Perspective projection object"), - .tp_basicsize = sizeof(struct ProjectionObject), - .tp_itemsize = 0, - .tp_flags = Py_TPFLAGS_DEFAULT, - .tp_new = PyType_GenericNew, - .tp_init = Projection_init, - .tp_methods = Projection_methods, -}; - -static PyModuleDef orthographicModule = { - PyModuleDef_HEAD_INIT, - .m_name = "vstarstack.library.projection.perspective", - .m_doc = "Perspective projection module", - .m_size = -1, -}; - -PyMODINIT_FUNC -PyInit_perspective(void) -{ - PyObject *m; - if (PyType_Ready(&Projection) < 0) - return NULL; - - m = PyModule_Create(&orthographicModule); - if (m == NULL) - return NULL; - - Py_INCREF(&Projection); - if (PyModule_AddObject(m, "Projection", (PyObject *)&Projection) < 0) - { - Py_DECREF(&Projection); - Py_DECREF(m); - return NULL; - } - - return m; -} diff --git a/src/vstarstack/library/projection/projections/lib/equirectangular.c b/src/vstarstack/library/projection/projections/lib/equirectangular.c new file mode 100644 index 00000000..f4f67d3c --- /dev/null +++ b/src/vstarstack/library/projection/projections/lib/equirectangular.c @@ -0,0 +1,47 @@ +/* + * Copyright (c) 2023 Vladislav Tsendrovskii + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, version 3 of the License. + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU General Public License for more details. + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include + +#include "equirectangular.h" + +bool equirectangular_projection_init(struct EquirectangularProjection *self, int w, int h) +{ + if (h <= 0 || w <= 0) + return false; + + self->w = w; + self->h = h; + return true; +} + +bool equirectangular_projection_project(void *self, + double y, double x, + double *lat, double *lon) +{ + const struct EquirectangularProjection *proj = (const struct EquirectangularProjection *)self; + *lon = (1 - 2*x/proj->w) * M_PI; + *lat = (1 - 2*y/proj->h) * M_PI_2; + return true; +} + +bool equirectangular_projection_reverse(void *self, + double lat, double lon, + double *y, double *x) +{ + const struct EquirectangularProjection *proj = (const struct EquirectangularProjection *)self; + *x = (1 - lon / M_PI)/2*proj->w; + *y = (1 - lat / M_PI_2)/2*proj->h; + return true; +} diff --git a/src/vstarstack/library/projection/projections/lib/equirectangular.h b/src/vstarstack/library/projection/projections/lib/equirectangular.h new file mode 100644 index 00000000..df985ed1 --- /dev/null +++ b/src/vstarstack/library/projection/projections/lib/equirectangular.h @@ -0,0 +1,33 @@ +/* + * Copyright (c) 2023 Vladislav Tsendrovskii + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, version 3 of the License. + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU General Public License for more details. + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#pragma once + +#include "projection_common.h" + +struct EquirectangularProjection +{ + int h; // height of image in pixels + int w; // width of image in pixels +}; + +bool equirectangular_projection_init(struct EquirectangularProjection *self, int w, int h); + +bool equirectangular_projection_project(void *self, + double y, double x, + double *lat, double *lon); + +bool equirectangular_projection_reverse(void *self, + double lat, double lon, + double *y, double *x); diff --git a/src/vstarstack/library/projection/projections/lib/orthographic.c b/src/vstarstack/library/projection/projections/lib/orthographic.c new file mode 100644 index 00000000..780f2836 --- /dev/null +++ b/src/vstarstack/library/projection/projections/lib/orthographic.c @@ -0,0 +1,82 @@ +/* + * Copyright (c) 2023 Vladislav Tsendrovskii + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, version 3 of the License. + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU General Public License for more details. + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include + +#include "orthographic.h" +#include + +bool orthographic_projection_init(struct OrthographicProjection *self, + int w, int h, + double a, double b, + double angle, double rot) +{ + if (h <= 0 || w <= 0 || a <= 0 || b <= 0) + return false; + + self->a = a; + self->b = b; + self->angle = angle; + self->rot = rot; + self->w = w; + self->h = h; + return true; +} + +bool orthographic_projection_project(void *self, + double y, double x, + double *lat, double *lon) +{ + const struct OrthographicProjection *proj = (const struct OrthographicProjection *)self; + + x -= proj->w/2; + y -= proj->h/2; + double X = x*cos(proj->angle) - y*sin(proj->angle); + double Y = x*sin(proj->angle) + y*cos(proj->angle); + + double sin_lat = Y / (proj->b/2); + if (fabs(sin_lat) > 1) + return false; + + *lat = -asin(sin_lat); + if (fabs(cos(*lat)) < 1e-6) + { + if (fabs(X) < 1e-4) + { + *lon = 0; + return true; + } + return false; + } + + double sin_lon = X / (proj->a/2) / cos(*lat); + if (fabs(sin_lon) > 1) + return false; + + *lon = asin(sin_lon) + proj->rot; + return true; +} + +bool orthographic_projection_reverse(void *self, + double lat, double lon, + double *y, double *x) +{ + const struct OrthographicProjection *proj = (const struct OrthographicProjection *)self; + + double X = proj->a/2 * sin(lon - proj->rot) * cos(lat); + double Z = -proj->b/2 * sin(lat); + *x = X*cos(proj->angle) + Z*sin(proj->angle) + proj->w/2; + *y = -X*sin(proj->angle) + Z*cos(proj->angle) + proj->h/2; + return true; +} diff --git a/src/vstarstack/library/projection/projections/lib/orthographic.h b/src/vstarstack/library/projection/projections/lib/orthographic.h new file mode 100644 index 00000000..c618cb62 --- /dev/null +++ b/src/vstarstack/library/projection/projections/lib/orthographic.h @@ -0,0 +1,40 @@ +/* + * Copyright (c) 2023 Vladislav Tsendrovskii + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, version 3 of the License. + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU General Public License for more details. + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#pragma once + +#include "projection_common.h" + +struct OrthographicProjection +{ + int w; // image width + int h; // image height + double a; // planet ellipse major axis + double b; // planet ellipse minor axis + double angle; // planet ellipse slope + double rot; // planet rotation angle +}; + +bool orthographic_projection_init(struct OrthographicProjection *self, + int w, int h, + double a, double b, + double angle, double rot); + +bool orthographic_projection_project(void *self, + double y, double x, + double *lat, double *lon); + +bool orthographic_projection_reverse(void *self, + double lat, double lon, + double *y, double *x); diff --git a/src/vstarstack/library/projection/projections/lib/perspective.c b/src/vstarstack/library/projection/projections/lib/perspective.c new file mode 100644 index 00000000..94d404c5 --- /dev/null +++ b/src/vstarstack/library/projection/projections/lib/perspective.c @@ -0,0 +1,64 @@ +/* + * Copyright (c) 2023 Vladislav Tsendrovskii + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, version 3 of the License. + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU General Public License for more details. + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include + +#include "perspective.h" + +bool perspective_projection_init(struct PerspectiveProjection *self, + double W, double H, double F, + double w, double h) +{ + if (h <= 0 || w <= 0 || + W <= 0 || H <= 0 || + F <= 0) + { + return false; + } + + self->W = W; + self->H = H; + self->F = F; + self->w = w; + self->h = h; + + self->kx = self->W / self->w; + self->ky = self->H / self->h; + return true; +} + +bool perspective_projection_project(void *self, + double y, double x, + double *lat, double *lon) +{ + const struct PerspectiveProjection *proj = (const struct PerspectiveProjection *)self; + double X = (proj->w / 2 - x) * proj->kx; + double Y = (proj->h / 2 - y) * proj->ky; + *lon = atan(X / proj->F); + *lat = atan(Y * cos(*lon) / proj->F); + return true; +} + +bool perspective_projection_reverse(void *self, + double lat, double lon, + double *y, double *x) +{ + const struct PerspectiveProjection *proj = (const struct PerspectiveProjection *)self; + double X = proj->F * tan(lon); + double Y = proj->F * tan(lat) / cos(lon); + *x = proj->w / 2 - X / proj->kx; + *y = proj->h / 2 - Y / proj->ky; + return true; +} diff --git a/src/vstarstack/library/projection/projections/lib/perspective.h b/src/vstarstack/library/projection/projections/lib/perspective.h new file mode 100644 index 00000000..d1bd0d64 --- /dev/null +++ b/src/vstarstack/library/projection/projections/lib/perspective.h @@ -0,0 +1,40 @@ +/* + * Copyright (c) 2023 Vladislav Tsendrovskii + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, version 3 of the License. + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU General Public License for more details. + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#pragma once + +#include "projection_common.h" + +struct PerspectiveProjection +{ + double H; // height of image in mm + double W; // width of image in mm + double F; // focal length in mm + int h; // height of image in pixels + int w; // width of image in pixels + double kx; // transformation from pixels to mm coefficient + double ky; // transformation from pixels to mm coefficient +}; + +bool perspective_projection_init(struct PerspectiveProjection *self, + double W, double H, double F, + double w, double h); + +bool perspective_projection_project(void *self, + double y, double x, + double *lat, double *lon); + +bool perspective_projection_reverse(void *self, + double lat, double lon, + double *y, double *x); diff --git a/src/vstarstack/library/projection/projections/lib/projection.h b/src/vstarstack/library/projection/projections/lib/projection.h new file mode 100644 index 00000000..57b41708 --- /dev/null +++ b/src/vstarstack/library/projection/projections/lib/projection.h @@ -0,0 +1,23 @@ +/* + * Copyright (c) 2023 Vladislav Tsendrovskii + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, version 3 of the License. + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU General Public License for more details. + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#pragma once + +#include + +#include "projection_common.h" + +#include "perspective.h" +#include "orthographic.h" +#include "equirectangular.h" diff --git a/src/vstarstack/library/projection/projections/lib/projection_common.h b/src/vstarstack/library/projection/projections/lib/projection_common.h new file mode 100644 index 00000000..82ce2cd2 --- /dev/null +++ b/src/vstarstack/library/projection/projections/lib/projection_common.h @@ -0,0 +1,25 @@ +/* + * Copyright (c) 2023 Vladislav Tsendrovskii + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, version 3 of the License. + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU General Public License for more details. + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#pragma once + +#include + +typedef bool (*forward_project_f)(void *self, + double y, double x, + double *lat, double *lon); + +typedef bool (*reverse_project_f)(void *self, + double lat, double lon, + double *y, double *x); diff --git a/src/vstarstack/library/projection/projections/module.c b/src/vstarstack/library/projection/projections/module.c new file mode 100644 index 00000000..dd111a89 --- /dev/null +++ b/src/vstarstack/library/projection/projections/module.c @@ -0,0 +1,291 @@ +/* + * Copyright (c) 2023 Vladislav Tsendrovskii + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, version 3 of the License. + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU General Public License for more details. + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include + +#include "projection_module.h" + +#define BASENAME "vstarstack.library.projection.projections" + +// arguments: W, H, F, w, h +static int PerspectiveProjection_init(PyObject *_self, PyObject *args, PyObject *kwds) +{ + double W, H, F; + int w, h; + struct PerspectiveProjectionObject *self = (struct PerspectiveProjectionObject *)_self; + static char *kwlist[] = {"w", "h", "W", "H", "F", NULL}; + if (!PyArg_ParseTupleAndKeywords(args, kwds, "iiddd", kwlist, + &w, &h, + &W, &H, &F)) + return -1; + + if (!perspective_projection_init(&self->proj, W, H, F, w, h)) + return -1; + + return 0; +} + +static int OrthographicProjection_init(PyObject *_self, PyObject *args, PyObject *kwds) +{ + double a, b, angle, rot; + int w, h; + struct OrthographicProjectionObject *self = (struct OrthographicProjectionObject *)_self; + static char *kwlist[] = {"w", "h", "a", "b", "angle", "rot", NULL}; + if (!PyArg_ParseTupleAndKeywords(args, kwds, "iidddd", kwlist, + &w, &h, + &a, &b, &angle, &rot)) + return -1; + + if (!orthographic_projection_init(&self->proj, w, h, a, b, angle, rot)) + return -1; + + return 0; +} + +static int EquirectangularProjection_init(PyObject *_self, PyObject *args, PyObject *kwds) +{ + int w, h; + struct EquirectangularProjectionObject *self = (struct EquirectangularProjectionObject *)_self; + static char *kwlist[] = {"w", "h", NULL}; + if (!PyArg_ParseTupleAndKeywords(args, kwds, "ii", kwlist, + &w, &h)) + return -1; + + if (!equirectangular_projection_init(&self->proj, w, h)) + return -1; + + return 0; +} + +// Common methods + +static PyObject *Projection_project(void *proj, + PyObject *args, + PyObject *kwds, + forward_project_f function) +{ + double x, y; + double lon, lat; + static char *kwlist[] = {"x", "y", NULL}; + if (!PyArg_ParseTupleAndKeywords(args, kwds, "dd", kwlist, + &x, &y)) + return Py_BuildValue("(OO)", Py_None, Py_None); + + if (!function(proj, y, x, &lat, &lon)) + return Py_BuildValue("(OO)", Py_None, Py_None); + return Py_BuildValue("(dd)", lon, lat); +} + +static PyObject *Projection_reverse(void *proj, + PyObject *args, + PyObject *kwds, + reverse_project_f function) +{ + double x, y; + double lon, lat; + static char *kwlist[] = {"lon", "lat", NULL}; + if (!PyArg_ParseTupleAndKeywords(args, kwds, "dd", kwlist, + &lon, &lat)) + return Py_BuildValue("(OO)", Py_None, Py_None); + + //printf("lat %lf lon %lf\n", lat, lon); + if (!function(proj, lat, lon, &y, &x)) + return Py_BuildValue("(OO)", Py_None, Py_None); + //printf("x %lf y %lf\n", x, y); + return Py_BuildValue("(dd)", x, y); +} + +// Perspective + +static PyObject *Perspective_forward(PyObject *_self, + PyObject *args, + PyObject *kwds) +{ + struct PerspectiveProjectionObject *self = (struct PerspectiveProjectionObject *)_self; + return Projection_project(&self->proj, args, kwds, perspective_projection_project); +} + + +static PyObject *Perspective_reverse(PyObject *_self, + PyObject *args, + PyObject *kwds) +{ + struct PerspectiveProjectionObject *self = (struct PerspectiveProjectionObject *)_self; + return Projection_reverse(&self->proj, args, kwds, perspective_projection_reverse); +} + +static PyMethodDef _PerspectiveProjection_methods[] = { + {"project", (PyCFunction)Perspective_forward, METH_VARARGS | METH_KEYWORDS, + "Project x,y to lat,lon"}, + {"reverse", (PyCFunction)Perspective_reverse, METH_VARARGS | METH_KEYWORDS, + "Project lon,lat to x,y"}, + {NULL} /* Sentinel */ +}; + +static PyTypeObject PerspectiveProjection = { + PyVarObject_HEAD_INIT(NULL, 0) + .tp_name = BASENAME ".PerspectiveProjection", + .tp_doc = PyDoc_STR("Perspective projection object"), + .tp_basicsize = sizeof(struct PerspectiveProjectionObject), + .tp_itemsize = 0, + .tp_flags = Py_TPFLAGS_DEFAULT, + .tp_new = PyType_GenericNew, + .tp_init = PerspectiveProjection_init, + .tp_methods = _PerspectiveProjection_methods, +}; + + +// Orthographic + +static PyObject *Orthographic_forward(PyObject *_self, + PyObject *args, + PyObject *kwds) +{ + struct OrthographicProjectionObject *self = (struct OrthographicProjectionObject *)_self; + return Projection_project(&self->proj, args, kwds, orthographic_projection_project); +} + +static PyObject *Orthographic_reverse(PyObject *_self, + PyObject *args, + PyObject *kwds) +{ + struct OrthographicProjectionObject *self = (struct OrthographicProjectionObject *)_self; + return Projection_reverse(&self->proj, args, kwds, orthographic_projection_reverse); +} + +static PyMethodDef _OrthographicProjection_methods[] = { + {"project", (PyCFunction)Orthographic_forward, METH_VARARGS | METH_KEYWORDS, + "Project x,y to lon,lat"}, + {"reverse", (PyCFunction)Orthographic_reverse, METH_VARARGS | METH_KEYWORDS, + "Project lon,lat to x,y"}, + {NULL} /* Sentinel */ +}; + +static PyTypeObject OrthographicProjection = { + PyVarObject_HEAD_INIT(NULL, 0) + .tp_name = BASENAME ".OrthographicProjection", + .tp_doc = PyDoc_STR("Orthographic projection object"), + .tp_basicsize = sizeof(struct OrthographicProjectionObject), + .tp_itemsize = 0, + .tp_flags = Py_TPFLAGS_DEFAULT, + .tp_new = PyType_GenericNew, + .tp_init = OrthographicProjection_init, + .tp_methods = _OrthographicProjection_methods, +}; + + +// Equirectangular + +static PyObject *Equirectangular_forward(PyObject *_self, + PyObject *args, + PyObject *kwds) +{ + struct EquirectangularProjectionObject *self = (struct EquirectangularProjectionObject *)_self; + return Projection_project(&self->proj, args, kwds, equirectangular_projection_project); +} + +static PyObject *Equirectangular_reverse(PyObject *_self, + PyObject *args, + PyObject *kwds) +{ + struct EquirectangularProjectionObject *self = (struct EquirectangularProjectionObject *)_self; + return Projection_reverse(&self->proj, args, kwds, equirectangular_projection_reverse); +} + +static PyMethodDef _EquirectangularProjection_methods[] = { + {"project", (PyCFunction)Equirectangular_forward, METH_VARARGS | METH_KEYWORDS, + "Project x,y to lon,lat"}, + {"reverse", (PyCFunction)Equirectangular_reverse, METH_VARARGS | METH_KEYWORDS, + "Project lon,lat to x,y"}, + {NULL} /* Sentinel */ +}; + +static PyTypeObject EquirectangularProjection = { + PyVarObject_HEAD_INIT(NULL, 0) + .tp_name = BASENAME ".EquirectangularProjection", + .tp_doc = PyDoc_STR("Equirectangular projection object"), + .tp_basicsize = sizeof(struct EquirectangularProjectionObject), + .tp_itemsize = 0, + .tp_flags = Py_TPFLAGS_DEFAULT, + .tp_new = PyType_GenericNew, + .tp_init = EquirectangularProjection_init, + .tp_methods = _EquirectangularProjection_methods, +}; + +// Module + +static PyModuleDef projectionModule = { + PyModuleDef_HEAD_INIT, + .m_name = BASENAME, + .m_doc = "Projection module", + .m_size = -1, +}; + +PyMODINIT_FUNC +PyInit_projections(void) +{ + PyObject *m = PyModule_Create(&projectionModule); + if (m == NULL) + return NULL; + + if (PyType_Ready(&PerspectiveProjection) < 0) + { + fprintf(stderr, "Bad type _PerspectiveProjection\n"); + return NULL; + } + + if (PyType_Ready(&OrthographicProjection) < 0) + { + fprintf(stderr, "Bad type _OrthographicProjection\n"); + return NULL; + } + + if (PyType_Ready(&EquirectangularProjection) < 0) + { + fprintf(stderr, "Bad type _EquirectangularProjection\n"); + return NULL; + } + + Py_INCREF(&PerspectiveProjection); + if (PyModule_AddObject(m, "PerspectiveProjection", (PyObject *)&PerspectiveProjection) < 0) + { + Py_DECREF(&PerspectiveProjection); + Py_DECREF(m); + fprintf(stderr, "Can not add PerspectiveProjection\n"); + return NULL; + } + + Py_INCREF(&OrthographicProjection); + if (PyModule_AddObject(m, "OrthographicProjection", (PyObject *)&OrthographicProjection) < 0) + { + Py_DECREF(&PerspectiveProjection); + Py_DECREF(&OrthographicProjection); + Py_DECREF(m); + fprintf(stderr, "Can not add OrthographicProjection\n"); + return NULL; + } + + Py_INCREF(&EquirectangularProjection); + if (PyModule_AddObject(m, "EquirectangularProjection", (PyObject *)&EquirectangularProjection) < 0) + { + Py_DECREF(&PerspectiveProjection); + Py_DECREF(&OrthographicProjection); + Py_DECREF(&EquirectangularProjection); + Py_DECREF(m); + fprintf(stderr, "Can not add EquirectangularProjection\n"); + return NULL; + } + + return m; +} diff --git a/src/vstarstack/library/projection/projections/projection_module.h b/src/vstarstack/library/projection/projections/projection_module.h new file mode 100644 index 00000000..15a19ade --- /dev/null +++ b/src/vstarstack/library/projection/projections/projection_module.h @@ -0,0 +1,39 @@ +/* + * Copyright (c) 2023 Vladislav Tsendrovskii + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, version 3 of the License. + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU General Public License for more details. + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#pragma once + +#define PY_SSIZE_T_CLEAN +#include +#include + +#include "lib/projection.h" + +struct PerspectiveProjectionObject +{ + PyObject_HEAD + struct PerspectiveProjection proj; +}; + +struct OrthographicProjectionObject +{ + PyObject_HEAD + struct OrthographicProjection proj; +}; + +struct EquirectangularProjectionObject +{ + PyObject_HEAD + struct EquirectangularProjection proj; +}; diff --git a/src/vstarstack/library/projection/tools.py b/src/vstarstack/library/projection/tools.py index 9eb4162b..1c48dde5 100644 --- a/src/vstarstack/library/projection/tools.py +++ b/src/vstarstack/library/projection/tools.py @@ -14,7 +14,7 @@ # import vstarstack.library.data -import vstarstack.library.projection.perspective +import vstarstack.library.projection.projections def add_description(dataframe : vstarstack.library.data.DataFrame, projection : str, **argv): """Add projection description to dataframe""" @@ -27,6 +27,15 @@ def add_description(dataframe : vstarstack.library.data.DataFrame, projection : dataframe.add_parameter(kh, "perspective_kh") elif projection == "equirectangular": pass + elif projection == "orthographic": + a = argv["a"] + b = argv["b"] + angle = argv["angle"] + rot = argv["rot"] + dataframe.add_parameter(a, "orthographic_a") + dataframe.add_parameter(b, "orthographic_b") + dataframe.add_parameter(angle, "orthographic_angle") + dataframe.add_parameter(rot, "orthographic_rot") dataframe.add_parameter(projection, "projection") @@ -40,11 +49,20 @@ def get_projection(dataframe : vstarstack.library.data.DataFrame): h = dataframe.params["h"] W = w * dataframe.params["perspective_kw"] H = h * dataframe.params["perspective_kh"] - return vstarstack.library.projection.perspective.Projection(W, H, F, w, h) + return vstarstack.library.projection.projections.PerspectiveProjection(w, h, W, H, F) if dataframe.params["projection"] == "equirectangular": w = dataframe.params["w"] h = dataframe.params["h"] - return vstarstack.library.projection.equirectangular.Projection(w, h) + return vstarstack.library.projection.projections.EquirectangularProjection(w, h) + + if dataframe.params["projection"] == "orthographic": + w = dataframe.params["w"] + h = dataframe.params["h"] + a = dataframe.params["a"] + b = dataframe.params["b"] + angle = dataframe.params["angle"] + rot = dataframe.params["rot"] + return vstarstack.library.projection.projections.OrthographicProjection(w, h, a, b, angle, rot) raise Exception("Unknown projection") diff --git a/tests/test_equirectangular.py b/tests/test_equirectangular.py index 744f8c9d..d85d074b 100644 --- a/tests/test_equirectangular.py +++ b/tests/test_equirectangular.py @@ -14,7 +14,7 @@ import math -import vstarstack.library.projection.equirectangular +import vstarstack.library.projection.projections thr = 1e-6 @@ -28,8 +28,8 @@ def test_center_f(): lon_expected = 0 lat_expected = 0 - proj = vstarstack.library.projection.equirectangular.Projection(w, h) - lat, lon = proj.project(y, x) + proj = vstarstack.library.projection.projections.EquirectangularProjection(w, h) + lon, lat = proj.project(x, y) assert abs(lat - lat_expected) < thr assert abs(lon - lon_expected) < thr @@ -44,8 +44,8 @@ def test_center_r(): x_expected = 750 y_expected = 500 - proj = vstarstack.library.projection.equirectangular.Projection(w, h) - y, x = proj.reverse(lat, lon) + proj = vstarstack.library.projection.projections.EquirectangularProjection(w, h) + x, y = proj.reverse(lon, lat) assert abs(y - y_expected) < thr assert abs(x - x_expected) < thr @@ -60,8 +60,8 @@ def test_left_f(): lon_expected = math.pi lat_expected = 0 - proj = vstarstack.library.projection.equirectangular.Projection(w, h) - lat, lon = proj.project(y, x) + proj = vstarstack.library.projection.projections.EquirectangularProjection(w, h) + lon, lat = proj.project(x, y) assert abs(lat - lat_expected) < thr assert abs(lon - lon_expected) < thr @@ -76,8 +76,8 @@ def test_left_r(): lon = math.pi lat = 0 - proj = vstarstack.library.projection.equirectangular.Projection(w, h) - y, x = proj.reverse(lat, lon) + proj = vstarstack.library.projection.projections.EquirectangularProjection(w, h) + x, y = proj.reverse(lon, lat) assert abs(y - y_expected) < thr assert abs(x - x_expected) < thr @@ -93,8 +93,8 @@ def test_top_f(): lon_expected = 0 lat_expected = math.pi/2 - proj = vstarstack.library.projection.equirectangular.Projection(w, h) - lat, lon = proj.project(y, x) + proj = vstarstack.library.projection.projections.EquirectangularProjection(w, h) + lon, lat = proj.project(x, y) assert abs(lat - lat_expected) < thr assert abs(lon - lon_expected) < thr @@ -109,8 +109,8 @@ def test_top_r(): lon = 0 lat = math.pi/2 - proj = vstarstack.library.projection.equirectangular.Projection(w, h) - y, x = proj.reverse(lat, lon) + proj = vstarstack.library.projection.projections.EquirectangularProjection(w, h) + x, y = proj.reverse(lon, lat) assert abs(y - y_expected) < thr assert abs(x - x_expected) < thr diff --git a/tests/test_movement_flat.py b/tests/test_movement_flat.py index 68cc4f6a..555e204a 100644 --- a/tests/test_movement_flat.py +++ b/tests/test_movement_flat.py @@ -13,7 +13,6 @@ # import math -import vstarstack.library.projection.perspective as perspective from vstarstack.library.movement import flat thr = 1e-6 @@ -46,21 +45,20 @@ def test_no_movement_forward(): s2y2 = h/2 s2x2 = w/2+1 - s1pos1 = (s1y1, s1x1) - s1pos2 = (s1y2, s1x2) + s1pos1 = (s1x1, s1y1) + s1pos2 = (s1x2, s1y2) - s2pos1 = (s2y1, s2x1) - s2pos2 = (s2y2, s2x2) + s2pos1 = (s2x1, s2y1) + s2pos2 = (s2x2, s2y2) - positions = [(h/2, w/2)] + positions = [(w/2, h/2)] - proj = perspective.Projection(W, H, F, w, h) movement = flat.Movement.build(s1pos1, s2pos1, s1pos2, s2pos2) - shifted = movement.apply(positions, proj) + shifted = movement.apply(positions) assert len(shifted) == 1 - assert abs(shifted[0][0] - h/2) < pixthr - assert abs(shifted[0][1] - w/2) < pixthr + assert abs(shifted[0][0] - w/2) < pixthr + assert abs(shifted[0][1] - h/2) < pixthr def test_no_movement_reverse(): # star1 on frame1 @@ -79,21 +77,20 @@ def test_no_movement_reverse(): s2y2 = h/2 s2x2 = w/2+1 - s1pos1 = (s1y1, s1x1) - s1pos2 = (s1y2, s1x2) + s1pos1 = (s1x1, s1y1) + s1pos2 = (s1x2, s1y2) - s2pos1 = (s2y1, s2x1) - s2pos2 = (s2y2, s2x2) + s2pos1 = (s2x1, s2y1) + s2pos2 = (s2x2, s2y2) - positions = [(h/2, w/2)] + positions = [(w/2, h/2)] - proj = perspective.Projection(W, H, F, w, h) movement = flat.Movement.build(s1pos1, s2pos1, s1pos2, s2pos2) - shifted = movement.reverse(positions, proj) + shifted = movement.reverse(positions) assert len(shifted) == 1 - assert abs(shifted[0][0] - h/2) < pixthr - assert abs(shifted[0][1] - w/2) < pixthr + assert abs(shifted[0][0] - w/2) < pixthr + assert abs(shifted[0][1] - h/2) < pixthr def test_rotation_forward(): # star1 on frame1 @@ -112,25 +109,24 @@ def test_rotation_forward(): s2y2 = h/2+1 s2x2 = w/2 - s1pos1 = (s1y1, s1x1) - s1pos2 = (s1y2, s1x2) + s1pos1 = (s1x1, s1y1) + s1pos2 = (s1x2, s1y2) - s2pos1 = (s2y1, s2x1) - s2pos2 = (s2y2, s2x2) + s2pos1 = (s2x1, s2y1) + s2pos2 = (s2x2, s2y2) - positions = [(h/2, w/2), (h/2, w/2+1), (h/2+1, w/2)] + positions = [(w/2, h/2), (w/2+1, h/2), (w/2, h/2+1)] - proj = perspective.Projection(W, H, F, w, h) movement = flat.Movement.build(s1pos1, s2pos1, s1pos2, s2pos2) - shifted = movement.apply(positions, proj) + shifted = movement.apply(positions) assert len(shifted) == 3 - assert abs(shifted[0][0] - h/2) < pixthr - assert abs(shifted[0][1] - w/2) < pixthr - assert abs(shifted[1][0] - (h/2+1)) < pixthr - assert abs(shifted[1][1] - w/2) < pixthr - assert abs(shifted[2][0] - h/2) < pixthr - assert abs(shifted[2][1] - (w/2-1)) < pixthr + assert abs(shifted[0][1] - h/2) < pixthr + assert abs(shifted[0][0] - w/2) < pixthr + assert abs(shifted[1][1] - (h/2+1)) < pixthr + assert abs(shifted[1][0] - w/2) < pixthr + assert abs(shifted[2][1] - h/2) < pixthr + assert abs(shifted[2][0] - (w/2-1)) < pixthr def test_rotation_reverse(): # star1 on frame1 @@ -149,25 +145,24 @@ def test_rotation_reverse(): s2y2 = h/2+1 s2x2 = w/2 - s1pos1 = (s1y1, s1x1) - s1pos2 = (s1y2, s1x2) + s1pos1 = (s1x1, s1y1) + s1pos2 = (s1x2, s1y2) - s2pos1 = (s2y1, s2x1) - s2pos2 = (s2y2, s2x2) + s2pos1 = (s2x1, s2y1) + s2pos2 = (s2x2, s2y2) - positions = [(h/2, w/2), (h/2, w/2+1), (h/2+1, w/2)] + positions = [(w/2, h/2), (w/2+1, h/2), (w/2, h/2+1)] - proj = perspective.Projection(W, H, F, w, h) movement = flat.Movement.build(s1pos1, s2pos1, s1pos2, s2pos2) - shifted = movement.reverse(positions, proj) + shifted = movement.reverse(positions) assert len(shifted) == 3 - assert abs(shifted[0][0] - h/2) < pixthr - assert abs(shifted[0][1] - w/2) < pixthr - assert abs(shifted[1][0] - (h/2-1)) < pixthr - assert abs(shifted[1][1] - w/2) < pixthr - assert abs(shifted[2][0] - h/2) < pixthr - assert abs(shifted[2][1] - (w/2+1)) < pixthr + assert abs(shifted[0][1] - h/2) < pixthr + assert abs(shifted[0][0] - w/2) < pixthr + assert abs(shifted[1][1] - (h/2-1)) < pixthr + assert abs(shifted[1][0] - w/2) < pixthr + assert abs(shifted[2][1] - h/2) < pixthr + assert abs(shifted[2][0] - (w/2+1)) < pixthr def test_shift_forward(): # star1 on frame1 @@ -186,25 +181,24 @@ def test_shift_forward(): s2y2 = h/2 s2x2 = w/2+2 - s1pos1 = (s1y1, s1x1) - s1pos2 = (s1y2, s1x2) + s1pos1 = (s1x1, s1y1) + s1pos2 = (s1x2, s1y2) - s2pos1 = (s2y1, s2x1) - s2pos2 = (s2y2, s2x2) + s2pos1 = (s2x1, s2y1) + s2pos2 = (s2x2, s2y2) - positions = [(h/2, w/2), (h/2, w/2+1), (h/2+1, w/2)] + positions = [(w/2, h/2), (w/2+1, h/2), (w/2, h/2+1)] - proj = perspective.Projection(W, H, F, w, h) movement = flat.Movement.build(s1pos1, s2pos1, s1pos2, s2pos2) - shifted = movement.apply(positions, proj) + shifted = movement.apply(positions) assert len(shifted) == 3 - assert abs(shifted[0][0] - h/2) < pixthr - assert abs(shifted[0][1] - (w/2+1)) < pixthr - assert abs(shifted[1][0] - h/2) < pixthr - assert abs(shifted[1][1] - (w/2+2)) < pixthr - assert abs(shifted[2][0] - (h/2+1)) < pixthr - assert abs(shifted[2][1] - (w/2+1)) < pixthr + assert abs(shifted[0][1] - h/2) < pixthr + assert abs(shifted[0][0] - (w/2+1)) < pixthr + assert abs(shifted[1][1] - h/2) < pixthr + assert abs(shifted[1][0] - (w/2+2)) < pixthr + assert abs(shifted[2][1] - (h/2+1)) < pixthr + assert abs(shifted[2][0] - (w/2+1)) < pixthr def test_shift_reverse(): # star1 on frame1 @@ -223,25 +217,24 @@ def test_shift_reverse(): s2y2 = h/2 s2x2 = w/2+2 - s1pos1 = (s1y1, s1x1) - s1pos2 = (s1y2, s1x2) + s1pos1 = (s1x1, s1y1) + s1pos2 = (s1x2, s1y2) - s2pos1 = (s2y1, s2x1) - s2pos2 = (s2y2, s2x2) + s2pos1 = (s2x1, s2y1) + s2pos2 = (s2x2, s2y2) - positions = [(h/2, w/2), (h/2, w/2+1), (h/2+1, w/2)] + positions = [(w/2, h/2), (w/2+1, h/2), (w/2, h/2+1)] - proj = perspective.Projection(W, H, F, w, h) movement = flat.Movement.build(s1pos1, s2pos1, s1pos2, s2pos2) - shifted = movement.reverse(positions, proj) + shifted = movement.reverse(positions) assert len(shifted) == 3 - assert abs(shifted[0][0] - h/2) < pixthr - assert abs(shifted[0][1] - (w/2-1)) < pixthr - assert abs(shifted[1][0] - h/2) < pixthr - assert abs(shifted[1][1] - (w/2)) < pixthr - assert abs(shifted[2][0] - (h/2+1)) < pixthr - assert abs(shifted[2][1] - (w/2-1)) < pixthr + assert abs(shifted[0][1] - h/2) < pixthr + assert abs(shifted[0][0] - (w/2-1)) < pixthr + assert abs(shifted[1][1] - h/2) < pixthr + assert abs(shifted[1][0] - (w/2)) < pixthr + assert abs(shifted[2][1] - (h/2+1)) < pixthr + assert abs(shifted[2][0] - (w/2-1)) < pixthr def test_rotate_shift_forward(): # star1 on frame1 @@ -260,23 +253,22 @@ def test_rotate_shift_forward(): s2y2 = h/2+1 s2x2 = w/2+1 - s1pos1 = (s1y1, s1x1) - s1pos2 = (s1y2, s1x2) + s1pos1 = (s1x1, s1y1) + s1pos2 = (s1x2, s1y2) - s2pos1 = (s2y1, s2x1) - s2pos2 = (s2y2, s2x2) + s2pos1 = (s2x1, s2y1) + s2pos2 = (s2x2, s2y2) - positions = [(h/2, w/2), (h/2, w/2+1)] + positions = [(w/2, h/2), (w/2+1, h/2)] - proj = perspective.Projection(W, H, F, w, h) movement = flat.Movement.build(s1pos1, s2pos1, s1pos2, s2pos2) - shifted = movement.apply(positions, proj) + shifted = movement.apply(positions) assert len(shifted) == 2 - assert abs(shifted[0][0] - h/2) < pixthr - assert abs(shifted[0][1] - (w/2+1)) < pixthr - assert abs(shifted[1][0] - (h/2+1)) < pixthr - assert abs(shifted[1][1] - (w/2+1)) < pixthr + assert abs(shifted[0][1] - h/2) < pixthr + assert abs(shifted[0][0] - (w/2+1)) < pixthr + assert abs(shifted[1][1] - (h/2+1)) < pixthr + assert abs(shifted[1][0] - (w/2+1)) < pixthr def test_multiply_1(): # star1 on frame1 @@ -303,22 +295,21 @@ def test_multiply_1(): s2y3 = h/2+1 s2x3 = w/2+3 - s1pos1 = (s1y1, s1x1) - s1pos2 = (s1y2, s1x2) - s1pos3 = (s1y3, s1x3) + s1pos1 = (s1x1, s1y1) + s1pos2 = (s1x2, s1y2) + s1pos3 = (s1x3, s1y3) - s2pos1 = (s2y1, s2x1) - s2pos2 = (s2y2, s2x2) - s2pos3 = (s2y3, s2x3) + s2pos1 = (s2x1, s2y1) + s2pos2 = (s2x2, s2y2) + s2pos3 = (s2x3, s2y3) movement1 = flat.Movement.build(s1pos1, s2pos1, s1pos2, s2pos2) movement2 = flat.Movement.build(s1pos2, s2pos2, s1pos3, s2pos3) movement = movement2 * movement1 - proj = perspective.Projection(W, H, F, w, h) positions = [s1pos1, s2pos1] - shifted = movement.apply(positions, proj) + shifted = movement.apply(positions) assert len(shifted) == 2 compare_points(shifted[0], s1pos3) @@ -349,22 +340,21 @@ def test_multiply_2(): s2y3 = h/2-1 s2x3 = w/2 - s1pos1 = (s1y1, s1x1) - s1pos2 = (s1y2, s1x2) - s1pos3 = (s1y3, s1x3) + s1pos1 = (s1x1, s1y1) + s1pos2 = (s1x2, s1y2) + s1pos3 = (s1x3, s1y3) - s2pos1 = (s2y1, s2x1) - s2pos2 = (s2y2, s2x2) - s2pos3 = (s2y3, s2x3) + s2pos1 = (s2x1, s2y1) + s2pos2 = (s2x2, s2y2) + s2pos3 = (s2x3, s2y3) movement1 = flat.Movement.build(s1pos1, s2pos1, s1pos2, s2pos2) movement2 = flat.Movement.build(s1pos2, s2pos2, s1pos3, s2pos3) movement = movement2 * movement1 - proj = perspective.Projection(W, H, F, w, h) positions = [s1pos1, s2pos1] - shifted = movement.apply(positions, proj) + shifted = movement.apply(positions) assert len(shifted) == 2 compare_points(shifted[0], s1pos3) @@ -395,22 +385,21 @@ def test_multiply_3(): s2y3 = h/2+1 s2x3 = w/2 - s1pos1 = (s1y1, s1x1) - s1pos2 = (s1y2, s1x2) - s1pos3 = (s1y3, s1x3) + s1pos1 = (s1x1, s1y1) + s1pos2 = (s1x2, s1y2) + s1pos3 = (s1x3, s1y3) - s2pos1 = (s2y1, s2x1) - s2pos2 = (s2y2, s2x2) - s2pos3 = (s2y3, s2x3) + s2pos1 = (s2x1, s2y1) + s2pos2 = (s2x2, s2y2) + s2pos3 = (s2x3, s2y3) movement1 = flat.Movement.build(s1pos1, s2pos1, s1pos2, s2pos2) movement2 = flat.Movement.build(s1pos2, s2pos2, s1pos3, s2pos3) movement = movement2 * movement1 - proj = perspective.Projection(W, H, F, w, h) positions = [s1pos1, s2pos1] - shifted = movement.apply(positions, proj) + shifted = movement.apply(positions) assert len(shifted) == 2 compare_points(shifted[0], s1pos3) @@ -433,22 +422,21 @@ def test_multiply_4(): s2y2 = h/2 + 1/2**0.5 s2x2 = w/2 + 1/2**0.5 + 1 - s1pos1 = (s1y1, s1x1) - s1pos2 = (s1y2, s1x2) + s1pos1 = (s1x1, s1y1) + s1pos2 = (s1x2, s1y2) - s2pos1 = (s2y1, s2x1) - s2pos2 = (s2y2, s2x2) + s2pos1 = (s2x1, s2y1) + s2pos2 = (s2x2, s2y2) - movement1 = flat.Movement(0, -h/2, -w/2) + movement1 = flat.Movement(0, -w/2, -h/2) movement2 = flat.Movement(math.pi/4, 0, 0) - movement3 = flat.Movement(0, 0, 1) - movement4 = flat.Movement(0, h/2, w/2) + movement3 = flat.Movement(0, 1, 0) + movement4 = flat.Movement(0, w/2, h/2) movement = movement4 * movement3 * movement2 * movement1 - proj = perspective.Projection(W, H, F, w, h) positions = [s1pos1, s2pos1] - shifted = movement.apply(positions, proj) + shifted = movement.apply(positions) assert len(shifted) == 2 compare_points(shifted[0], s1pos2) @@ -471,22 +459,21 @@ def test_inverse(): s2y2 = h/2+1 s2x2 = w/2+1 - s1pos1 = (s1y1, s1x1) - s1pos2 = (s1y2, s1x2) + s1pos1 = (s1x1, s1y1) + s1pos2 = (s1x2, s1y2) - s2pos1 = (s2y1, s2x1) - s2pos2 = (s2y2, s2x2) + s2pos1 = (s2x1, s2y1) + s2pos2 = (s2x2, s2y2) positions = [s1pos1, s2pos1] - proj = perspective.Projection(W, H, F, w, h) movement = flat.Movement.build(s1pos1, s2pos1, s1pos2, s2pos2) rev_movement = movement.inverse() mov1 = movement * rev_movement mov2 = rev_movement * movement - shifted1 = mov1.apply(positions, proj) - shifted2 = mov2.apply(positions, proj) + shifted1 = mov1.apply(positions) + shifted2 = mov2.apply(positions) assert len(shifted1) == 2 assert len(shifted2) == 2 compare_points(shifted1[0], s1pos1) diff --git a/tests/test_movement_sphere.py b/tests/test_movement_sphere.py index 28a5daee..9a6b50cf 100644 --- a/tests/test_movement_sphere.py +++ b/tests/test_movement_sphere.py @@ -14,8 +14,9 @@ import sys import math +import numpy as np -import vstarstack.library.projection.perspective as perspective +from vstarstack.library.projection.projections import PerspectiveProjection from vstarstack.library.movement import sphere sys.path.append('../') @@ -49,22 +50,21 @@ def test_no_rotation_forward(): s2lat2 = 0 s2lon2 = 1 - s1pos1 = (s1lat1, s1lon1) - s1pos2 = (s1lat2, s1lon2) + s1pos1 = (s1lon1, s1lat1) + s1pos2 = (s1lon2, s1lat2) - s2pos1 = (s2lat1, s2lon1) - s2pos2 = (s2lat2, s2lon2) + s2pos1 = (s2lon1, s2lat1) + s2pos2 = (s2lon2, s2lat2) - positions = [(h/2, w/2)] + positions = np.array([[w/2, h/2]], dtype='double') - proj = perspective.Projection(W, H, F, w, h) + proj = PerspectiveProjection(w, h, W, H, F) movement = sphere.Movement.build(s1pos1, s2pos1, s1pos2, s2pos2) - shifted = movement.apply(positions, proj) + shifted = movement.apply(positions, proj, proj) assert len(shifted) == 1 - assert abs(shifted[0][0] - h/2) < pixthr - assert abs(shifted[0][1] - w/2) < pixthr - + assert abs(shifted[0][0] - w/2) < pixthr + assert abs(shifted[0][1] - h/2) < pixthr def test_no_rotation_reverse(): s1lat1 = 0 @@ -79,21 +79,21 @@ def test_no_rotation_reverse(): s2lat2 = 0 s2lon2 = 1 - s1pos1 = (s1lat1, s1lon1) - s1pos2 = (s1lat2, s1lon2) + s1pos1 = (s1lon1, s1lat1) + s1pos2 = (s1lon2, s1lat2) - s2pos1 = (s2lat1, s2lon1) - s2pos2 = (s2lat2, s2lon2) + s2pos1 = (s2lon1, s2lat1) + s2pos2 = (s2lon2, s2lat2) - positions = [(h/2, w/2)] + positions = np.array([[w/2, h/2]], dtype='double') - proj = perspective.Projection(W, H, F, w, h) + proj = PerspectiveProjection(w, h, W, H, F) movement = sphere.Movement.build(s1pos1, s2pos1, s1pos2, s2pos2) - shifted = movement.reverse(positions, proj) + shifted = movement.reverse(positions, proj, proj) assert len(shifted) == 1 - assert abs(shifted[0][0] - h/2) < pixthr - assert abs(shifted[0][1] - w/2) < pixthr + assert abs(shifted[0][0] - w/2) < pixthr + assert abs(shifted[0][1] - h/2) < pixthr def test_lon_rotation_forward(): @@ -111,27 +111,27 @@ def test_lon_rotation_forward(): s2lat2 = 0 s2lon2 = 1 + dlon - s1pos1 = (s1lat1, s1lon1) - s1pos2 = (s1lat2, s1lon2) + s1pos1 = (s1lon1, s1lat1) + s1pos2 = (s1lon2, s1lat2) - s2pos1 = (s2lat1, s2lon1) - s2pos2 = (s2lat2, s2lon2) + s2pos1 = (s2lon1, s2lat1) + s2pos2 = (s2lon2, s2lat2) x = w/2 y = h/2 x_moved_expected = x - w * F/W*math.tan(dlon) y_moved_expected = h/2 - positions = [(y, x)] + positions = np.array([[x, y]], dtype='double') - proj = perspective.Projection(W, H, F, w, h) + proj = PerspectiveProjection(w, h, W, H, F) movement = sphere.Movement.build(s1pos1, s2pos1, s1pos2, s2pos2) - shifted = movement.apply(positions, proj) - assert len(shifted) == 1 - assert abs(shifted[0][0] - y_moved_expected) < pixthr - assert abs(shifted[0][1] - x_moved_expected) < pixthr + shifted = movement.apply(positions, proj, proj) + assert len(shifted) == 1 + assert abs(shifted[0][0] - x_moved_expected) < pixthr + assert abs(shifted[0][1] - y_moved_expected) < pixthr def test_lon_rotation_reverse(): dlon = 0.1 @@ -148,26 +148,26 @@ def test_lon_rotation_reverse(): s2lat2 = 0 s2lon2 = 1 + dlon - s1pos1 = (s1lat1, s1lon1) - s1pos2 = (s1lat2, s1lon2) + s1pos1 = (s1lon1, s1lat1) + s1pos2 = (s1lon2, s1lat2) - s2pos1 = (s2lat1, s2lon1) - s2pos2 = (s2lat2, s2lon2) + s2pos1 = (s2lon1, s2lat1) + s2pos2 = (s2lon2, s2lat2) x = w/2 y = h/2 x_moved_expected = x + w * F/W*math.tan(dlon) y_moved_expected = h/2 - positions = [(y, x)] + positions = np.array([[x, y]], dtype='double') - proj = perspective.Projection(W, H, F, w, h) + proj = PerspectiveProjection(w, h, W, H, F) movement = sphere.Movement.build(s1pos1, s2pos1, s1pos2, s2pos2) - shifted = movement.reverse(positions, proj) + shifted = movement.reverse(positions, proj, proj) assert len(shifted) == 1 - assert abs(shifted[0][0] - y_moved_expected) < pixthr - assert abs(shifted[0][1] - x_moved_expected) < pixthr + assert abs(shifted[0][0] - x_moved_expected) < pixthr + assert abs(shifted[0][1] - y_moved_expected) < pixthr def test_neglon_rotation_forward(): @@ -185,28 +185,30 @@ def test_neglon_rotation_forward(): s2lat2 = 0 s2lon2 = 1 + dlon - s1pos1 = (s1lat1, s1lon1) - s1pos2 = (s1lat2, s1lon2) + s1pos1 = (s1lon1, s1lat1) + s1pos2 = (s1lon2, s1lat2) - s2pos1 = (s2lat1, s2lon1) - s2pos2 = (s2lat2, s2lon2) + s2pos1 = (s2lon1, s2lat1) + s2pos2 = (s2lon2, s2lat2) x = w/2 y = h/2 x_moved_expected = x - w * F/W*math.tan(dlon) y_moved_expected = h/2 - positions = [(y, x)] + positions = np.array([[x, y]], dtype='double') - proj = perspective.Projection(W, H, F, w, h) + proj = PerspectiveProjection(w, h, W, H, F) movement = sphere.Movement.build(s1pos1, s2pos1, s1pos2, s2pos2) - shifted = movement.apply(positions, proj) + shifted = movement.apply(positions, proj, proj) assert len(shifted) == 1 - assert abs(shifted[0][0] - y_moved_expected) < pixthr - assert abs(shifted[0][1] - x_moved_expected) < pixthr + assert abs(shifted[0][0] - x_moved_expected) < pixthr + assert abs(shifted[0][1] - y_moved_expected) < pixthr def test_multiply_1(): + proj = PerspectiveProjection(w, h, W, H, F) + s1lat1 = 0 s1lon1 = 0 @@ -219,11 +221,11 @@ def test_multiply_1(): s2lat2 = 0.01 s2lon2 = 0 - s1pos1 = (s1lat1, s1lon1) - s1pos2 = (s1lat2, s1lon2) + s1pos1 = (s1lon1, s1lat1) + s1pos2 = (s1lon2, s1lat2) - s2pos1 = (s2lat1, s2lon1) - s2pos2 = (s2lat2, s2lon2) + s2pos1 = (s2lon1, s2lat1) + s2pos2 = (s2lon2, s2lat2) movement1 = sphere.Movement.build(s1pos1, s2pos1, s1pos2, s2pos2) s1lat1 = 0 @@ -238,20 +240,19 @@ def test_multiply_1(): s2lat2 = 0 s2lon2 = 1.01 - s1pos1 = (s1lat1, s1lon1) - s1pos2 = (s1lat2, s1lon2) + s1pos1 = (s1lon1, s1lat1) + s1pos2 = (s1lon2, s1lat2) - s2pos1 = (s2lat1, s2lon1) - s2pos2 = (s2lat2, s2lon2) + s2pos1 = (s2lon1, s2lat1) + s2pos2 = (s2lon2, s2lat2) movement2 = sphere.Movement.build(s1pos1, s2pos1, s1pos2, s2pos2) movement = movement2 * movement1 - proj = perspective.Projection(W, H, F, w, h) + positions = np.array([[w/2, h/2]], dtype='double') + expected = proj.reverse(0.01, 0) - positions = [(h/2, w/2)] - expected = proj.reverse(0, 0.01) + shifted = movement.apply(positions, proj, proj) - shifted = movement.apply(positions, proj) assert len(shifted) == 1 compare_points(shifted[0], expected) @@ -268,25 +269,29 @@ def test_inverse(): s2lat2 = 0.01 s2lon2 = 0 - s1pos1 = (s1lat1, s1lon1) - s1pos2 = (s1lat2, s1lon2) + s1pos1 = (s1lon1, s1lat1) + s1pos2 = (s1lon2, s1lat2) + + s2pos1 = (s2lon1, s2lat1) + s2pos2 = (s2lon2, s2lat2) + + proj = PerspectiveProjection(w, h, W, H, F) - s2pos1 = (s2lat1, s2lon1) - s2pos2 = (s2lat2, s2lon2) + pnt1 = proj.reverse(s1lon1, s1lat1) + pnt2 = proj.reverse(s2lon1, s2lat1) - positions = [s1pos1, s2pos1] + positions = np.array([pnt1, pnt2], dtype='double') - proj = perspective.Projection(W, H, F, w, h) movement = sphere.Movement.build(s1pos1, s2pos1, s1pos2, s2pos2) rev_movement = movement.inverse() mov1 = movement * rev_movement mov2 = rev_movement * movement - shifted1 = mov1.apply(positions, proj) - shifted2 = mov2.apply(positions, proj) + shifted1 = mov1.apply(positions, proj, proj) + shifted2 = mov2.apply(positions, proj, proj) assert len(shifted1) == 2 assert len(shifted2) == 2 - compare_points(shifted1[0], s1pos1) - compare_points(shifted1[1], s2pos1) - compare_points(shifted2[0], s1pos1) - compare_points(shifted2[1], s2pos1) + compare_points(shifted1[0], pnt1) + compare_points(shifted1[1], pnt2) + compare_points(shifted2[0], pnt1) + compare_points(shifted2[1], pnt2) diff --git a/tests/test_projection_orthographic.py b/tests/test_projection_orthographic.py index 83944958..269d40fd 100644 --- a/tests/test_projection_orthographic.py +++ b/tests/test_projection_orthographic.py @@ -13,7 +13,7 @@ # import math -from vstarstack.library.projection.orthographic import Projection +from vstarstack.library.projection.projections import OrthographicProjection def equal(val1, val2): """Test that values are same""" @@ -21,36 +21,36 @@ def equal(val1, val2): def test_1(): # sphere - proj = Projection(100, 100, 50, 50, 0, 0) - lat, lon = proj.project(50, 50) - print(lat, lon) + proj = OrthographicProjection(100, 100, 50, 50, 0, 0) + lon, lat = proj.project(50, 50) + print(lon, lat) assert equal(lon, 0) assert equal(lat, 0) def test_2(): - proj = Projection(100, 100, 50, 50, 0, 0) - lat, lon = proj.project(50, 75) + proj = OrthographicProjection(100, 100, 50, 50, 0, 0) + lon, lat = proj.project(75, 50) print(lat, lon) assert equal(lon, math.pi/2) assert equal(lat, 0) def test_3(): - proj = Projection(100, 100, 50, 50, 0, 0) - lat, lon = proj.project(50, 25) - print(lat, lon) + proj = OrthographicProjection(100, 100, 50, 50, 0, 0) + lon, lat = proj.project(25, 50) + print(lon, lat) assert equal(lon, -math.pi/2) assert equal(lat, 0) def test_4(): - proj = Projection(100, 100, 50, 50, 0, 0) - lat, lon = proj.project(75, 50) - print(lat, lon) + proj = OrthographicProjection(100, 100, 50, 50, 0, 0) + lon, lat = proj.project(50, 75) + print(lon, lat) assert equal(lon, 0) assert equal(lat, -math.pi/2) def test_5(): - proj = Projection(100, 100, 50, 50, 0, 0) - lat, lon = proj.project(75, 75) + proj = OrthographicProjection(100, 100, 50, 50, 0, 0) + lon, lat = proj.project(75, 75) print(lat, lon) assert lon is None assert lat is None diff --git a/tests/test_projection_perspective.py b/tests/test_projection_perspective.py index f110b3c1..78319dbd 100644 --- a/tests/test_projection_perspective.py +++ b/tests/test_projection_perspective.py @@ -13,8 +13,7 @@ # import math - -import vstarstack.library.projection.perspective +from vstarstack.library.projection.projections import PerspectiveProjection thr = 1e-6 @@ -32,8 +31,8 @@ def test_center_f(): lon_expected = 0 lat_expected = 0 - proj = vstarstack.library.projection.perspective.Projection(W, H, F, w, h) - lat, lon = proj.project(y, x) + proj = PerspectiveProjection(w, h, W, H, F) + lon, lat = proj.project(x, y) assert abs(lat - lat_expected) < thr assert abs(lon - lon_expected) < thr @@ -51,8 +50,8 @@ def test_center_r(): x_expected = 750 y_expected = 500 - proj = vstarstack.library.projection.perspective.Projection(W, H, F, w, h) - y, x = proj.reverse(lat, lon) + proj = PerspectiveProjection(w, h, W, H, F) + x, y = proj.reverse(lon, lat) assert abs(y - y_expected) < thr assert abs(x - x_expected) < thr @@ -70,8 +69,8 @@ def test_left_f(): lon_expected = math.atan(W/2/F) lat_expected = 0 - proj = vstarstack.library.projection.perspective.Projection(W, H, F, w, h) - lat, lon = proj.project(y, x) + proj = PerspectiveProjection(w, h, W, H, F) + lon, lat = proj.project(x, y) assert abs(lat - lat_expected) < thr assert abs(lon - lon_expected) < thr @@ -89,8 +88,8 @@ def test_left_r(): lon = math.atan(W/2/F) lat = 0 - proj = vstarstack.library.projection.perspective.Projection(W, H, F, w, h) - y, x = proj.reverse(lat, lon) + proj = PerspectiveProjection(w, h, W, H, F) + x, y = proj.reverse(lon, lat) assert abs(y - y_expected) < thr assert abs(x - x_expected) < thr @@ -109,8 +108,8 @@ def test_top_f(): lon_expected = 0 lat_expected = math.atan(H/2/F) - proj = vstarstack.library.projection.perspective.Projection(W, H, F, w, h) - lat, lon = proj.project(y, x) + proj = PerspectiveProjection(w, h, W, H, F) + lon, lat = proj.project(x, y) assert abs(lat - lat_expected) < thr assert abs(lon - lon_expected) < thr @@ -128,8 +127,8 @@ def test_top_r(): lon = 0 lat = math.atan(H/2/F) - proj = vstarstack.library.projection.perspective.Projection(W, H, F, w, h) - y, x = proj.reverse(lat, lon) + proj = PerspectiveProjection(w, h, W, H, F) + x, y = proj.reverse(lon, lat) assert abs(y - y_expected) < thr assert abs(x - x_expected) < thr diff --git a/tests/test_stars_pipeline.py b/tests/test_stars_pipeline.py index fd1e9a45..596d3645 100644 --- a/tests/test_stars_pipeline.py +++ b/tests/test_stars_pipeline.py @@ -17,7 +17,6 @@ import matplotlib.pyplot as plt import vstarstack.library.common -import vstarstack.library.projection.perspective import vstarstack.library.loaders.classic import vstarstack.library.stars.detect import vstarstack.library.stars.describe @@ -28,6 +27,8 @@ import vstarstack.library.movement.move_image import vstarstack.library.merge +from vstarstack.library.projection.projections import PerspectiveProjection + from vstarstack.library.movement.sphere import Movement dir_path = os.path.dirname(os.path.realpath(__file__)) @@ -38,7 +39,7 @@ def test_1(): F = 10000 w = 1270 h = 1270 - proj = vstarstack.library.projection.perspective.Projection(W, H, F, w, h) + proj = PerspectiveProjection(w, h, W, H, F) names = ["1.png", "2.png", "3.png", "4.png"] images = [] diff --git a/tests/test_usage.py b/tests/test_usage.py index e9368c03..6189bcc7 100644 --- a/tests/test_usage.py +++ b/tests/test_usage.py @@ -23,7 +23,7 @@ def test_path_autocompletion_1(): assert variants[0] == ("stars", True) def test_path_autocompletion_2(): - variants = complete_path_in_dir(dir_path, "s") + variants = complete_path_in_dir(dir_path, "st") assert len(variants) == 1 assert variants[0] == ("stars", True) From 811be3962e98e1bb6d2aec6c92180727eb3a709b Mon Sep 17 00:00:00 2001 From: Vladislav Tsendrovskii Date: Mon, 14 Aug 2023 16:57:26 +0500 Subject: [PATCH 2/7] Make sphere Movement pickeable --- src/vstarstack/library/movement/sphere.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/vstarstack/library/movement/sphere.py b/src/vstarstack/library/movement/sphere.py index 31e16d9b..95c5f5c2 100644 --- a/src/vstarstack/library/movement/sphere.py +++ b/src/vstarstack/library/movement/sphere.py @@ -15,6 +15,7 @@ import math import json +from typing import Any import numpy as np from scipy.spatial.transform import Rotation @@ -63,6 +64,9 @@ def __init__(self, rot): q = self.rot.as_quat() self.mov = SphereMovement(q[3], q[0], q[1], q[2]) + def __reduce__(self) -> str | tuple[Any, ...]: + return (self.__class__, (self.rot, )) + def magnitude(self): """Magnitude of movement""" rotvec = self.rot.as_rotvec() From 1d02bbc69c47f0f1daa2e1ef4256b426a30cf0b0 Mon Sep 17 00:00:00 2001 From: Vladislav Tsendrovskii Date: Mon, 14 Aug 2023 18:12:09 +0500 Subject: [PATCH 3/7] Fix memory leak --- src/vstarstack/library/movement/move_image.py | 4 +- .../library/movement/movements/module.c | 39 ++++++++++--------- .../library/projection/projections/module.c | 7 ++-- src/vstarstack/tool/shift.py | 10 ++++- 4 files changed, 35 insertions(+), 25 deletions(-) diff --git a/src/vstarstack/library/movement/move_image.py b/src/vstarstack/library/movement/move_image.py index 2d40d3d8..52912013 100644 --- a/src/vstarstack/library/movement/move_image.py +++ b/src/vstarstack/library/movement/move_image.py @@ -52,7 +52,9 @@ def move_image(image: np.ndarray, image_weight_layer = np.ones(shape)*image_weight positions = _generate_points(h, w) - original_positions = transformation.reverse(positions.astype('double'), input_proj, output_proj) + original_positions = transformation.reverse(positions.astype('double'), + input_proj, + output_proj) num = positions.shape[0] transform_array = np.zeros([h, w, 2], dtype='double') for index in range(num): diff --git a/src/vstarstack/library/movement/movements/module.c b/src/vstarstack/library/movement/movements/module.c index e061395a..d53cc4dd 100644 --- a/src/vstarstack/library/movement/movements/module.c +++ b/src/vstarstack/library/movement/movements/module.c @@ -42,10 +42,13 @@ static bool generic_forward_project(void *self, PyObject *_lat = PyTuple_GetItem(res, 1); PyObject *_lon = PyTuple_GetItem(res, 0); if (_lat == NULL || _lon == NULL) + { + Py_DECREF(res); return false; + } *lat = PyFloat_AsDouble(_lat); *lon = PyFloat_AsDouble(_lon); - //printf("xy %lf:%lf -> lonlat %lf:%lf\n", x, y, *lon, *lat); + Py_DECREF(res); return true; } @@ -58,10 +61,13 @@ static bool generic_reverse_project(void *self, PyObject *_y = PyTuple_GetItem(res, 1); PyObject *_x = PyTuple_GetItem(res, 0); if (_y == NULL || _x == NULL) + { + Py_DECREF(res); return false; + } *y = PyFloat_AsDouble(_y); *x = PyFloat_AsDouble(_x); - //printf("lonlat %lf:%lf -> xy %lf:%lf\n", lon, lat, *x, *y); + Py_DECREF(res); return true; } @@ -87,9 +93,17 @@ static int SphereMovements_init(PyObject *_self, PyObject *args, PyObject *kwds) return 0; } -bool is_perspective(PyObject *obj) +static void SphereMovements_finalize(PyObject *_self) { - return false; + PyObject *error_type, *error_value, *error_traceback; + + /* Save the current exception, if any. */ + PyErr_Fetch(&error_type, &error_value, &error_traceback); + + struct SphereMovementObject *self = (struct SphereMovementObject *)_self; + + /* Restore the saved exception. */ + PyErr_Restore(error_type, error_value, error_traceback); } typedef void (*action_f)(struct SphereMovement *mov, @@ -139,6 +153,7 @@ static PyObject *apply_action(PyObject *_self, } size_t num = dims[0]; + printf("Processing %i points\n", (int)num); const double *posi = PyArray_DATA(points); @@ -158,21 +173,6 @@ static PyObject *apply_action(PyObject *_self, .forward = generic_forward_project, .reverse = generic_reverse_project, }; -/* - if (is_perspective(input_proj)) - { - in_proj.projection = &(((struct PerspectiveProjectionObject *)input_proj)->proj); - in_proj.forward = perspective_projection_project; - in_proj.reverse = perspective_projection_reverse; - } - - if (is_perspective(output_proj)) - { - out_proj.projection = &(((struct PerspectiveProjectionObject *)output_proj)->proj); - out_proj.forward = perspective_projection_project; - out_proj.reverse = perspective_projection_reverse; - } - */ PyArrayObject *output_points = (PyArrayObject *)PyArray_ZEROS(2, dims, NPY_DOUBLE, 0); if (output_points == NULL) @@ -218,6 +218,7 @@ static PyTypeObject _SphereMovement = { .tp_flags = Py_TPFLAGS_DEFAULT, .tp_new = PyType_GenericNew, .tp_init = SphereMovements_init, + .tp_finalize = SphereMovements_finalize, .tp_methods = _SphereMovements_methods, }; diff --git a/src/vstarstack/library/projection/projections/module.c b/src/vstarstack/library/projection/projections/module.c index dd111a89..5eae033d 100644 --- a/src/vstarstack/library/projection/projections/module.c +++ b/src/vstarstack/library/projection/projections/module.c @@ -95,14 +95,13 @@ static PyObject *Projection_reverse(void *proj, double x, y; double lon, lat; static char *kwlist[] = {"lon", "lat", NULL}; - if (!PyArg_ParseTupleAndKeywords(args, kwds, "dd", kwlist, - &lon, &lat)) + + if (!PyArg_ParseTupleAndKeywords(args, kwds, "dd", kwlist, &lon, &lat)) return Py_BuildValue("(OO)", Py_None, Py_None); - //printf("lat %lf lon %lf\n", lat, lon); if (!function(proj, lat, lon, &y, &x)) return Py_BuildValue("(OO)", Py_None, Py_None); - //printf("x %lf y %lf\n", x, y); + return Py_BuildValue("(dd)", x, y); } diff --git a/src/vstarstack/tool/shift.py b/src/vstarstack/tool/shift.py index 91bf40fd..e2ff5ecb 100644 --- a/src/vstarstack/tool/shift.py +++ b/src/vstarstack/tool/shift.py @@ -29,6 +29,9 @@ import vstarstack.tool.common +ncpu = vstarstack.tool.cfg.nthreads +#ncpu = 1 + def select_shift(project: vstarstack.tool.cfg.Project, argv: list[str]): """Select optimal shift source""" if len(argv) >= 2: @@ -52,9 +55,12 @@ def select_shift(project: vstarstack.tool.cfg.Project, argv: list[str]): def _make_shift(name : str, filename : str, shift : Movement, outfname : str): print(f"Processing {name}") dataframe = vstarstack.library.data.DataFrame.load(filename) + print("Loaded") result = vstarstack.library.movement.move_image.move_dataframe(dataframe, shift) + print("Transformed") vstarstack.tool.common.check_dir_exists(outfname) result.store(outfname) + print("Saved") def apply_shift(project: vstarstack.tool.cfg.Project, argv: list[str]): """Apply shifts to images""" @@ -76,8 +82,10 @@ def apply_shift(project: vstarstack.tool.cfg.Project, argv: list[str]): images = vstarstack.tool.common.listfiles(npy_dir, ".zip") args = [(name, filename, shifts[name], os.path.join(shifted_dir, name + ".zip")) for name, filename in images] - with mp.Pool(vstarstack.tool.cfg.nthreads) as pool: + with mp.Pool(ncpu) as pool: pool.starmap(_make_shift, args) + #for arg in args: + # _make_shift(*arg) commands = { "select-shift": (select_shift, From 702d0133c979e6f7abfcbe8c0ae715297257c104 Mon Sep 17 00:00:00 2001 From: Vladislav Tsendrovskii Date: Mon, 14 Aug 2023 18:34:30 +0500 Subject: [PATCH 4/7] Debug messaging --- src/vstarstack/library/movement/movements/module.c | 2 +- src/vstarstack/tool/shift.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/vstarstack/library/movement/movements/module.c b/src/vstarstack/library/movement/movements/module.c index d53cc4dd..f5994dd6 100644 --- a/src/vstarstack/library/movement/movements/module.c +++ b/src/vstarstack/library/movement/movements/module.c @@ -153,7 +153,7 @@ static PyObject *apply_action(PyObject *_self, } size_t num = dims[0]; - printf("Processing %i points\n", (int)num); + //printf("Processing %i points\n", (int)num); const double *posi = PyArray_DATA(points); diff --git a/src/vstarstack/tool/shift.py b/src/vstarstack/tool/shift.py index e2ff5ecb..3c17fbf6 100644 --- a/src/vstarstack/tool/shift.py +++ b/src/vstarstack/tool/shift.py @@ -55,12 +55,12 @@ def select_shift(project: vstarstack.tool.cfg.Project, argv: list[str]): def _make_shift(name : str, filename : str, shift : Movement, outfname : str): print(f"Processing {name}") dataframe = vstarstack.library.data.DataFrame.load(filename) - print("Loaded") + print(f"Loaded {name}") result = vstarstack.library.movement.move_image.move_dataframe(dataframe, shift) - print("Transformed") + print(f"Transformed {name}") vstarstack.tool.common.check_dir_exists(outfname) result.store(outfname) - print("Saved") + print(f"Saved {name}") def apply_shift(project: vstarstack.tool.cfg.Project, argv: list[str]): """Apply shifts to images""" From b373d4265c875c958256147f1bebec1907635e73 Mon Sep 17 00:00:00 2001 From: Vladislav Tsendrovskii Date: Mon, 14 Aug 2023 19:12:18 +0500 Subject: [PATCH 5/7] Fix projection --- src/vstarstack/library/stars/describe.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/vstarstack/library/stars/describe.py b/src/vstarstack/library/stars/describe.py index 126835e2..cc9234bd 100644 --- a/src/vstarstack/library/stars/describe.py +++ b/src/vstarstack/library/stars/describe.py @@ -18,7 +18,7 @@ def project_star(star : dict, proj): """Project star to sphere""" - lat, lon = proj.project(star["y"], star["x"]) + lon, lat = proj.project(star["x"], star["y"]) star["lon"] = lon star["lat"] = lat return star From 3efa7b89cada61db7314b15dde408d0abf9530c4 Mon Sep 17 00:00:00 2001 From: Vladislav Tsendrovskii Date: Mon, 14 Aug 2023 20:49:31 +0500 Subject: [PATCH 6/7] Add test --- tests/test_projection_perspective.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/tests/test_projection_perspective.py b/tests/test_projection_perspective.py index 78319dbd..b889aeab 100644 --- a/tests/test_projection_perspective.py +++ b/tests/test_projection_perspective.py @@ -132,3 +132,17 @@ def test_top_r(): assert abs(y - y_expected) < thr assert abs(x - x_expected) < thr + +def test_custom_1(): + w = 4640 + h = 3506 + F = 420 + W = w*0.0038 + H = h*0.0038 + + x = 1725 + y = 2224 + proj = PerspectiveProjection(w, h, W, H, F) + lon, lat = proj.project(x, y) + + assert abs(lon - 0.0053832813307391) < thr From 503326e0043bf2dc898da77127938a3aeefe767f Mon Sep 17 00:00:00 2001 From: Vladislav Tsendrovskii Date: Mon, 14 Aug 2023 20:53:48 +0500 Subject: [PATCH 7/7] Fix --- src/vstarstack/tool/stars/detect.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/vstarstack/tool/stars/detect.py b/src/vstarstack/tool/stars/detect.py index 0a87d699..830a0716 100644 --- a/src/vstarstack/tool/stars/detect.py +++ b/src/vstarstack/tool/stars/detect.py @@ -28,7 +28,7 @@ def detect_stars(projection, """Detect stars on image""" stars = detect.detect_stars(gray) for star in stars: - star["lat"], star["lon"] = projection.project(star["y"], star["x"]) + star["lon"], star["lat"] = projection.project(star["x"], star["y"]) return sorted(stars, key=lambda x: x["size"], reverse=True) def _process_file(fname, jsonfile):