Picamera и фотограмметрия

В прошлом году ко мне несколько раз обращались за консультацией по вопросу использования Raspberry Pi для получения 3D моделей. Хотя поставившие вопрос исчезли не попрощавшись, задача мне показалась интересной, я начал экспериментировать. В общем виде проблему можно сформулировать как подбор оптимальных объективов и точек съемки для получения серии снимков, пригодных для фотограмметрии. Предполагалось, что будет использовано несколько камер и система для их перемещения. Меня в первую очередь интересовала установка и методика для поиска оптимального расположения камер и целесообразность использования объективов рыбий глаз для уменьшения габаритов установки. Сколько камер и как автоматизировать процесс. Поскольку потенциальные заказчики исчезли, то для работы предполагалось использовать только имеющееся оборудование и программное обеспечение с открытым исходным кодом. В результате, экспериментальная установка была создана, и ниже приведенная статья является для меня шпаргалкой по решению возникавших проблем по ее эксплуатации.

В предыдущих двух описанных мной сканерах в статье «3D сканеры Ciclop и Piclop» использовался вращающийся стол. Однако, это не самое удобное устройство для получения серии фотографий для фотограмметрии. С ним у нас получаются снимки объекта с разных сторон, но фон при этом остается постоянным и это сбивает программу с толку. Необходимо либо делать очень ярко освещенный фон, который будет восприниматься на фотографии как равномерно белый, либо его ретушировать, обрабатывая каждую фотографию. Поэтому для получения серии снимков была создана из имеющихся фото кубиков для панорамной съемки установка, при которой камера вращалась вокруг объекта. Установка управлялась по Bluetooth с телефона с помощью программы на De Re BASIC! (с тех пор у языка сменился главный разработчик и в Play Маркет это просто BASIC!) приведенной в моей старой статье. Управление камерой может осуществляться любой из программ, описанных и приведенных в статье «Picamera и дистанционная съемка с живой картинкой».

Объективов для Raspberry Pi у меня много, однако только на некоторых из них написано фокусное расстояние. Поэтому чтобы приблизительно определить фокусное расстояние, я воспользовался возможностями калибровки по снимкам шахматной доски из пакета opencv. Тут возможны два варианта: можно воспользоваться функцией cv2.fisheye.calibrate и получить фокусное расстояние в пикселях и 4 коэффициента дисторсии относительно идеальной геометрии рыбий глаз, или функцией cv2.calibrateCamera, которая кроме фокусного расстояния дает нам 5 коэффициентов дисторсии по модели Брауна-Конради. Мы знаем размер матрицы в мм и размер снимка в пикселях и, следовательно, можем легко вычислить фокусное расстояние в мм. Для большинства объективов обе программы дают близкий результат для фокусного расстояния и, естественно, совершенно разные коэффициенты.

# -*- coding: utf-8 -*-
import yaml
import cv2
assert cv2.__version__[0] == '3', 'The fisheye module requires opencv version >= 3.0.0'
import numpy as np
import glob
CHECKERBOARD = (7,9)
subpix_criteria = (cv2.TERM_CRITERIA_EPS+cv2.TERM_CRITERIA_MAX_ITER, 30, 0.1)
calibration_flags = cv2.fisheye.CALIB_RECOMPUTE_EXTRINSIC+cv2.fisheye.CALIB_CHECK_COND+cv2.fisheye.CALIB_FIX_SKEW
objp = np.zeros((1, CHECKERBOARD[0]*CHECKERBOARD[1], 3), np.float32)
objp[0,:,:2] = np.mgrid[0:CHECKERBOARD[0], 0:CHECKERBOARD[1]].T.reshape(-1, 2)
_img_shape = None
objpoints = [] # 3d point in real world space
imgpoints = [] # 2d points in image plane.
images = glob.glob('./data/mebful/*.*')
for fname in images:
    img = cv2.imread(fname)
    if _img_shape == None:
        _img_shape = img.shape[:2]
    else:
        assert _img_shape == img.shape[:2], "All images must share the same size."
    gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
    # Find the chess board corners
    ret, corners = cv2.findChessboardCorners(gray, CHECKERBOARD, cv2.CALIB_CB_ADAPTIVE_THRESH+cv2.CALIB_CB_FAST_CHECK+cv2.CALIB_CB_NORMALIZE_IMAGE)
    # If found, add object points, image points (after refining them)
    if ret == True:
        objpoints.append(objp)
        cv2.cornerSubPix(gray,corners,(3,3),(-1,-1),subpix_criteria)
        imgpoints.append(corners)
N_OK = len(objpoints)
K = np.zeros((3, 3))
D = np.zeros((4, 1))
rvecs = [np.zeros((1, 1, 3), dtype=np.float64) for i in range(N_OK)]
tvecs = [np.zeros((1, 1, 3), dtype=np.float64) for i in range(N_OK)]
rms, _, _, _, _ = \
    cv2.fisheye.calibrate(
        objpoints,
        imgpoints,
        gray.shape[::-1],
        K,
        D,
        rvecs,
        tvecs,
        calibration_flags,
        (cv2.TERM_CRITERIA_EPS+cv2.TERM_CRITERIA_MAX_ITER, 30, 1e-6)
    )
print("Found " + str(N_OK) + " valid images for calibration")
print("DIM=" + str(_img_shape[::-1]))
#print("K=np.array(" + str(K.tolist()) + ")")
print("D=np.array(" + str(D.tolist()) + ")")
print(K)

f= (K[0,0]+K[1,1])/2
print(f)
#RaspberryPi;RP_OV5647;3.76;
print(_img_shape[1])
fmm=3.76/_img_shape[1]*f
print("F="+str(fmm)+" mm")

#!/usr/bin/env python

'''
camera calibration for distorted images with chess board samples
reads distorted images, calculates the calibration 

usage:
    calibrate.py [--debug ] [--square_size] []

'''

# Python 2/3 compatibility
from __future__ import print_function

import numpy as np
import cv2

# local modules
from common import splitfn

# built-in modules
import os

if __name__ == '__main__':
    import sys
    import getopt
    from glob import glob

    args, img_mask = getopt.getopt(sys.argv[1:], '', ['debug=', 'square_size='])
    args = dict(args)
    args.setdefault('--debug', './output/')
    args.setdefault('--square_size', 1.0)
    if not img_mask:
        img_mask = './data/mebful/*.*'  # default
    else:
        img_mask = img_mask[0]

    img_names = glob(img_mask)
    debug_dir = args.get('--debug')
    if not os.path.isdir(debug_dir):
        os.mkdir(debug_dir)
    square_size = float(args.get('--square_size'))

    pattern_size = (7, 9)
    pattern_points = np.zeros((np.prod(pattern_size), 3), np.float32)
    pattern_points[:, :2] = np.indices(pattern_size).T.reshape(-1, 2)
    pattern_points *= square_size

    obj_points = []
    img_points = []
    h, w = 0, 0
    img_names_undistort = []
    for fn in img_names:
        print('processing %s... ' % fn, end='')
        img = cv2.imread(fn, 0)
        if img is None:
            print("Failed to load", fn)
            continue

        h, w = img.shape[:2]
        found, corners = cv2.findChessboardCorners(img, pattern_size)
        if found:
            term = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_COUNT, 30, 0.1)
            cv2.cornerSubPix(img, corners, (5, 5), (-1, -1), term)

        if not found:
            print('chessboard not found')
            continue

        img_points.append(corners.reshape(-1, 2))
        obj_points.append(pattern_points)

        print('ok')

    # calculate camera distortion
    rms, camera_matrix, dist_coefs, rvecs, tvecs = cv2.calibrateCamera(obj_points, img_points, (w, h), None, None)

    print("\nRMS:", rms)
    print("camera matrix:\n", camera_matrix)
    print("distortion coefficients: ", dist_coefs.ravel())
    K=camera_matrix
    f= (K[0,0]+K[1,1])/2
    print(f)
    #RaspberryPi;RP_OV5647;3.76;
    print(img.shape[1])
    fmm=3.76/img.shape[1]*f
    print("F="+str(fmm)+" mm")

В программе meshroom тоже есть возможность провести калибровку камеры выбрав в редакторе графов, при нажатой правой кнопке мыши, пункт CameraCalibration. Далее надо задать папку, в которой лежат снимки шахматной доски и число клеток в ней. Можно задать и число коэффициентов дисторсии, но в моей версии программа на это не реагирует и вычисляет только 3 коэффициента. Результаты калибровки будут записаны в файл /MeshroomCache/CameraCalibration/XXXX/cameraCalibration.cal.txt.

В дальнейшем, фокусное расстояние и коэффициенты надо будет занести вручную в CameraInit. AliceVision использует фокусное расстояние как ориентир для дальнейшего уточнения и по умолчанию попытается взять фокусное расстояние в мм из заголовка Exif. По умолчанию Raspberry Pi записывает в Exif данные о матрицы, а объектив считает стандартным с фокусным расстоянием 3,6 мм. Реальное значение объектива можно записать программой exiftool, запустив ее в папке с фотографиями. Если используется несколько камер, то надо задать каждой свой серийный номер.

Поскольку моему компьютеру уже более 10 лет: (Processor : 2x Intel(R) Core(TM)2 Duo CPU E8500 @ 3.16GHz, Memory : 6109MB, GPU : GeForce GTS 450), то предполагалось использовать программу Meshroom/AliceVision версию от 2019 года, поскольку версия 2020 года уже отказывался работать с моей видеокартой. На сайте Хабр есть хороший перевод статьи John Hable "AliceVision: фотограмметрия из командной строки". Там описывается работа с версией 2018 года, она мало отличается по интерфейсу от следующей. Некоторые замеченные мной отличия я отмечу ниже. Программа Meshroom является графическим интерфейсом, написанным на Python для пакета AliceVision. Запускается командой LD_LIBRARY_PATH=/usr/lib/nvidia-384 ./Meshroom, где nvidia-ХХХ папка с драйверами на вашем компьютере. Выполняемые команды отображаются в виде прямоугольников в окне программы. Щелчок по прямоугольнику открывает окно настройки параметров. Можно нажать кнопку Start и начнется последовательное выполнение команд. Команду можно выполнить отдельно нажав правую клавишу мыши на соответствующем прямоугольнике и выбрав в меню пункт Compute. Дополнительные задачи могут быть добавлены, если нажать правую клавишу мыши на пустом месте и из выпадающего меню выбрать необходимую команду, например, CameraCalibration.

CameraInit. Интересен для нас тем, что позволяет задать уточнить фокусное расстояние, тип объектива, количество параметров дисторсии.

FeatureExtraction. Извлекает из изображения характерные черты. По-английски это SIFT (Scale invariant feature transform) или по-русски масштабно-инвариантная трансформация признаков. После выполнения операции, если нажать на кнопочку Display Features, то можно увидеть и точки, которые программа выбрала как характерные. Можно выбрать способ отображения в виде точек или повернутых квадратов разных размеров. И таким образом оценить, где могут понадобится дополнительные снимки и, возможно, проекция сетки.

ImageMatching. Вычисляет, какие из изображений логично сопоставлять друг с другом.

FeatureMatching. Ищет соответствия между изображениями.

StructureFromMotion. Важный пункт для оценки пригодности добавленных снимков. После его выполнения программа показывает рассчитанное положение камер и положения характерных точек в пространстве.

PrepareDenseScene. Устраняет дисторсию в снимках и складывает в папку, где их можно посмотреть и оценить. Адрес папки в последней строке свойств.

DepthMap. Вычисляет карту глубин для каждого снимка. Этот процесс выполняется GPU, и температура моей видеокарты при этом поднимается на 30 градусов до 73 градусов.

В результате получается два снимка в формате EXR: ХХХdepthMap.exr и ХХХ simMap.exr. При преобразовании в JPEG файла depthMap.exr я отредактировал уровни для восприятия глазом.

DepthMapFilter. Карты глубин разных снимков перекрываются. В этом пункте осуществляется принудительное согласование.

Eсли у нас достаточно много характерных точек и они равномерно распределены, можно обойтись и без карты глубин и сразу попытаться построить трехмерную модель. В этом случае области, для которых могла быть построена корректная карта глубин, но для которых нет характерных точек, пропадут. Однако, для предварительной быстрой оценки это может быть полезно.

Meshing. Генерируется полигональная сетка, называемая на жаргоне меш (от английского polygon mesh). Дважды шелкнув по прямоугольнику, ее можно рассмотреть и оценить.

MeshFiltering. Сглаживание, устранение больших треугольников.

MeshDecimate. Дополнительная операция, упрощающая сетку.

Сравним финальные результаты в случаях, когда мы использовали карту глубин и когда нет.

Таким образом, есть установка для съемки и программа для получения 3D модели. Я сделал и обработал более 100 снимков. Теперь, убирая снимки, надо выяснить, какие из них избыточны и практически не влияют на качество получаемой модели. Кроме того, вероятно, не все заданные параметры обработки оптимальны. Для выяснения всего этого нужны большие вычислительные мощности. Если я не обновлю свою технику, то скорых результатов и продолжения этой статьи ждать не приходится. В моих дальнейших планах также эксперименты с проецируемой сеткой для объектов, на поверхности которых трудно выделить характерные точки.