initial

.gitignore

@@ -0,0 +1,2 @@
+**/*.egg-info
+/dist

pyproject.toml

@@ -0,0 +1,21 @@
+[build-system]
+requires = ["setuptools>=61.0"]
+build-backend = "setuptools.build_meta"
+[project]
+name = "pyskyline"
+version = "0.0.1"
+authors = [
+  { name="Alexandre FOUCHER", email="foucher@univ-ubs.fr" },
+]
+description = "A small example package"
+readme = "README.md"
+requires-python = ">=3.8"
+classifiers = [
+    "Programming Language :: Python :: 3",
+    "License :: OSI Approved :: MIT License",
+    "Operating System :: OS Independent",
+]
+
+# [project.urls]
+# Homepage = "https://github.com/pypa/sampleproject"
+# Issues = "https://github.com/pypa/sampleproject/issues"

requirements.txt

@@ -0,0 +1,4 @@
+perlin-numpy @ git+https://github.com/pvigier/perlin-numpy
+numpy
+scipy
+opencv-python

src/pyskyline/__init__.py

@@ -0,0 +1,18 @@
+"""
+Copyright (C) 2024 University of South Brittany, Lab-STICC UMR 6285 All Rights Reserved.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+"""
+
+from .raycast import bresenham_line
+from .terrain import Terrain

src/pyskyline/localization/mosse_correlation.py

@@ -0,0 +1,26 @@
+import numpy as np
+import math
+from scipy.fft import fft, ifft
+from scipy.signal import gaussian
+
+g = gaussian(256, std=2.25,sym=True)
+G = fft(g)
+
+def compute_score(F : np.ndarray, h : np.ndarray) -> float:
+    """ Compute the correlation score of the two skyline based on a MOSSE correlation filter.
+
+    Args:
+        F (np.ndarray): Real skyline in fourier domain.
+        h (np.ndarray): Simulated skyline in time domain.
+
+    Returns:
+        float: Correlation score, higher is the score higher the correlation is.
+    """
+
+    H = fft(h)
+    r = abs(ifft(H*G/F))
+    shift = np.argmax(r)
+    peak = r[shift]
+    r[shift-5:shift+5] = 0
+    score = peak/math.sqrt(np.sum(np.power(r,2)))
+    return math.log(score)

src/pyskyline/raycast.py

@@ -0,0 +1,57 @@
+#!/usr/bin/env python
+
+"""
+Copyright (C) 2024 University of South Brittany, Lab-STICC UMR 6285 All Rights Reserved.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+"""
+
+import numpy as np
+
+def bresenham_line(x:int, y:int, dist:float, angle:float, limit:int) -> list:
+    """Generate list of points of all grid cell where raycast hit.
+    see https://zingl.github.io/bresenham.html
+
+    Args:
+        x (int): X coordinate.
+        y (int): Y coordinate.
+        dist (float): Ray distance.
+        angle (float): Ray angle.
+        limit (int): coordinate limit.
+
+    Returns:
+        list: List of points [(x,y),...].
+    """
+    points = []
+
+    x1    = int(x + np.cos(angle)* dist)
+    y1    = int(y + np.sin(angle)* dist)
+    dx    =  abs(x1 - x)
+    sx    = 1 if x < x1 else -1
+    dy    = -abs(y1 - y)
+    sy    = 1 if y < y1 else -1
+    error = dx + dy
+
+    while (x != x1 and y != y1) and (0 <= x < limit and 0 <= y < limit):
+        points.append((x,y))
+
+        e2 = 2 * error
+        if e2 >= dy:
+            error += dy
+            x     += sx
+        if e2 <= dx:
+            error += dx
+            y     += sy
+
+
+    return points

src/pyskyline/terrain.py

@@ -0,0 +1,147 @@
+#!/usr/bin/env python
+
+"""
+Copyright (C) 2024 University of South Brittany, Lab-STICC UMR 6285 All Rights Reserved.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+"""
+
+import os
+import math
+import numpy as np
+import cv2 as cv
+import perlin_numpy as perlin
+from .raycast import bresenham_line
+
+print(__file__)
+
+class Terrain:
+    DEFAULT_MASK       = os.path.dirname(__file__) + '/../../resources/mask.jpg'
+    DEFAULT_RAYCOUNT   = 256
+    DEFAULT_RAYDIST    = 1000.0
+    DEFAULT_LAYERS     = {
+        -1 : (169, 166,  97), # Water
+        18 : (175, 214, 238), # Sand
+        36 : ( 34, 139,  34), # Grass
+        82 : ( 20, 100,  20)  # Tree
+    }
+
+    def __init__(self, seed:int, size:int, height:float, scale:float) -> None:
+        self.seed           = seed
+        self.size           = size
+        self.height         = height
+        self.scale          = scale
+        self.elevation_map  = None
+        self.color_map      = None
+        self.skylines       = None
+
+    def generate_elevation_map(self, mask_path:str) -> None:
+        """ Generate grayscale elevation map from 2D perlin noise
+
+        Args:
+            mask_path (str): Terrain grayscale mask path.
+        """
+        np.random.seed(self.seed)
+        perlin_noise = perlin.generate_fractal_noise_2d(
+            shape   = (512,)*2,
+            res     = (2**3,)*2,
+            octaves = 5
+        )
+
+        # Set noise in [0;1]
+        p_min = np.min(perlin_noise)
+        p_max = np.max(perlin_noise)
+        perlin_noise = (perlin_noise - p_min) / (p_max - p_min)
+
+        if self.size != 512:
+            perlin_noise = cv.resize(perlin_noise, (self.size,)*2, interpolation=cv.INTER_CUBIC)
+
+        # Open and rescale mask
+        mask = cv.imread(mask_path, cv.IMREAD_GRAYSCALE)
+        mask = cv.resize(mask, (self.size,)*2, interpolation=cv.INTER_CUBIC)
+
+        elevation_map = perlin_noise * mask
+
+        self.elevation_map = elevation_map.astype(dtype=np.uint8)
+
+    def generate_colored_map(self, layers:'dict|None'=None) -> None:
+        """ Generate RGB colored map based on elevation map.
+
+        Args:
+            layers (dict): Layer dictionary where the key is the height and value the color in BGR.
+        """
+        if layers is None:
+            layers = self.DEFAULT_LAYERS
+
+        h,w       = self.elevation_map.shape
+        color_map = np.full((h,w,3), layers[-1], dtype=np.uint8)
+
+        for height, color in list(layers.items())[1:]:
+            mask            = self.elevation_map > height
+            color_map[mask] = color
+
+        self.color_map = cv.cvtColor(color_map, cv.COLOR_BGR2RGB)
+
+    def compute_skyline(self, x:int, y:int, ray_count:int=256, ray_dist:float=1000.0) -> np.ndarray:
+        """_summary_
+
+        Args:
+            x (int): X coordinate
+            y (int): Y coordinate
+            ray_count (int, optional): Raycast count. Defaults to 256.
+            ray_dist (float, optional): Raycast distance in meter. Defaults to 1000.0.
+
+        Returns:
+            np.ndarray: skyline array.
+        """
+        ray_step  = 2 * np.pi / (ray_count-1)
+        line      = np.zeros((ray_count), dtype=np.float16)
+
+        for i in range(ray_count):
+            max_height  = -10
+            max_v_angle = 0
+            points      = bresenham_line(x, y, ray_dist/self.scale, ray_step*i, self.size)
+            for x0, y0 in points:
+                height = self.elevation_map[y0,x0]
+                if height > max_height:
+                    max_height = height
+                    dist = math.sqrt((x0-x)**2+(y0-y)**2) * self.scale
+                    if dist != 0.0:
+                        v_angle    = np.arctan((height * self.height/255)/dist)
+                        if v_angle > max_v_angle:
+                            max_v_angle = v_angle
+
+            line[i] = np.rad2deg(max_v_angle)
+
+        return line
+
+    @staticmethod
+    def generate(seed:int=0x7B16, size:int=512, height:float=25.0, scale:float=1.0,
+                 mask_path:str=DEFAULT_MASK) -> 'Terrain':
+        """ Generate a procedural terrain using seed and mask, and precompute all horizon lines
+
+        Args:
+            seed (int, optional): Random seed for the terrain generation. Defaults to 0x7B16.
+            size (int, optional): Size of the generated elevation map in pixel. Defaults to 512.
+            height (float, optional): Maximum altitude of the elevation map. Defaults to 40.0.
+            scale (float, optional): Scale of a pixel in meter. Defaults to 1.0.
+            mask_path (str, optional): Terrain mask path. Defaults to DEFAULT_MASK.
+
+        Returns:
+            Terrain: _description_
+        """
+        terrain = Terrain(seed,size,height,scale)
+        terrain.generate_elevation_map(mask_path)
+        terrain.generate_colored_map()
+
+        return terrain

test/test_terrain.py

@@ -0,0 +1,13 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+import pyskyline
+
+terrain = pyskyline.Terrain.generate()
+h = terrain.compute_skyline(330,330,360)
+
+plt.figure()
+plt.imshow(terrain.elevation_map)
+plt.figure()
+plt.plot(np.arange(0,360,360/h.shape[0]), h)
+plt.show()