В этой статье представлена реализация на Python алгоритма распознавания источников освещения на картах окружения (LDR или HDR) при помощи равнопромежуточной проекции (equirectangular projection). Однако после внесения незначительных изменений её также можно использовать с простыми фоновыми изображениями или кубическими картами. Примеры возможного применения алгоритма: программы трассировки лучей, в которых требуется распознавать первичные источники освещения для испускания из них лучей; в растеризованных рендерерах он может применяться для отбрасывания теней, использующих карту окружения; кроме того, алгоритм также можно применять в программах устранения засветов, например в AR.
Алгоритм состоит из следующих этапов:
- Снижение разрешения исходного изображения, например, до 1024.
- Преобразование изображения в яркость (luminance), при необходимости с размытием изображения.
- Применение метода квази-Монте-Карло.
- Преобразование из сферических координат в равнопромежуточные.
- Фильтрация сэмплов на основании яркости соседа.
- Сортировка сэмплов на основании их яркости.
- Фильтрация сэмплов на основании евклидовой метрики.
- Слияние сэмплов при помощи алгоритма Брезенхэма.
- Вычисление позиции кластера освещения на основании его яркости.
Существует множество алгоритмов снижения разрешения изображений. Билинейная фильтрация — самый быстрый или простой в реализации, к тому же он лучше всего подходит в большинстве случаев. Для преобразования яркости и в LDR-, и HDR-изображениях можно использовать стандартную формулу:
lum = img[:, :, 0] * 0.2126 + img[:, :, 1] * 0.7152 + img[:, :, 2] * 0.0722
Дополнительно можно применить к изображению яркости небольшое размытие, например, в 1-2 пикселя для изображения разрешением 1024, для устранения всех высокочастотных деталей (в частности, вызванных снижением разрешения).
Равнопромежуточная проекция
Самая распространённая проекция в картах окружения — это равнопромежуточная проекция3. Мой алгоритм может работать и с другими проекциями, например, с панорамной и с кубическими картами, однако в статье мы рассмотрим только равнопромежуточную проекцию. Для начала нужно нормализовать координаты изображения:
pos[0] = x / width
pos[1] = y / height Затем нам нужно выполнить преобразование из и в декартовы координаты при помощи сферических координат, т.е. θ и φ, где θ = x * 2π, а φ = y * π.
def sphereToEquirectangular(pos):
invAngles = (0.1591, 0.3183)
xy = (math.atan2(pos[1], pos[0]), math.asin(pos[2]))
xy = (xy[0] * invAngles[0], xy[1] * invAngles[1])
return (xy[0] + 0.5, xy[1] + 0.5)
def equirectangularToSphere(pos):
angles = (1 / 0.1591, 1 / 0.3183)
thetaPhi = (pos[0] - 0.5, pos[1] - 0.5)
thetaPhi = (thetaPhi[0] * angles[0], thetaPhi[1] * angles[1])
length = math.cos(thetaPhi[1])
return (math.cos(thetaPhi[0]) * length, math.sin(thetaPhi[0]) * length, math.sin(thetaPhi[1]))Сэмплирование Хаммерсли
Следующим этапом будет применение к сфере метода квази-Монте-Карло, например, сэмплирования Хаммерсли2:
Можно использовать и другие методы сэмплирования, например Холтона4, однако Хаммерсли быстрее и обеспечивает хорошее распределение сэмплов по сфере. Холтон будет хорошим выбором для сэмплов плоскости, если вместо карты окружения используется простое изображение. Обязательным требованием для сэмплирования Хаммерсли является инверсия корней (ряд) ван дер Корпута, подробнее см. по ссылкам2. Вот его быстрая реализация:
def vdcSequence(bits):
bits = (bits << 16) | (bits >> 16)
bits = ((bits & 0x55555555) << 1) | ((bits & 0xAAAAAAAA) >> 1)
bits = ((bits & 0x33333333) << 2) | ((bits & 0xCCCCCCCC) >> 2)
bits = ((bits & 0x0F0F0F0F) << 4) | ((bits & 0xF0F0F0F0) >> 4)
bits = ((bits & 0x00FF00FF) << 8) | ((bits & 0xFF00FF00) >> 8)
return float(bits) * 2.3283064365386963e-10 # / 0x100000000
def hammersleySequence(i, N):
return (float(i) / float(N), vdcSequence(i))Затем мы используем равномерное наложение на сферу:
def sphereSample(u, v):
PI = 3.14159265358979
phi = v * 2.0 * PI
cosTheta = 2.0 * u - 1.0 # map to -1,1
sinTheta = math.sqrt(1.0 - cosTheta * cosTheta);
return (math.cos(phi) * sinTheta, math.sin(phi) * sinTheta, cosTheta)Для сэмплирования Хаммерсли мы применяем фиксированное количество сэмплов, зависящее от разрешения изображения, и преобразуем из сферических координат в декартовы, а потом в равнопромежуточные:
samplesMultiplier = 0.006
samples = int(samplesMultiplier * width * height)
samplesList = []
# apply hammersley sampling
for i in range(0, samples):
xi = hammersleySequence(i, samples)
xyz = sphereSample(xi[0], xi[1]) # to cartesian
imagePos = sphereToEquirectangular(xyz)
luminance = lum[imagePos[0] * width, imagePos[1] * height]Это даст нам хорошее распределение сэмплов, которые будут проверяться на наличие источников освещения:
Фильтрация источников освещения
В первом проходе фильтрации мы игнорируем все сэмплы, не превосходящие порогового значения яркости (для HDR-карт оно может быть выше), а затем сортируем все сэмплы по их яркости:
localSize = int(float(12) * (width / 1024.0)) + 1
samplesList = []
# apply hammersley sampling
for i in range(0, samples):
xi = hammersleySequence(i, samples)
xyz = sphereSample(xi[0], xi[1]) # to cartesian
imagePos = sphereToEquirectangular(xyz)
luminance = lum[imagePos [0] * width, imagePos [1] * height]
sample = Sample(luminance, imagePos , xyz)
luminanceThreshold = 0.8
#do a neighbour search for the maximum luminance
nLum = computeNeighborLuminance(lum, width, height, sample.imagePos, localSize)
if nLum > luminanceThreshold:
samplesList.append(sample)
samplesList = sorted(samplesList, key=lambda obj: obj.luminance, reverse=True)Следующий проход будет выполнять фильтрацию на основании евклидовой метрики и порогового значения расстояния между пикселями (зависящего от разрешения изображения) — это пространственная структура данных, которую можно использовать, чтобы избавиться от сложности O(N2):
euclideanThreshold = int(float(euclideanThresholdPixel) * (width / 2048.0))
# filter based euclidian distance
filteredCount = len(samplesList)
localIndices = np.empty(filteredCount); localIndices.fill(-1)
for i in range(0, filteredCount):
cpos = samplesList[i].pos
if localIndices[i] == -1:
localIndices[i] = i
for j in range(0, filteredCount):
if i != j and localIndices[j] == -1 and distance2d(cpos, samplesList[j].pos) < euclideanThreshold:
localIndices[j] = iПолучившиеся сэмплы проходят этап слияния, чтобы ещё больше снизить количество источников освещения:
Слияние источников освещения
На последнем этапе выполняется слияние сэмплов, принадлежащих к одному кластеру освещения. Для этого мы можем использовать алгоритм Брезенхэма и начать с сэмплов с самой большой яркости, потому что они уже упорядочены. Когда мы находим источник освещения, удовлетворяющий проверке Брезенхэма, то используем его позицию для изменения позиции источника на основании веса прогона:
# apply bresenham check and compute position of the light clusters
lights = []
finalIndices = np.empty(filteredCount); finalIndices.fill(-1)
for i in localIndices:
sample = samplesList[i]
startPos = sample.pos
if finalIndices[i] == -1:
finalIndices[i] = i
light = Light()
light.originalPos = np.array(sample.pos) # position of the local maxima
light.worldPos = np.array(sample.worldPos)
light.pos = np.array(sample.pos)
light.luminance = sample.luminance
for j in localIndices:
if i != j and finalIndices[j] == -1:
endPos = samplesList[j].pos
if bresenhamCheck(lum, width, height, startPos[0], startPos[1], endPos[0], endPos[1]):
finalIndices[j] = i
# compute final position of the light source
sampleWeight = samplesList[j].luminance / sample.luminance
light.pos = light.pos + np.array(endPos) * sampleWeight
light.pos = light.pos / (1.0 + sampleWeight)
imagePos = light.pos * np.array([1.0 / width, 1.0 / height)
light.worldPos = equirectangularToSphere(imagePos)
lights.append(light)Функция Брезенхэма проверяет наличие непрерывной линии, имеющей одинаковую яркость. Если дельта в текущем пикселе превышает определённый порог, то проверка завершается неудачно:
def bresenhamCheck(lum, imageSize, x0, y0, x1, y1):
dX = int(x1 - x0)
stepX = int((dX > 0) - (dX < 0))
dX = abs(dX) << 1
dY = int(y1 - y0)
stepY = int((dY > 0) - (dY < 0))
dY = abs(dY) << 1
luminanceThreshold = 0.15
prevLum = lum[x0][y0]
sumLum = 0.0
c = 0
if (dX >= dY):
# delta may go below zero
delta = int (dY - (dX >> 1))
while (x0 != x1):
# reduce delta, while taking into account the corner case of delta == 0
if ((delta > 0) or (delta == 0 and (stepX > 0))):
delta -= dX
y0 += stepY
delta += dY
x0 += stepX
sumLum = sumLum + min(lum[x0][y0], 1.25)
c = c + 1
if(abs(sumLum / c - prevLum) > luminanceThreshold and (sumLum / c) < 1.0):
return 0
else:
# delta may go below zero
delta = int(dX - (dY >> 1))
while (y0 != y1):
# reduce delta, while taking into account the corner case of delta == 0
if ((delta > 0) or (delta == 0 and (stepY > 0))):
delta -= dY
x0 += stepX
delta += dX
y0 += stepY
sumLum = sumLum + min(lum[x0][y0], 1.25)
c = c + 1
if(abs(sumLum / c - prevLum) > luminanceThreshold and (sumLum / c) < 1.0):
return 0
return 1Следует учесть, что при необходимости в проверку Брезенхэма можно внести улучшения, которые приведут к более качественному слиянию сэмплов, например, она может учитывать горизонтальный перенос источников освещения, находящихся по краям изображения. Кроме того, функцию можно легко расширить, чтобы она аппроксимировала площадь источников освещения. Ещё одно улучшение: можно добавить порог расстояния, чтобы не объединять сэмплы, расположенные слишком далеко. Окончательные результаты:
Синим обозначены локальные максимумы кластеров освещения, голубым — окончательные позиции источников освещения, а красным — сэмплы, являющиеся частью того же кластера освещения и соединённые с ним линиями.
Другие примеры результатов:
- Detection of light sources in digital photographs by Maciej Laskowski
- Hammersley Points on the Hemisphere by Holger Dammertz
- Equirectangular Projection by Paul Reed
- Sampling with Hammersley and Halton Points by Tien-Tsin Wong
Автор: PatientZero
