foucher/pyskyline
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

3/10/24

Browse source
This commit is contained in:
Alexandre Foucher 2024-10-03 15:47:52 +02:00
parent bc69ad483b
commit f85f23776a
7 changed files with 208 additions and 108 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

74
evaluate/eval.py Normal file
View file

@ -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)

180
example/mosse_viz.py
View file

@ -1,11 +1,9 @@
import numpy as np import numpy as np
import scipy.signal as signal import scipy.signal as signal
from scipy.fftpack import fft, fftshift, ifft from scipy.fftpack import fft, ifft
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from matplotlib.widgets import Slider from matplotlib.widgets import Slider
N = 256
def gaussian_filter1d(size,sigma): def gaussian_filter1d(size,sigma):
filter_range = np.linspace(-int(size/2),int(size/2),size) 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] 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) y[i] = np.random.normal(scale=1) + (y[i-1] if i > 1 else 0)
return np.convolve(y,gaussian_filter1d(N,1),'same') 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) self.update()
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)
F = fft(f) def init_ui(self) -> None:
F_ = np.conjugate(F) 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]})
G = fft(g)
K_ = (G*F_)/(F*F_) ax1.set_title('Terrain')
self.plt_f, = ax1.plot([])
self.plt_h, = ax1.plot([])
H = fft(h) ax2.set_title('MOSSE response signal')
r = ifft(H*K_) 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') self._slider1 = Slider(ax4, 'sigma', 0.1, 5, valinit=self._sigma)
plt_f, = ax1.plot(x,f) self._slider2 = Slider(ax5, 'shift', -180, 180, valinit=self._shift, valstep=1)
plt_h, = ax1.plot(x,h) 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') self._slider1.on_changed(self.update_sigma)
ax2.set_ylim([0, 1.2]) self._slider2.on_changed(self.update_shift)
line_r = ax2.axvline(x=-N//2+np.argmax(abs(r)), color='r') self._slider3.on_changed(self.update_seed)
plt_r, = ax2.plot(x,abs(r)) self._slider4.on_changed(self.update_size_f)
plt_r2, = ax2.plot(x,abs(r)) self._slider5.on_changed(self.update_size_h)
ax3.set_title('Gaussian') fig.subplots_adjust(hspace=0.75)
plt_g, = ax3.plot(x,g)
ax1.set_xlim([-180,180]) return fig, [ax1, ax2, ax3, ax4, ax5, ax6, ax7, ax8]
ax2.set_xlim([-180,180])
ax3.set_xlim([-180,180])
slider1 = Slider(ax4, 'sigma', 0.3, 10, valinit=0.1) def update(self):
slider2 = Slider(ax5, 'shift', -N//2 , N//2, valinit=0, valstep=1) np.random.seed(self._seed)
slider3 = Slider(ax6, 'seed', 0 , 50, valinit=0, valstep=1) g = signal.windows.gaussian(self._size_F, std=self._sigma,sym=True)
slider4 = Slider(ax7, 'N', 128 , 1024, valinit=256, valstep=8) f = generate_signal(self._size_F) * signal.windows.hamming(self._size_F)
h = np.zeros((self._size_F))
sigma = 0.3 s = self._size_F//2-self._size_H//2
shift = 0 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))
def update(): F = fft(f)
# K_ = (G*F_)/(F*F_) G = fft(g)
window = np.ones((N)) #signal.windows.hamming(N) H = fft(h)
H = fft(h*window)
F = fft(f*window)
R = H*G/F R = H*G/F
r = ifft(R) r = ifft(R)
s = np.argmax(abs(r)) s = np.argmax(abs(r))
r2 = np.copy(r) r2 = np.copy(r)
r2[s-5:s+5] = 0 r2[s-5:s+5] = 0
plt_g.set_data(x,g) x = np.linspace(-180,180,self._size_F)
plt_r.set_data(x,abs(r)) self.plt_g.set_data(x, g)
plt_r2.set_data(x,abs(r2)) self.plt_r.set_data(x,abs(r))
plt_h.set_data(x,h) self.plt_r2.set_data(x,abs(r2))
plt_f.set_data(x,f) self.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]) self.plt_h.set_data(x,h)
ax2.set_ylim([0, np.max(r)+0.2]) self._axs[0].set_ylim([min(np.min(h),np.min(f))-1, max(np.max(h),np.max(f))+1])
line_r.set_xdata(round((np.argmax(abs(r))/N-0.5)*360)) self._axs[1].set_ylim([0, np.max(r)+0.2])
fig.canvas.draw_idle() 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): def update_sigma(self, val):
global g, G, sigma, N self._sigma = val
sigma = val self.update()
g = signal.gaussian(N, std=sigma,sym=True)
G = fft(g)
update()
def update_shift(val): def update_size_f(self, val):
global shift, H, h if val < self._size_H:
shift = -val self._slider4.eventson = False
h = f_full[1024+round(shift*(N/360)):1024+N+round(shift*(N/360))] self._slider4.set_val(max(val, self._size_H))
noise = np.random.normal(0,0.5, N) self._slider4.eventson = True
h = h+noise
update()
def update_seed(val): self._size_F = max(val, self._size_H)
global f_full, f, h, F, H, F_, shift self.update()
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_n(val): def update_size_h(self, val):
global g, G, N, f_full, f, h, F, H, F_, shift, x, sigma, slider2 if val > self._size_F:
N = val self._slider5.eventson = False
x = np.arange(-180,180,360/N) self._slider5.set_val(min(val, self._size_F))
g = signal.gaussian(N, std=sigma,sym=True) self._slider5.eventson = 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
update() self._size_H = min(val, self._size_F)
self.update()
def update_seed(self, val):
self._seed = val
self.update()
slider1.on_changed(update_sigma) def update_shift(self, val):
slider2.on_changed(update_shift) self._shift = val
slider3.on_changed(update_seed) self.update()
slider4.on_changed(update_n)
plt.subplots_adjust(hspace=0.5) def show(self) -> None:
plt.show() plt.show()
if __name__ == '__main__':
ui = UI()
ui.show()

2
requirements.txt
View file

@ -1,2 +1,4 @@
noise==1.2.2 noise==1.2.2
scipy==1.14.1 scipy==1.14.1
matplotlib
pyqt5

4
src/cppskyline/include/Terrain.hpp
View file

@ -17,7 +17,7 @@ namespace skyline
std::uint16_t m_raycount; std::uint16_t m_raycount;
float m_raydist; float m_raydist;
float m_raystep; float m_raystep;
std::map<std::uint64_t, float*> m_skylines; std::map<std::uint32_t, float*> m_skylines;
float *compute_skyline(std::int32_t x, std::int32_t y); 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 size, std::uint8_t *dem, std::uint8_t ocena_height, float max_altitude, float scale, std::uint16_t raycount, float raydist);
~Terrain(); ~Terrain();
std::size_t get_size(); 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); bool is_water(std::int32_t x, std::int32_t y);
void precompute_all_skylines(); void precompute_all_skylines();
}; };

31
src/cppskyline/src/Terrain.cpp
View file

@ -1,6 +1,7 @@
# include "Terrain.hpp" # include "Terrain.hpp"
# include <math.h> # include <math.h>
# include <iostream> # include <iostream>
# include <vector>
// ============================================================================ // ============================================================================
// P R I V A T E // 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 = y0<y1 ? 1 : -1; int dy = -abs(y1-y0), sy = y0<y1 ? 1 : -1;
int err = dx+dy, e2; int err = dx+dy, e2;
std::uint8_t height = 0; std::uint8_t height = 0;
float max_angle = 0; float max_angle = 0;
float dist, angle; float dist, angle;
@ -78,9 +78,9 @@ bool skyline::Terrain::is_water(int32_t x, int32_t y)
return this->m_dem[y * this->m_size + x] <= this->m_ocean_height; return this->m_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 // Check if skyline is already computed
if (this->is_water(x,y)) 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() void skyline::Terrain::precompute_all_skylines()
{ {
std::vector<std::uint32_t> sharedList;
#pragma omp parallel for collapse(2) #pragma omp parallel for collapse(2)
for( std::uint32_t x=0; x<this->m_size; x++) for( int x=0; x<this->m_size; x++)
{ {
for( std::uint32_t y=0; y<this->m_size; y++) for( int y=0; y<this->m_size; y++)
{ {
if (this->is_water(x,y)) if (this->is_water(x,y))
{ {
std::uint64_t id = y * this->m_size + x; #pragma omp critical
if (this->m_skylines.find(id) == this->m_skylines.end()) { sharedList.push_back(y << 16 | x);
this->m_skylines.insert({id, this->compute_skyline(x, y)});
}
} }
} }
} }
#pragma omp parallel for
for( int i=0; i<sharedList.size(); i++)
{
std::uint32_t id = sharedList.at(i);
std::uint16_t y = (0xFFFF0000 & id) >> 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" { extern "C" {

2
src/pyskyline/mosse_correlation.py
View file

@ -54,4 +54,4 @@ class MosseCorrelation:
if print_progress: if print_progress:
progress_bar.stop() progress_bar.stop()
return score_map return score_map/np.max(score_map)

17
src/pyskyline/terrain.py
View file

@ -1,6 +1,8 @@
import numpy as np import numpy as np
import noise import noise
import ctypes import ctypes
import random
from scipy import ndimage
from . import DEFAULT_RAYCOUNT, DEFAULT_RAYDIST from . import DEFAULT_RAYCOUNT, DEFAULT_RAYDIST
@ -73,12 +75,25 @@ class Terrain():
_LIB.Terrain_precompute_all_skylines(self.obj) _LIB.Terrain_precompute_all_skylines(self.obj)
print(f"{(time.time()-a):.1f}s") 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: def rgb(self) -> np.ndarray:
rgb = np.full((self._size, self._size, 3), COLOR_WATER, dtype=np.uint8) 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['ocean']] = COLOR_SAND
rgb[self._dem > GROUND_LAYER['sand']] = COLOR_GRASS rgb[self._dem > GROUND_LAYER['sand']] = COLOR_GRASS
rgb[self._dem > GROUND_LAYER['grass']] = COLOR_TREE rgb[self._dem > GROUND_LAYER['grass']] = COLOR_TREE
return rgb return np.asarray(rgb)
def grayscale(self) -> np.ndarray: def grayscale(self) -> np.ndarray:
return self._dem return self._dem
def water_tile_count(self) -> int:
return (self._dem <= GROUND_LAYER['ocean']).sum()