[Из песочницы] Стабилизация таймлапс-видео на калькуляторе (IPython+OpenCV)20.08.2015 15:48

Подобно многим стихийным и сезонным любителям астрофотографии, в этом августе я ловил ночью Персеиды. Улов небольшой есть, но сейчас не о нём, а о том, что побочным результатом такого лова стала серия фотографий, которые напрашивались на то, чтобы сделать из них таймлапс. Но вот незадача: установка камеры оказалась не столь уж жесткой, как хотелось бы, и между кадрами появилось небольшое смещение. Попытался исправить его плагином дешейкинга в VirtualDub, но результаты не порадовали. Тогда было решено сделать свой велосипед: подробнее о результатах и том как они получены — под катом.
Традиционное «до» и «после» (тут показан небольшой фрагмент). Картинка уменьшена, но даже тут видно «дрожание камеры»:

После обработки:

На чём всё будет делаться: IPython notebook+NumPy+OpenCV.

Необходимое предупреждение: в посте не будет ни нового слова в обработке сигналов, ни нового про означенные язык и библиотеки; разве что новички найдут тут пример «как не надо программировать», но зато «как можно быстро придумать и опробовать свой алгоритм в IPython notebook». Профессионалам же предлагаю полюбоваться на звёздочки.

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

• повышена яркость фона и контраст между ним и тёмными неподвижными объектами
• обрезкой кадра исключены деревья- они уж совсем не образец неподвижности
• имена файлов остались без изменений- просто для удобства

Вот так выглядели предварительно обработанные кадры (справа) по сравнению с исходными (слева):

Чтение и предварительная обработка кадров
Прочитаем исходные файлы (которые были предварительно подготовлены для оценки по ним смещения). Рабочим каталогом для простоты был выбран тот же самый где лежат эти кадры, так что просто прочитаем их имена в массив sampledata. Подробности такого вот обслуживающего кода я не буду приводить в посте, чтобы не загромождать. Их, а также некоторые написанные в процессе работы функции-макросы, можно посмотреть в исходном документе IPython notebook.

Возьмём пару кадров для того чтобы на них всё тестировать (пусть 0 и 4). Покажем эти кадры и разницу между ними простым вычитанием cv2.absdiff ():

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

qx,qy=4,4
k=ones((qx,qy))
im1g=cv2.erode(im1g,k)
im2g=cv2.erode(im2g,k)
show1(im1g,u"после обработки")

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

Оценка сдвига между изображениями
Поиск соответствующих друг другу точек на соседних изображениях сделаем, как показано в примере find_obj.py из opencv: найдём различимые фрагменты с помощью Scale Invariant Feature Transform (SIFT), а потом сопоставим и отфильтруем полученный массив.

detector = cv2.SIFT()
norm = cv2.NORM_L2
matcher = cv2.BFMatcher(norm)

def filter_matches(kp1, kp2, matches, ratio = 0.75):
    mkp1, mkp2 = [], []
    for m in matches:
        if len(m) == 2 and m[0].distance < m[1].distance * ratio:
            m = m[0]
            mkp1.append( kp1[m.queryIdx] )
            mkp2.append( kp2[m.trainIdx] )
    p1 = np.float32([kp.pt for kp in mkp1])
    p2 = np.float32([kp.pt for kp in mkp2])
    kp_pairs = zip(mkp1, mkp2)
    return p1, p2, kp_pairs

def detectandselectmatches(fr1a,fr2a):
    kp1, desc1 = detector.detectAndCompute(fr1a, None)
    kp2, desc2 = detector.detectAndCompute(fr2a, None)
    raw_matches = matcher.knnMatch(desc1, trainDescriptors = desc2, k = 2) #2
    p1, p2, kp_pairs = filter_matches(kp1, kp2, raw_matches)
    return p1, p2

p1, p2 = detectandselectmatches(im1g,im2g)

Полученные массивы p1, p2 представляют собой наборы координат x, y совпадающих точек на 1 и 2 кадре соответственно.

[665.927307129,17.939201355]     [668.513000488,19.468919754]
[744.969177246,60.6581344604]    [747.49786377,61.8129844666]
[746.388549805,77.1945953369]    [749.15411377,78.5462799072]
[892.944763184,169.295532227]    [895.570373535,170.530929565]
[906.57824707,185.634231567]     [908.093933105,186.593307495]
...

Если взять между ними разницу, то получим массив смещений кадра по версии каждой из точек. На этом этапе можно сделать ещё одну фильтрацию (иногда некоторые точки из-за ошибочного соответствия весьма сильно выбиваются из ряда), но если использовать медиану вместо среднего, то все эти ошибочные значения просто не важны. Это хорошо видно на картинке: синим представлено смещение для кажой из пар точек, выбранное значение обозначенно красным, а для сравнения зелёным — простое среднее.

dp=p2-p1
mx,my=np.median(dp[:,0]),np.median(dp[:,1])

dx=2.68393 dy=1.34059 

Выполнение обратного сдвига
Осталось произвести собственно смещение. Построим матрицу афинного преобразования вручную, просто внеся соответствующие величины на нужное место. Можно было бы воспользоваться для этого специальными функциями из opencv, но для вращения и смещения матрица выглядит совсем просто:

где  — угол поворота, а  и  — соответствующие величины смещения

def getshiftmatrix( (dx,dy)):
    return array([[ 1.,  0,   dx],
                  [ 0,  1.,   dy]])

и собственно запустим преобразование

def shiftimg(im2,shift):
    tr1=getshiftmatrix(shift)
    return cv2.warpAffine(im2,tr1,tuple(reversed(im2.shape[:2]) ))
im2r= shiftimg(im2,tuple(-array(shift)))

На результат посмотрим, сравнив, как и в начале, простым вычитанием:

Вуаля! То, что требовалось — на месте всех неподвижных объектов чернота. Значит кадры совместились.

Обработка всех кадров
Всё работает. Можно оформлять обработку… или ещё нет? Остались несколько деталей:

• смещение-то было подсчитано на специально обработанном и обрезанном кадре, а производить его надо уже на основном;
• после смещения по краям кадра останутся черные полосы там, откуда картинку сдвинули. Их можно замазать из соседнего кадра.

Учтём всё это и напишем.

# в переменной basepath4orig сообщим где брать основые кадры для коррекции

arshifts=[]

im1= cv2.imread(sampledata[0]) # base frame
im1g = preprocess(im1)
kp1, desc1 = detector.detectAndCompute(im1g, None)

imgprev=cv2.imread(basepath4orig+sampledata[0]) #base original frame

for i,x in enumerate(sampledata):
    print x,
    im2g = preprocess(cv2.imread(x))    
    kp2, desc2 = detector.detectAndCompute(im2g, None)
    
    raw_matches = matcher.knnMatch(desc1, trainDescriptors = desc2, k = 2) #2
    p1, p2, kp_pairs = filter_matches(kp1, kp2, raw_matches)
    dp=p2-p1
    
    if len(dp)<=0: shift=0,0 
    else: dx,dy=np.median(dp[:,0]),np.median(dp[:,1])
    print dx,dy    
    
    #process original frame
    imgr= shiftimg(cv2.imread(basepath4orig+x),(-dx,-dy))

    if -dy>0: imgr[:int(ceil(abs(dy))),:,:] = imgprev[:int(ceil(abs(dy))),:,:]
    if -dy<0: imgr[-int(ceil(abs(dy))):,:,:] = imgprev[-int(ceil(abs(dy))):,:,:]
    if -dx>0: imgr[:,:int(ceil(abs(dx))),:] = imgprev[:,:int(ceil(abs(dx))),:]
    if -dx<0: imgr[:,-int(ceil(abs(dx))):,:] = imgprev[:,-int(ceil(abs(dx))):,:]
    
    imgprev=imgr
    
    cv2.imwrite('shifted_'+x+'.JPG',imgr)

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

Результат на youtube.com:

Ссылка на копию документа IPython notebook, в котором всё и делалось

Использованные материалы и средства:
Также мои благодарности

© Habrahabr.ru