[Перевод] Распознавание источников освещения на картах окружения

image


В этой статье представлена реализация на Python алгоритма распознавания источников освещения на картах окружения (LDR или HDR) при помощи равнопромежуточной проекции (equirectangular projection). Однако после внесения незначительных изменений её также можно использовать с простыми фоновыми изображениями или кубическими картами. Примеры возможного применения алгоритма: программы трассировки лучей, в которых требуется распознавать первичные источники освещения для испускания из них лучей; в растеризованных рендерерах он может применяться для отбрасывания теней, использующих карту окружения; кроме того, алгоритм также можно применять в программах устранения засветов, например в AR.

Алгоритм состоит из следующих этапов:

  1. Снижение разрешения исходного изображения, например, до 1024.
  2. Преобразование изображения в яркость (luminance), при необходимости с размытием изображения.
  3. Применение метода квази-Монте-Карло.
  4. Преобразование из сферических координат в равнопромежуточные.
  5. Фильтрация сэмплов на основании яркости соседа.
  6. Сортировка сэмплов на основании их яркости.
  7. Фильтрация сэмплов на основании евклидовой метрики.
  8. Слияние сэмплов при помощи алгоритма Брезенхэма.
  9. Вычисление позиции кластера освещения на основании его яркости.


Существует множество алгоритмов снижения разрешения изображений. Билинейная фильтрация — самый быстрый или простой в реализации, к тому же он лучше всего подходит в большинстве случаев. Для преобразования яркости и в 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:

4f9b0770f72e69ede59e56a25258343d.png


Можно использовать и другие методы сэмплирования, например Холтона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]


Это даст нам хорошее распределение сэмплов, которые будут проверяться на наличие источников освещения:

aa80a64703ef6e38cd0b73f9b02a5852.png


Фильтрация источников освещения


В первом проходе фильтрации мы игнорируем все сэмплы, не превосходящие порогового значения яркости (для 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


Получившиеся сэмплы проходят этап слияния, чтобы ещё больше снизить количество источников освещения:

d39152e1e587c40c3ed29888faf02120.png


Слияние источников освещения


На последнем этапе выполняется слияние сэмплов, принадлежащих к одному кластеру освещения. Для этого мы можем использовать алгоритм Брезенхэма и начать с сэмплов с самой большой яркости, потому что они уже упорядочены. Когда мы находим источник освещения, удовлетворяющий проверке Брезенхэма, то используем его позицию для изменения позиции источника на основании веса прогона:

  # 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


Следует учесть, что при необходимости в проверку Брезенхэма можно внести улучшения, которые приведут к более качественному слиянию сэмплов, например, она может учитывать горизонтальный перенос источников освещения, находящихся по краям изображения. Кроме того, функцию можно легко расширить, чтобы она аппроксимировала площадь источников освещения. Ещё одно улучшение: можно добавить порог расстояния, чтобы не объединять сэмплы, расположенные слишком далеко. Окончательные результаты:

4224a8fd98cec4203235bb4678ee6826.png


Синим обозначены локальные максимумы кластеров освещения, голубым — окончательные позиции источников освещения, а красным — сэмплы, являющиеся частью того же кластера освещения и соединённые с ним линиями.

Другие примеры результатов:

a44e484510a596be61f83f9c303f25fc.png



  1. Detection of light sources in digital photographs by Maciej Laskowski
  2. Hammersley Points on the Hemisphere by Holger Dammertz
  3. Equirectangular Projection by Paul Reed
  4. Sampling with Hammersley and Halton Points by Tien-Tsin Wong

© Habrahabr.ru