From f85f23776ad368c4a4a0d562aa53f851b06e2025 Mon Sep 17 00:00:00 2001 From: Alexandre Foucher Date: Thu, 3 Oct 2024 15:47:52 +0200 Subject: [PATCH] 3/10/24 --- evaluate/eval.py | 74 ++++++++++++ example/mosse_viz.py | 184 ++++++++++++++--------------- requirements.txt | 4 +- src/cppskyline/include/Terrain.hpp | 4 +- src/cppskyline/src/Terrain.cpp | 31 +++-- src/pyskyline/mosse_correlation.py | 2 +- src/pyskyline/terrain.py | 17 ++- 7 files changed, 208 insertions(+), 108 deletions(-) create mode 100644 evaluate/eval.py diff --git a/evaluate/eval.py b/evaluate/eval.py new file mode 100644 index 0000000..85040c8 --- /dev/null +++ b/evaluate/eval.py @@ -0,0 +1,74 @@ +#!/usr/bin/env python3 + +import math +import cv2 as cv +import pyskyline +import numpy as np +import argparse +import matplotlib.pyplot as plt +from matplotlib.widgets import Slider + +def show_result(terrain, skyline, coord, scoremap): + fig = plt.figure() + gs = fig.add_gridspec(2,3) + fig.suptitle("Mosse correlation", fontsize=16) + ax1 = fig.add_subplot(gs[0, 0]) + ax2 = fig.add_subplot(gs[0, 1]) + ax3 = fig.add_subplot(gs[0, 2]) + ax4 = fig.add_subplot(gs[1, :]) + + im1 = ax1.imshow(scoremap, cmap='hot', interpolation='nearest') + ax1.set_title('Correlation map') + cbar1 = ax1.figure.colorbar(im1) + cbar1.ax.set_ylabel("Correlation score", rotation=-90, va="bottom") + + rgb = cv.cvtColor(terrain.rgb(), cv.COLOR_BGR2RGB) + rgb = cv.circle(rgb, coord, 2, (255,0,0), -1) + ax2.imshow(rgb) + ax2.set_title('Colored terrain') + + im2 = ax3.imshow((terrain.grayscale().astype(float) - 135) * (40.0/120.0)) + ax3.set_title('DEM') + cbar = ax1.figure.colorbar(im2) + cbar.ax.set_ylabel("Altitude", rotation=-90, va="bottom") + + x = np.linspace(-180,180,skyline.shape[0]) + ax4.plot(x,skyline,label="Agent vision") + ax4.plot(x,terrain.get_skyline(coord[0], coord[1]), label="DEM data") + ax4.legend() + + plt.show() + +def eval_score(correlation_map, coord, scale, alpha=10): + alpha = alpha / scale + score = 0 + h,w = correlation_map.shape + x,y = coord + for j in range(h): + for i in range(w): + dist = math.sqrt((x-i)**2+(y-j)**2) + score += correlation_map[j,i] * 1 #(math.pow(0.5,math.pow(dist/alpha,2)-1)-1) + return score + +if __name__ == '__main__': + parser = argparse.ArgumentParser() + + parser.add_argument("--size", help="Terrain size (default=128)", type=int, default=128) + parser.add_argument("--seed", help="Terrain seed (If not set seed is random)", type=int) + + args = parser.parse_args() + + seed = np.random.randint(0,1024) if args.seed is None else args.seed + terrain = pyskyline.Terrain.generate(size=args.size,seed=seed,scale=8.0,max_altitude=40.0) + + terrain.precompute_all_skylines() + + x, y = terrain.get_random_water_coordinate() + skyline = terrain.get_skyline(x, y) + angle = np.random.randint(1,40) + skyline = np.roll(skyline, angle) + correlation_map = pyskyline.MosseCorrelation.process(terrain, skyline) + score = eval_score(correlation_map,(x,y), terrain._scale,100) + score = score / (2 * terrain.water_tile_count()) + 0.5 + print(f"{score:0.2f} with angle of {angle}") + show_result(terrain, skyline, (x,y), correlation_map) diff --git a/example/mosse_viz.py b/example/mosse_viz.py index 63902d4..c281cac 100755 --- a/example/mosse_viz.py +++ b/example/mosse_viz.py @@ -1,11 +1,9 @@ import numpy as np import scipy.signal as signal -from scipy.fftpack import fft, fftshift, ifft +from scipy.fftpack import fft, ifft import matplotlib.pyplot as plt from matplotlib.widgets import Slider -N = 256 - def gaussian_filter1d(size,sigma): filter_range = np.linspace(-int(size/2),int(size/2),size) gaussian_filter = [1 / (sigma * np.sqrt(2*np.pi)) * np.exp(-x**2/(2*sigma**2)) for x in filter_range] @@ -18,117 +16,115 @@ def generate_signal(N:int) -> np.ndarray: y[i] = np.random.normal(scale=1) + (y[i-1] if i > 1 else 0) return np.convolve(y,gaussian_filter1d(N,1),'same') -if __name__ == '__main__': +class UI: + def __init__(self) -> None: + self._size_F = 256 + self._size_H = 256 + self._sigma = 1.1 + self._seed = 0 + self._shift = 0 + self._fig, self._axs = self.init_ui() - np.random.seed(42) - f_full = generate_signal(2048) - x = np.arange(-180,180,360/N) - f = f_full[1024:1024+N] - h = f_full[1024:1024+N] - g = signal.gaussian(N, std=10,sym=True) + self.update() - F = fft(f) - F_ = np.conjugate(F) - G = fft(g) + def init_ui(self) -> None: + fig, (ax1, ax2, ax3, ax4, ax5, ax6, ax7, ax8) = plt.subplots(8, 1, gridspec_kw={'height_ratios':[4,4,4,1,1,1,1,1]}) - K_ = (G*F_)/(F*F_) + ax1.set_title('Terrain') + self.plt_f, = ax1.plot([]) + self.plt_h, = ax1.plot([]) - H = fft(h) - r = ifft(H*K_) + ax2.set_title('MOSSE response signal') + ax2.set_ylim([0, 1.2]) + self.line_r = ax2.axvline() + self.plt_r, = ax2.plot([]) + self.plt_r2, = ax2.plot([]) - # ========================================== + ax3.set_title('Gaussian') + self.plt_g, = ax3.plot([]) - fig, (ax1, ax2, ax3, ax4, ax5, ax6, ax7) = plt.subplots(7, 1, gridspec_kw={'height_ratios':[4,4,4,1,1,1,1]}) + ax1.set_xlim([-180,180]) + ax2.set_xlim([-180,180]) + ax3.set_xlim([-180,180]) - ax1.set_title('Terrain') - plt_f, = ax1.plot(x,f) - plt_h, = ax1.plot(x,h) + self._slider1 = Slider(ax4, 'sigma', 0.1, 5, valinit=self._sigma) + self._slider2 = Slider(ax5, 'shift', -180, 180, valinit=self._shift, valstep=1) + self._slider3 = Slider(ax6, 'seed', 0, 50, valinit=self._seed, valstep=1) + self._slider4 = Slider(ax7, 'Size f', 64, 1024, valinit=self._size_F, valstep=8) + self._slider5 = Slider(ax8, 'Size h', 64, 1024, valinit=self._size_H, valstep=8) - ax2.set_title('MOSSE response signal') - ax2.set_ylim([0, 1.2]) - line_r = ax2.axvline(x=-N//2+np.argmax(abs(r)), color='r') - plt_r, = ax2.plot(x,abs(r)) - plt_r2, = ax2.plot(x,abs(r)) + self._slider1.on_changed(self.update_sigma) + self._slider2.on_changed(self.update_shift) + self._slider3.on_changed(self.update_seed) + self._slider4.on_changed(self.update_size_f) + self._slider5.on_changed(self.update_size_h) - ax3.set_title('Gaussian') - plt_g, = ax3.plot(x,g) + fig.subplots_adjust(hspace=0.75) - ax1.set_xlim([-180,180]) - ax2.set_xlim([-180,180]) - ax3.set_xlim([-180,180]) + return fig, [ax1, ax2, ax3, ax4, ax5, ax6, ax7, ax8] - slider1 = Slider(ax4, 'sigma', 0.3, 10, valinit=0.1) - slider2 = Slider(ax5, 'shift', -N//2 , N//2, valinit=0, valstep=1) - slider3 = Slider(ax6, 'seed', 0 , 50, valinit=0, valstep=1) - slider4 = Slider(ax7, 'N', 128 , 1024, valinit=256, valstep=8) - - sigma = 0.3 - shift = 0 - - def update(): - # K_ = (G*F_)/(F*F_) - window = np.ones((N)) #signal.windows.hamming(N) - H = fft(h*window) - F = fft(f*window) + def update(self): + np.random.seed(self._seed) + g = signal.windows.gaussian(self._size_F, std=self._sigma,sym=True) + f = generate_signal(self._size_F) * signal.windows.hamming(self._size_F) + h = np.zeros((self._size_F)) + s = self._size_F//2-self._size_H//2 + h[s:s+self._size_H] = (f[s:s+self._size_H] + np.random.normal(0,0.5, self._size_H)) * signal.windows.hamming(self._size_H) + h = np.roll(h,int(self._shift/180*self._size_F//2)) + F = fft(f) + G = fft(g) + H = fft(h) R = H*G/F r = ifft(R) s = np.argmax(abs(r)) r2 = np.copy(r) r2[s-5:s+5] = 0 - plt_g.set_data(x,g) - plt_r.set_data(x,abs(r)) - plt_r2.set_data(x,abs(r2)) - plt_h.set_data(x,h) - plt_f.set_data(x,f) - ax1.set_ylim([min(np.min(h),np.min(f))-1, max(np.max(h),np.max(f))+1]) - ax2.set_ylim([0, np.max(r)+0.2]) - line_r.set_xdata(round((np.argmax(abs(r))/N-0.5)*360)) - fig.canvas.draw_idle() + x = np.linspace(-180,180,self._size_F) + self.plt_g.set_data(x, g) + self.plt_r.set_data(x,abs(r)) + self.plt_r2.set_data(x,abs(r2)) + self.plt_f.set_data(x,f) + self.plt_h.set_data(x,h) + self._axs[0].set_ylim([min(np.min(h),np.min(f))-1, max(np.max(h),np.max(f))+1]) + self._axs[1].set_ylim([0, np.max(r)+0.2]) + self._axs[2].set_ylim([0, np.max(g)*1.1]) + self.line_r.set_xdata([round((np.argmax(abs(r))/self._size_F-0.5)*360)]) + self._fig.canvas.draw_idle() - def update_sigma(val): - global g, G, sigma, N - sigma = val - g = signal.gaussian(N, std=sigma,sym=True) - G = fft(g) - update() + def update_sigma(self, val): + self._sigma = val + self.update() + + def update_size_f(self, val): + if val < self._size_H: + self._slider4.eventson = False + self._slider4.set_val(max(val, self._size_H)) + self._slider4.eventson = True + + self._size_F = max(val, self._size_H) + self.update() - def update_shift(val): - global shift, H, h - shift = -val - h = f_full[1024+round(shift*(N/360)):1024+N+round(shift*(N/360))] - noise = np.random.normal(0,0.5, N) - h = h+noise - update() + def update_size_h(self, val): + if val > self._size_F: + self._slider5.eventson = False + self._slider5.set_val(min(val, self._size_F)) + self._slider5.eventson = True + + self._size_H = min(val, self._size_F) + self.update() - def update_seed(val): - global f_full, f, h, F, H, F_, shift - np.random.seed(val) - f_full = generate_signal(2048) - f = f_full[1024:1024+N] - h = f_full[1024+round(shift*(N/360)):1024+N+round(shift*(N/360))] - noise = np.random.normal(0,0.5, N) - h = h+noise - update() + def update_seed(self, val): + self._seed = val + self.update() - def update_n(val): - global g, G, N, f_full, f, h, F, H, F_, shift, x, sigma, slider2 - N = val - x = np.arange(-180,180,360/N) - g = signal.gaussian(N, std=sigma,sym=True) - G = fft(g) - f = f_full[1024:1024+N] - h = f_full[1024+shift:1024+N+shift] - noise = np.random.normal(0,0.5, N) - h = h+noise + def update_shift(self, val): + self._shift = val + self.update() - update() + def show(self) -> None: + plt.show() - - slider1.on_changed(update_sigma) - slider2.on_changed(update_shift) - slider3.on_changed(update_seed) - slider4.on_changed(update_n) - - plt.subplots_adjust(hspace=0.5) - plt.show() \ No newline at end of file +if __name__ == '__main__': + ui = UI() + ui.show() \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index 0c37a0d..048b2c4 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,2 +1,4 @@ noise==1.2.2 -scipy==1.14.1 \ No newline at end of file +scipy==1.14.1 +matplotlib +pyqt5 \ No newline at end of file diff --git a/src/cppskyline/include/Terrain.hpp b/src/cppskyline/include/Terrain.hpp index 232f59d..d25b2fd 100644 --- a/src/cppskyline/include/Terrain.hpp +++ b/src/cppskyline/include/Terrain.hpp @@ -17,7 +17,7 @@ namespace skyline std::uint16_t m_raycount; float m_raydist; float m_raystep; - std::map m_skylines; + std::map m_skylines; float *compute_skyline(std::int32_t x, std::int32_t y); @@ -25,7 +25,7 @@ namespace skyline Terrain(std::size_t size, std::uint8_t *dem, std::uint8_t ocena_height, float max_altitude, float scale, std::uint16_t raycount, float raydist); ~Terrain(); std::size_t get_size(); - float *get_skyline(std::int32_t x, std::int32_t y); + float *get_skyline(std::uint32_t x, std::uint32_t y); bool is_water(std::int32_t x, std::int32_t y); void precompute_all_skylines(); }; diff --git a/src/cppskyline/src/Terrain.cpp b/src/cppskyline/src/Terrain.cpp index daa75b5..f402d3e 100644 --- a/src/cppskyline/src/Terrain.cpp +++ b/src/cppskyline/src/Terrain.cpp @@ -1,6 +1,7 @@ # include "Terrain.hpp" # include # include +# include // ============================================================================ // P R I V A T E @@ -19,7 +20,6 @@ float *skyline::Terrain::compute_skyline(std::int32_t x, std::int32_t y) int dy = -abs(y1-y0), sy = y0m_dem[y * this->m_size + x] <= this->m_ocean_height; } -float *skyline::Terrain::get_skyline(std::int32_t x, std::int32_t y) +float *skyline::Terrain::get_skyline(std::uint32_t x, std::uint32_t y) { - std::uint64_t id = y * this->m_size + x; + std::uint64_t id = (y << 16) | x; // Check if skyline is already computed if (this->is_water(x,y)) @@ -99,20 +99,33 @@ float *skyline::Terrain::get_skyline(std::int32_t x, std::int32_t y) void skyline::Terrain::precompute_all_skylines() { + std::vector sharedList; + #pragma omp parallel for collapse(2) - for( std::uint32_t x=0; xm_size; x++) + for( int x=0; xm_size; x++) { - for( std::uint32_t y=0; ym_size; y++) + for( int y=0; ym_size; y++) { if (this->is_water(x,y)) { - std::uint64_t id = y * this->m_size + x; - if (this->m_skylines.find(id) == this->m_skylines.end()) { - this->m_skylines.insert({id, this->compute_skyline(x, y)}); - } + #pragma omp critical + sharedList.push_back(y << 16 | x); } } } + + #pragma omp parallel for + for( int i=0; i> 16; + std::uint16_t x = (0x0000FFFF & id); + + float *skyline = this->compute_skyline(x, y); + + #pragma omp critical + this->m_skylines.insert({id, skyline}); + } } extern "C" { diff --git a/src/pyskyline/mosse_correlation.py b/src/pyskyline/mosse_correlation.py index 9b50d03..2f385b3 100644 --- a/src/pyskyline/mosse_correlation.py +++ b/src/pyskyline/mosse_correlation.py @@ -54,4 +54,4 @@ class MosseCorrelation: if print_progress: progress_bar.stop() - return score_map \ No newline at end of file + return score_map/np.max(score_map) \ No newline at end of file diff --git a/src/pyskyline/terrain.py b/src/pyskyline/terrain.py index f50e41f..7d4dc5f 100755 --- a/src/pyskyline/terrain.py +++ b/src/pyskyline/terrain.py @@ -1,6 +1,8 @@ import numpy as np import noise import ctypes +import random +from scipy import ndimage from . import DEFAULT_RAYCOUNT, DEFAULT_RAYDIST @@ -73,12 +75,25 @@ class Terrain(): _LIB.Terrain_precompute_all_skylines(self.obj) print(f"{(time.time()-a):.1f}s") + def get_random_water_coordinate(self, coastal_range=100.0) -> tuple: + r = int(round(coastal_range/self._scale)) + valid_location = self._dem <= GROUND_LAYER['ocean'] + valid_location = ndimage.binary_erosion(valid_location, structure=np.ones((r,r))) + while True: + x = random.randint(0, self._size) + y = random.randint(0, self._size) + if valid_location[y,x] == 1: + return (x,y) + def rgb(self) -> np.ndarray: rgb = np.full((self._size, self._size, 3), COLOR_WATER, dtype=np.uint8) rgb[self._dem > GROUND_LAYER['ocean']] = COLOR_SAND rgb[self._dem > GROUND_LAYER['sand']] = COLOR_GRASS rgb[self._dem > GROUND_LAYER['grass']] = COLOR_TREE - return rgb + return np.asarray(rgb) def grayscale(self) -> np.ndarray: return self._dem + + def water_tile_count(self) -> int: + return (self._dem <= GROUND_LAYER['ocean']).sum()