from numpy import *
import numpy as np
import cv2
import cv2.cv as cv
import os
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline
import _NEMOpy as npy #небольшая личная библиотека шорткатов. не обращайте внимания :)
Подобно многим стихийным и сезонным любителям астрофотографии, в этом августе я, как обычно ловил ночью персеиды. Улов небольшой есть, но сейчас не о нём, а о том что побочным результатом такого лова стала серия фотографий, которые напрашивались на то чтобы сделать из них таймлапс. Но вот незадача: установка камеры оказалась не столь уж жесткой как хотелось бы и между кадрами появилось небольшое смещение. Попытался исправить его плагином дешейкинга в VirtualDub, но результаты не порадовали. Тогда было решено сделать свой велосипед
%%html
традиционное "до" и "после" (тут показан небольшой фрагмент).
картинка уменьшена в 4 раза, но даже тут видно "дрожание камеры"
<img src="https://habrastorage.org/files/c76/97d/6b6/c7697d6b6e9a4ff0bab94f0955341f9a.gif"/>
и после обработки:
<img src="https://habrastorage.org/files/f7a/bb7/a32/f7abb7a326ff436782f154f42dddb2cc.gif"/>
На чём всё будет делаться: IPython notebook+SciPy+OpenCV
Необходимое предупреждение: в посте не будет ни нового слова в обработке сигналов ни нового про означенные язык и библиотеки- разве что новички найдут тут пример "как не надо программировать", но зато "как можно быстро придумать и опробовать свой алгоритм в IPython notebook". Профессионалам же предлагаю полюбоваться на звёздочки.
Решив отказываться, по возможности, от платных программ для которых есть бесплатные аналоги, я начал искать, среди прочего, замену матлабу. И остановился на связке IPython+SciPy[+OpenCV]. Однако использую их именно в роли большого и очень удобного, но калькулятора: для быстрого прототипирования каких-либо идей и решений или для одноразовой обработки, когда проще самому объяснить компьютеру что мне от него требуется чем искать для этого подходящую программу, которая может ещё оказаться платной или делать немножко не то что мне надо- вот про этот случай я и хочу рассказать в посте.
Для улучшения различимости неподвижных объектов решено было создать специальную версию всех изображений и в дальнейшем находить смещение уже по ней. Что было для этого сделано:
plt.subplot(121)
npy.imshow_(cv2.imread('U:\\_CAM_\\astro\\20150813_perseidy\\ExportJPG77_1080\\DSC_1313.jpg'))
plt.subplot(122)
npy.imshow_(cv2.imread('.\\DSC_1313.jpg'))
npy.mplssz(16)
print "так выглядели предварительно обработанные кадры (справа) по сравнению с исходными (слева)"
так выглядели предварительно обработанные кадры (справа) по сравнению с исходными (слева)
прочитаем исходные файлы (те что были предварительно подготовлены для оценки по ним смещения). Рабочим каталогом для простоты был выбран тот же самый где лежат эти кадры так что...
basepath = ('./')
sampledata=[x for x in os.listdir(basepath) if x.endswith('.jpg') ]
print ", ".join(sampledata[:7]),'... и так всего', len(sampledata),' файла'
DSC_1251.jpg, DSC_1252.jpg, DSC_1253.jpg, DSC_1254.jpg, DSC_1255.jpg, DSC_1256.jpg, DSC_1257.jpg ... и так всего 84 файла
В процессе работы написались для сокращения некоторые функции-макросы. сложим их тут чтобы всегда были доступны. Сюда же инлайнем парочку функций из моей личной библиотеки подобных шорткатов.
def npy_imshow_(img,colormap='gray',szw=8,szh=7,adapt4screensize=(1024,1280),minQlabels=5):
"""unuversal imshow form matplotlib imshow. RGB for BGR and default Grayscale for others
adapt4screensize=(1024,1280) - image decimated to this+ size to fit screen and not overuse memory. Use None to not decimate
minQlabels - (use only for adapt4screensize<>None) - minQlabels for axis adaptive labels
"""
h,w= img.shape[:2]
if colormap=='gray':
if img.ndim==3:
if img.shape[2]==3:
plt.imshow(cv2.cvtColor(img,cv2.cv.CV_BGR2RGB))
else:
plt.imshow(img,colormap)
plt.gcf().set_size_inches(szw,szh)
def show12plus(im1g,im2g,plus,text1='',text2='',text3=''):
plt.subplot(221)
npy_imshow_(im1g)
plt.title(text1)
plt.subplot(222)
npy_imshow_(im2g)
plt.title(text2)
plt.subplot(2,2,(3,4))
npy_imshow_(plus)
plt.title(text3)
plt.gcf().set_size_inches(14,6)
def show1(im,text=''):
plt.subplot(111)
npy_imshow_(im)
plt.title(text)
plt.gcf().set_size_inches(14,7)
npy.mplRc_RusViaFonts()
Возьмём пару кадров для того чтобы на них всё тестировать (пусть 0 и 4). Покажем эти кадры и разницу между ними
im1=cv2.imread(sampledata[0])
im2=cv2.imread(sampledata[4])
im1g=cv2.cvtColor(im1,cv2.cv.CV_BGR2GRAY)
im2g=cv2.cvtColor(im2,cv2.cv.CV_BGR2GRAY)
d2c=cv2.absdiff(im1,im2)
show12plus(im1,im2,cv2.absdiff(im1,im2),sampledata[0],sampledata[4],u"разница кадров")
Сдвиг кадра виден по тому как проявились края неподвижных объектов, а вот движущиеся звёзды- не помощники в сдвига камеры. Так что избавимся по возможности них с помощью операции erode. Размер ядра подобран опытным путём
qx,qy=4,4
k=ones((qx,qy))
im1g=cv2.erode(im1g,k)
im2g=cv2.erode(im2g,k)
show1(im1g,u"после обработки")
npy.mplssz(10)
Видно что на кадре хорошо выделились неподвижные объекты, к которым можно будет привязатсья.
Готовые этапы обработаки удобно сразу оформлять в компактные процедуры
def preprocess(im1):
im1g=cv2.cvtColor(im1,cv2.cv.CV_BGR2GRAY)
qx,qy=4,4
k=ones((qx,qy))
return cv2.erode(im1g,k)
im1g =preprocess(im1)
im2g =preprocess(im2)
Поиск соответствующих друг другу точек на соседних изображениях сделаем как показано в примере из opencv find_obj.py: найдём различимые фрагменты с помощью Scale Invariant Feature Transform (SIFT), а потом сопоставим и отфильтруем полученный массив (© функция filter_matches
непосредственно использована из примера, detectandselectmatches
тоже во многом заимствует его функциональность. Все ©© на них- за соответствующими авторами). Подробно на их работе я останавливаться сейчас не буду- желающие всегда могут посмотреть хелп- там всё написано, и погонять пример - он довольно нагляден.
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 кадре соответственно
for x in array(zip(p1, p2)).reshape(-1,4)[:5]: print "[{0},{1}] \t [{2},{3}]".format(*x)
print "..."
[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])
print u"вычисленное смещение dx=%s dy=%s " % (mx,my)
вычисленное смещение dx=2.68393 dy=1.34059
plt.subplot(121)
plt.title(u"величина смещения каждой точки по x и по y")
plt.scatter(dp[:,0],dp[:,1],alpha=0.6)
plt.scatter([mx],[my],alpha=0.9,color='r')
npy.mplGridSzAndCross()
plt.subplot(122)
plt.scatter(dp[:,0],dp[:,1],alpha=0.25)
plt.scatter([mx],[my],alpha=0.9,color='r')
plt.scatter([np.mean(dp[:,0])],[np.mean(dp[:,0])],alpha=1,color='g')
plt.title(u"то же со вниманием к окрестности вычисленного смещения")
plt.xlim(-4,4); plt.ylim(-4,4);npy.mplGridSzAndCross(11,4)
plt.show()
print "синим представлено смещение для кажой из пар точек, а выбранное значение обозначенно красным."
print "для сравнения зелёным- простое среднее."
синим представлено смещение для кажой из пар точек, а выбранное значение обозначенно красным. для сравнения зелёным- простое среднее.
def detectshift(im1,im2):
p1, p2 = detectandselectmatches(im1,im2)
dp=p2-p1
if len(dp)<=0: return 0,0
mx,my=np.median(dp[:,0]),np.median(dp[:,1])
return mx,my
shift = detectshift(im1g,im2g)
Осталось произвести собственно смещение. Построим матрицу афинного преобразования вручную, просто внеся соответствующие величины на нужное место. Можно было бы воспользоваться для этого специальными функциями из opencv, но для вращения и смещения матрица выглядит совсем просто: $$ M = \begin{bmatrix} cos( \phi ) & -sin (\phi) & dx \\ sin( \phi ) & cos (\phi) & dy \end{bmatrix} $$ где $\phi$ - угол поворота, а $dx$ и $dy$ - соответствующие величины смещения
def getshiftmatrix( (dx,dy)):
return array([[ 1., 0, dx],
[ 0, 1., dy]])
print getshiftmatrix(shift)
[[ 1. 0. 2.68392944] [ 0. 1. 1.34059143]]
и собственно запустим преобразование
def shiftimg(im2,shift):
tr1=getshiftmatrix(shift)
return cv2.warpAffine(im2,tr1,tuple(reversed(im2.shape[:2]) ))
im2r= shiftimg(im2,tuple(-array(shift)))
# show1(im2r)
На результат посмотрим, сравнив как и в начале, простым вычитанием
show12plus(im1,im2r,cv2.absdiff(im1,im2r))
Вуаля! то что требовалось- на месте всех неподвижные объектов чернота. Значит кадры совместились
теперь- всё вместе и на какой-нибудь другой паре кадров- тот же 0 как опорный и случайный
im1= cv2.imread(sampledata[0])
im2= cv2.imread(sampledata[6])
im1g =preprocess(im1)
im2g =preprocess(im2)
shift = detectshift(im1g,im2g)
print shift
im2r= shiftimg(im2,tuple(- array(shift)))
show12plus(im1,im2r,cv2.absdiff(im1,im2r))
(2.7861176, 1.7378769)
работает. можно оформлять обработку... или ещё нет?
остались несколько деталей
сообщим где брать основые кадры для коррекции
basepath4orig = ('U:\\_CAM_\\astro\\20150813_perseidy\\ExportJPG77_1080\\')
Выберем один кадр опорым (пусть первый) и в цикле пройдём по всем кадрам вычисляя смещение и выполняя сдвиг. Не страшно что опорный кадр будет тоже обработан - сравние с собой даст нулевое смещение. По вычисленным параметрам произведём сдвиг уже оригинального кадра, а чёрные края после смещения заменим значениями с того же места в предыдущем кадре.
%%time
# в переменной 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)
# print arshifts
DSC_1251.jpg 0.0 0.0 DSC_1252.jpg 2.77911 0.627548 DSC_1253.jpg 2.18556 -2.26859 DSC_1254.jpg 2.31195 -3.30951 DSC_1255.jpg 2.68393 1.34059 DSC_1256.jpg 2.10768 -2.89339 DSC_1257.jpg 2.78612 1.73788 DSC_1258.jpg 1.98948 -2.46089 DSC_1259.jpg 1.97647 -2.32199 DSC_1260.jpg 1.97235 -2.29891 DSC_1261.jpg 1.91949 -2.29082 DSC_1262.jpg 2.6347 2.30543 DSC_1263.jpg 2.60498 2.51817 DSC_1264.jpg 2.43317 2.54622 DSC_1265.jpg 2.49084 2.45909 DSC_1266.jpg 1.75325 -1.89879 DSC_1267.jpg 2.63647 2.8098 DSC_1268.jpg 1.72272 -1.63324 DSC_1269.jpg 2.46857 2.83755 DSC_1270.jpg 2.42548 3.06012 DSC_1271.jpg 1.67932 -1.49174 DSC_1272.jpg 1.88031 -1.49158 DSC_1273.jpg 1.64722 -1.50623 DSC_1274.jpg 1.70966 -1.22985 DSC_1275.jpg 1.61435 -1.22524 DSC_1276.jpg 2.37256 3.36833 DSC_1277.jpg 1.60156 -1.14243 DSC_1278.jpg 1.3515 -0.94507 DSC_1279.jpg 2.18668 3.91653 DSC_1280.jpg 2.05737 3.70058 DSC_1281.jpg 1.39325 -0.930043 DSC_1282.jpg 1.28571 0.31342 DSC_1283.jpg 1.43927 -0.676718 DSC_1284.jpg 2.01013 4.19549 DSC_1285.jpg 1.88617 3.9442 DSC_1286.jpg 1.93124 3.85625 DSC_1287.jpg 1.98816 3.84153 DSC_1288.jpg 1.79617 4.03269 DSC_1289.jpg 1.32611 -0.558083 DSC_1290.jpg 1.33246 -0.55442 DSC_1291.jpg 1.18536 -0.529177 DSC_1292.jpg 1.83032 4.1209 DSC_1293.jpg 1.86185 3.94077 DSC_1294.jpg 2.06213 4.14894 DSC_1295.jpg 1.6629 4.11497 DSC_1296.jpg 1.09 -0.308632 DSC_1297.jpg 1.85065 4.5742 DSC_1298.jpg 1.1192 -0.244278 DSC_1299.jpg 1.11206 -0.150941 DSC_1300.jpg 1.65155 4.14029 DSC_1301.jpg 1.67245 4.1913 DSC_1302.jpg 0.948151 0.186799 DSC_1303.jpg 1.63101 4.44316 DSC_1304.jpg 1.93103 4.97454 DSC_1305.jpg 1.10663 0.179092 DSC_1306.jpg 1.0569 0.228849 DSC_1307.jpg 1.95453 5.04773 DSC_1308.jpg 1.69073 4.38213 DSC_1309.jpg 0.972473 -0.106834 DSC_1310.jpg 1.90158 4.84666 DSC_1311.jpg 1.54582 4.75343 DSC_1312.jpg 1.65399 5.24611 DSC_1313.jpg 1.57697 4.61638 DSC_1314.jpg 1.42847 5.01108 DSC_1315.jpg 0.992493 0.460233 DSC_1316.jpg 1.07025 0.457748 DSC_1317.jpg 0.999939 0.264526 DSC_1318.jpg 2.00232 5.29204 DSC_1319.jpg 0.795135 0.533727 DSC_1320.jpg 0.868713 0.388855 DSC_1321.jpg 2.10025 5.46939 DSC_1322.jpg 1.72119 5.45898 DSC_1323.jpg 1.58069 4.97898 DSC_1324.jpg 1.29068 5.49902 DSC_1325.jpg 0.780212 0.647081 DSC_1326.jpg 0.680725 0.909674 DSC_1327.jpg 1.52545 5.47639 DSC_1328.jpg 1.77603 5.55197 DSC_1329.jpg 1.32104 5.06866 DSC_1330.jpg 1.81339 5.39985 DSC_1331.jpg 1.43182 5.28482 DSC_1332.jpg 0.768524 0.941899 DSC_1333.jpg 1.64362 5.79778 DSC_1334.jpg 1.66406 5.69852 Wall time: 29.7 s
Вот и всё- результат получен, неподвижные объекты- как прикленнные, звёзды вертятся вокруг земли :) как им и положено, облака плывут по своим делам, а я раздумываю что бы ещё такого снять- уже специально для таймлапса.
%%html
<a href="http://www.youtube.com/watch?v=j4JJIajjrAc">результат на youtube.com</a>
<iframe width="640" height="480" src="https://www.youtube.com/embed/j4JJIajjrAc?rel=0&showinfo=0" frameborder="0" allowfullscreen></iframe>
самые внимательные могли заметить что в youtube попала предыдущая версия- без коррекции чёрных полос по краям- заменять уже не буду, в gifках в посте уже дана нормальная версия