|
| 1 | +import numpy as np |
| 2 | +from timeit import default_timer as timer |
| 3 | + |
| 4 | + |
| 5 | +def make_circular_index(pixel_radius, sample_factor=10): |
| 6 | + """ Creates a circular index around the point at index 0, 0 |
| 7 | + This method is sample-based and therefore not precise. |
| 8 | + Solution is found by sampling at regular interval along the circle |
| 9 | + and storing the resulting points in 2d grid of resolution 1. |
| 10 | + Points are returned in order, forming a 2d grid circular path.""" |
| 11 | + px_r = pixel_radius # radius of the circle in pixels |
| 12 | + |
| 13 | + # ensures sample_factor samples per angular pixel. |
| 14 | + kNSamples = max(12, 2 * np.pi * px_r * sample_factor) |
| 15 | + |
| 16 | + theta = np.arange(0, kNSamples) * 2 * np.pi / kNSamples |
| 17 | + x = px_r * np.cos(theta) |
| 18 | + y = px_r * np.sin(theta) |
| 19 | + i = np.round(x).astype(int) |
| 20 | + j = np.round(y).astype(int) |
| 21 | + ij = np.vstack([i,j]) |
| 22 | + unique_id = np.sort(np.unique(ij, axis=1, return_index=True)[1]) |
| 23 | + i = ij[0, unique_id] |
| 24 | + j = ij[1, unique_id] |
| 25 | + return i, j |
| 26 | + |
| 27 | +def index_to_grid(i, j, gridsize=None): |
| 28 | + """ Converts the circle from index format to 2d-grid format """ |
| 29 | + size = gridsize |
| 30 | + if gridsize is None: |
| 31 | + size_i = max(i) - min(i) + 3 |
| 32 | + size_j = max(j) - min(j) + 3 |
| 33 | + size = max(size_i, size_j) |
| 34 | + mid = int((size - 1) / 2) |
| 35 | + grid = np.zeros((size, size)) |
| 36 | + grid[(i+mid, j+mid)] = 1 |
| 37 | + return grid |
| 38 | + |
| 39 | +class CircularIndexCreator(object): |
| 40 | + """ Allows creating circular indices for given radii, |
| 41 | + and stores the results in a lookup table to avoid repeat computation """ |
| 42 | + def __init__(self): |
| 43 | + self.px_r_archive = {} |
| 44 | + self.indices_archive = {} |
| 45 | + |
| 46 | + def make_circular_index(self, radius, resolution=1.): |
| 47 | + r = radius # [m] |
| 48 | + res = resolution # [m / px] |
| 49 | + px_r = radius / resolution # pixel-radius of circle [px] |
| 50 | + |
| 51 | + if px_r in self.px_r_archive: |
| 52 | + i, j = self.px_r_archive[px_r] |
| 53 | + return i, j |
| 54 | + |
| 55 | + i, j = make_circular_index(px_r) |
| 56 | + # store result in lookup table |
| 57 | + self.px_r_archive[px_r] = (i, j) |
| 58 | + return i, j |
| 59 | + |
| 60 | +if __name__ == "__main__": |
| 61 | + from matplotlib import pyplot as plt |
| 62 | + N_CIRCLES = 1000 |
| 63 | + PLOT_RADII = False |
| 64 | + PLOT_EVOLUTION = True |
| 65 | + px_radii = np.linspace(0, 100, N_CIRCLES) |
| 66 | + cc = CircularIndexCreator() |
| 67 | + circle_lengths = [] |
| 68 | + t_st = timer() |
| 69 | + for px_r in px_radii: |
| 70 | + i, j = cc.make_circular_index(px_r, 1) |
| 71 | + circle_lengths.append(len(i)) |
| 72 | + t_end = timer() |
| 73 | + print("First-time computation of {} circles took {} seconds.".format(N_CIRCLES, t_end - t_st)) |
| 74 | + t_st = timer() |
| 75 | + for px_r in px_radii: |
| 76 | + i, j = cc.make_circular_index(px_r, 1) |
| 77 | + t_end = timer() |
| 78 | + print("Second-time computation of {} circles took {} seconds.".format(N_CIRCLES, t_end - t_st)) |
| 79 | + if PLOT_RADII: |
| 80 | + plt.figure() |
| 81 | + for px_r in px_radii: |
| 82 | + i, j = cc.make_circular_index(px_r, 1) |
| 83 | + plt.plot(i, j, '-x') |
| 84 | + plt.show() |
| 85 | + plt.figure() |
| 86 | + plt.plot(px_radii, circle_lengths, '-x') |
| 87 | + plt.show() |
| 88 | + if PLOT_EVOLUTION: |
| 89 | + print("Plotting evolution of circles with radius") |
| 90 | + N = 20 |
| 91 | + px_radii = np.logspace(np.log(0.5),np.log10(N-1),200) |
| 92 | + for px_r in px_radii: |
| 93 | + i, j = cc.make_circular_index(px_r, 1) |
| 94 | + plt.ion() |
| 95 | + plt.figure() |
| 96 | + for px_r in px_radii: |
| 97 | + plt.gca().clear() |
| 98 | + i, j = cc.make_circular_index(px_r, 1) |
| 99 | + grid = index_to_grid(i, j, 2*N+1) |
| 100 | + plt.imshow(1-grid, cmap='gray') |
| 101 | + plt.gca().add_patch(plt.Circle((N,N), radius=px_r, color='white', fill=False)) |
| 102 | + plt.pause(0.005) |
| 103 | + plt.ioff() |
0 commit comments