Непосредственно «отрезки» сегментов линии до массива numpy

Question

Один из моих первых проектов, реализованных в python, делает моделирование перколяции в Монте-Карло. Код постоянно рос. Первой частью была визуализация перколяции палочки. В области ширины * длина заданной плотности (палочки/площади) прямых палок с определенной длиной нанесены со случайными начальными координатами и направлением. Поскольку я часто использую gnuplot, я написал сгенерированные (x, y) начальные и конечные координаты в текстовый файл, чтобы потом gnuplot их.

Я нашел here хороший способ анализа данных изображения с использованием scipy.ndimage.measurements. Изображение читается с помощью ndimage.imread в оттенках серого. Получившийся массив numpy далее сводится к логическим значениям, поскольку меня интересуют только соединения между разными палочками. Затем полученные кластеры анализируются с помощью ndimage.measurements. Это позволяет мне узнать, есть ли пути, которые соединяются с одной стороны на другую или нет. Здесь приведен минимальный пример.

import random 
import math 
from scipy.ndimage import measurements 
from scipy.ndimage import imread 
import numpy as np 
import matplotlib.pyplot as plt 

#dimensions of plot 
width = 10 
length = 8 
stick_length = 1 
fig = plt.figure(frameon=False) 
ax = fig.add_axes([0, 0, 1, 1]) 
fig.set_figwidth(width) 
fig.set_figheight(length) 
ax.axis('off') 

file = open("coordinates.txt", "w") 

for i in range (300): 
    # randomly create (x,y) start coordinates in channel and direction 
    xstart = width * random.random()  # xstart = 18 
    ystart = length * random.random()  # ystart = 2 
    # randomly generate direction of stick from start coordinates and convert from GRAD in RAD 
    dirgrad = 360 * random.random() 
    dirrad = math.radians(dirgrad) 
    # calculate (x,y) end coordinates 
    xend = xstart + (math.cos(dirrad) * stick_length) 
    yend = ystart + (math.sin(dirrad) * stick_length) 
    # write start and end coordinates into text file for gnuplot plotting 
    file.write(str(i) + ":\t" + str(xstart) + "\t" + str(ystart) + "\t" + str(dirgrad) + ":\t" + str(xend) + "\t" + str(yend) + "\n") 
    file.write(str(i) + ":\t" + str(xend) + "\t" + str(yend) + "\n\n") 
    # or plot directly with matplotlib 
    ax.plot([xstart,xend],[ystart,yend],"black", lw=1) 
fig.savefig("testimage.png", dpi=100) 

# now read just saved image and do analysis with scipy.ndimage 
fig1, ax1 = plt.subplots(1,1) 
img_input = imread("testimage.png", flatten = True) # read image to np.ndarray in grey scales 
img_bw = img_input < 255 # convert grey scales to b/w (boolean) 
labeled_array, num_clusters = measurements.label(img_bw) #labeled_array: labeled clusters in array, num_clusters: number of clusters 
area = measurements.sum(img_bw, labeled_array, index=np.arange(labeled_array.max() + 1)) # area of each cluster 
areaImg = area[labeled_array] # label each cluster with labelnumber=area 
cax = ax1.imshow(areaImg, origin='upper', interpolation='nearest', cmap = 'rainbow') 
cbar = fig1.colorbar(cax) 
fig1.savefig("testimage_analyzed.png") 

Хотя это работает в основном только тонкие моделирования методом Монте-Карло с 1000 итераций для большего числа различных плотностей рукояти в конечном итоге работает 8 часов или больше. Отчасти это связано с тем, что созданные изображения & массивов довольно велики, и тысячи палочек построены для более высоких плотностей. Причина в том, что я хочу моделировать диапазон геометрий (например, длина между 500 и 20000 пикселей), одновременно сводя к минимуму ошибку из-за пикселизации.

Я думаю, лучший способ - не использовать данные изображения и рассматривать его как векторную проблему, но я понятия не имею, как начать алгоритм. И многие соединения могут привести к большим массивам данных.

Пребывание с описанным выше способом очевидно, что запись данных в файл и повторное чтение его не очень эффективны. Поэтому я ищу способы ускорить это. В качестве первого шага я использовал matplotlib для создания изображения, но, по крайней мере, при построении каждой палочки с помощью отдельного сюжетного вызова это до 10 раз медленнее для большего количества палочек. Создание списка палочных координат в массиве и построение полного списка с помощью одного сюжетного вызова может ускорить его, но все же оставляет узкое место в написании и чтении изображения.

Можете ли вы указать мне на эффективный метод прямого генерации массива numping типа boolean, представляющего черно-белое изображение палочек? Может ли построить список координат и каким-то образом преобразовать фигуру в массив? Я также нашел это интересное discussion, где линии нанесены на изображение PIL. Может ли это быть быстрее, чем matplotlib?

Answer 1

Рисование сегмента линии в массиве является фундаментальной возможностью любой графической библиотеки. Самый простой способ - это, вероятно, Bresenham's algorithm. Алгоритм прост и быстр - реализован на быстром языке. Я бы не рекомендовал его реализовать в чистом питоне. Недостатком простейшей версии алгоритма является то, что он не сглаживается. Линии показывают "jaggies". Найдите «алгоритмы рисования линий» для более продвинутых методов с лучшим сглаживанием.

У меня Cython implementation of Bresenham's algorithm в моем eyediagram package. Функция bres_segment_count увеличивает значения во входном массиве по прямой от (x0, y0) до (x1, y1). Модификация, которая просто устанавливает значения массива в 1, будет тривиальным изменением этого кода.

Например,

In [21]: dim = 250 

In [22]: num_sticks = 300 

Каждая строка sticks содержит [x0, y0, x1, y1], конечные точки на "палке":

In [23]: sticks = np.random.randint(0, dim, size=(num_sticks, 4)).astype(np.int32) 

In [24]: img = np.zeros((dim, dim), dtype=np.int32) 

bres_segments_count рисует каждую палочку, используя Алгоритм Брешенема. Обратите внимание, что вместо простой установки значения в строке, например, 1, значения в img по линии увеличиваются.

In [25]: from eyediagram._brescount import bres_segments_count 

In [26]: bres_segments_count(sticks, img) 

In [27]: plt.imshow(img, interpolation='nearest', cmap=cm.hot) 
Out[27]: <matplotlib.image.AxesImage at 0x10f94b110> 

Вот сюжет, который генерируется: