Rotate sky with module on C

.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,

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',

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:

src/vstarstack/library/movement/basic_movement.py

@@ -13,6 +13,7 @@
 # along with this program. If not, see <https://www.gnu.org/licenses/>.
 #
 
+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

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:

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)

src/vstarstack/library/movement/move_image.py

@@ -14,29 +14,31 @@
 #
 
 import numpy as np
+import scipy.ndimage
+
 import vstarstack.library.common
 import vstarstack.library.data
 from vstarstack.library.data import DataFrame
 
 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,28 @@ 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 +93,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")

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 <https://www.gnu.org/licenses/>.
+ */
+
+#include <math.h>
+#include "flat.h"
+

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 <https://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include <unistd.h>