#!/usr/bin/env python
# coding: utf-8

# In[2]:


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
get_ipython().run_line_magic('matplotlib', 'inline')
import _NEMOpy as npy #небольшая личная библиотека шорткатов. не обращайте внимания :)


# #Стабилизация таймлапс-видео калькуляторе (IPython+OpenCV)

# ### Постановка задачи

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

# In[23]:


get_ipython().run_cell_magic('html', '', 'традиционное "до" и "после" (тут показан небольшой фрагмент). \nкартинка уменьшена в 4 раза, но даже тут видно "дрожание камеры"\n<img src="https://habrastorage.org/files/c76/97d/6b6/c7697d6b6e9a4ff0bab94f0955341f9a.gif"/>\n\nи после обработки:\n<img src="https://habrastorage.org/files/f7a/bb7/a32/f7abb7a326ff436782f154f42dddb2cc.gif"/>\n')


# **На чём всё будет делаться**: IPython notebook+SciPy+OpenCV
# 
# Необходимое предупреждение: в посте не будет ни нового слова в обработке сигналов ни нового про означенные язык и библиотеки- разве что новички найдут тут пример "как не надо программировать", но зато "как можно быстро придумать и опробовать свой алгоритм в IPython notebook". Профессионалам же предлагаю полюбоваться на звёздочки.

# #####Почему "калькулятор",а также о том зачем мне было делать свой велосипед- лирическое отступление 
# Решив отказываться, по возможности, от платных программ для которых есть бесплатные аналоги, я начал искать, среди прочего, замену матлабу. И остановился на связке IPython+SciPy[+OpenCV]. Однако использую их именно в роли большого и очень удобного, но калькулятора: для быстрого прототипирования каких-либо идей и решений или для одноразовой обработки, когда проще самому объяснить компьютеру что мне от него требуется чем искать для этого подходящую программу, которая может ещё оказаться платной или делать немножко не то что мне надо- вот про этот случай я и хочу рассказать в посте.
# </spoiler>

# ### Предварительная подготовка данных

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

# In[4]:


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 "так выглядели предварительно обработанные кадры (справа) по сравнению с исходными (слева)"


# ### Чтение и подготовка данных

# прочитаем исходные файлы (те что были предварительно подготовлены для оценки по ним смещения). Рабочим каталогом для простоты был выбран тот же самый где лежат эти кадры так что...

# In[5]:


basepath = ('./')
sampledata=[x for x in os.listdir(basepath) if x.endswith('.jpg') ]
print ", ".join(sampledata[:7]),'... и так всего', len(sampledata),' файла'


# В процессе работы написались для сокращения некоторые функции-макросы. сложим их тут чтобы всегда были доступны. Сюда же инлайнем парочку функций из моей личной библиотеки подобных шорткатов.

# In[6]:


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). Покажем эти кадры и разницу между ними

# In[7]:


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](http://docs.opencv.org/doc/tutorials/imgproc/erosion_dilatation/erosion_dilatation.html). Размер ядра подобран опытным путём

# In[8]:


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


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

# In[9]:


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](https://github.com/Itseez/opencv/blob/master/samples/python2/find_obj.py): найдём различимые фрагменты с помощью [Scale Invariant Feature Transform (SIFT)](http://docs.opencv.org/modules/nonfree/doc/feature_detection.html?highlight=sift#sift), а потом сопоставим и отфильтруем полученный массив (© функция `filter_matches`  непосредственно использована из [примера](https://github.com/Itseez/opencv/blob/master/samples/python2/find_obj.py), `detectandselectmatches` тоже во многом заимствует его функциональность. Все ©© на них- за соответствующими авторами). Подробно на их работе я останавливаться сейчас не буду- желающие всегда могут посмотреть [хелп](http://docs.opencv.org/modules/nonfree/doc/feature_detection.html?highlight=sift#sift)- там всё написано, и погонять [пример](https://github.com/Itseez/opencv/blob/master/samples/python2/find_obj.py) - он довольно нагляден.

# In[10]:


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


# In[11]:


p1, p2 = detectandselectmatches(im1g,im2g)


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

# In[12]:


for x in array(zip(p1, p2)).reshape(-1,4)[:5]: print "[{0},{1}] \t [{2},{3}]".format(*x)
print "..."


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

# In[13]:


dp=p2-p1
mx,my=np.median(dp[:,0]),np.median(dp[:,1])
print u"вычисленное смещение dx=%s dy=%s " % (mx,my)


# In[14]:


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 "для сравнения зелёным- простое среднее."


# In[15]:


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)


# ### Выполнение обратного сдвига

# Осталось произвести собственно смещение. Построим матрицу [афинного преобразования](http://docs.opencv.org/modules/imgproc/doc/geometric_transformations.html?highlight=warpaffine) вручную, просто внеся соответствующие величины на нужное место. Можно было бы воспользоваться для этого специальными функциями из opencv, но для вращения и смещения матрица выглядит совсем просто:
# $$
# M = \begin{bmatrix}
#  cos( \phi ) &  -sin (\phi) &  dx  \\ 
#  sin( \phi ) &  cos (\phi) &  dy
# \end{bmatrix}
# $$
# где  $\phi$ - угол поворота, а $dx$ и $dy$ - соответствующие величины смещения

# In[16]:


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


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

# In[17]:


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)


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

# In[18]:


show12plus(im1,im2r,cv2.absdiff(im1,im2r))


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

# теперь- всё вместе и на какой-нибудь другой паре кадров- тот же 0 как опорный и случайный 

# In[19]:


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))


# работает. можно оформлять обработку... или ещё нет?

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

# ### Обработка всех кадров

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

# In[20]:


basepath4orig = ('U:\\_CAM_\\astro\\20150813_perseidy\\ExportJPG77_1080\\')


# Выберем один кадр опорым (пусть первый) и в цикле пройдём по всем кадрам вычисляя смещение и выполняя сдвиг. Не страшно что опорный кадр будет тоже обработан - сравние с собой даст нулевое смещение. По вычисленным параметрам произведём сдвиг уже оригинального кадра, а чёрные края после смещения заменим значениями с того же места в предыдущем кадре.

# In[21]:


get_ipython().run_cell_magic('time', '', "\n# в переменной basepath4orig сообщим где брать основые кадры для коррекции\n\narshifts=[]\n\nim1= cv2.imread(sampledata[0]) # base frame\nim1g = preprocess(im1)\nkp1, desc1 = detector.detectAndCompute(im1g, None)\n\nimgprev=cv2.imread(basepath4orig+sampledata[0]) #base original frame\n\nfor i,x in enumerate(sampledata):\n    print x,\n    im2g = preprocess(cv2.imread(x))    \n    kp2, desc2 = detector.detectAndCompute(im2g, None)\n    \n    raw_matches = matcher.knnMatch(desc1, trainDescriptors = desc2, k = 2) #2\n    p1, p2, kp_pairs = filter_matches(kp1, kp2, raw_matches)\n    dp=p2-p1\n    \n    if len(dp)<=0: shift=0,0 \n    else: dx,dy=np.median(dp[:,0]),np.median(dp[:,1])\n    print dx,dy    \n    \n    #process original frame\n    imgr= shiftimg(cv2.imread(basepath4orig+x),(-dx,-dy))\n\n    if -dy>0: imgr[:int(ceil(abs(dy))),:,:] = imgprev[:int(ceil(abs(dy))),:,:]\n    if -dy<0: imgr[-int(ceil(abs(dy))):,:,:] = imgprev[-int(ceil(abs(dy))):,:,:]\n    if -dx>0: imgr[:,:int(ceil(abs(dx))),:] = imgprev[:,:int(ceil(abs(dx))),:]\n    if -dx<0: imgr[:,-int(ceil(abs(dx))):,:] = imgprev[:,-int(ceil(abs(dx))):,:]\n    \n    imgprev=imgr\n    \n    cv2.imwrite('shifted_'+x+'.JPG',imgr)\n    \n# print arshifts \n")


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

# In[22]:


get_ipython().run_cell_magic('html', '', '<a href="http://www.youtube.com/watch?v=j4JJIajjrAc">результат на youtube.com</a>\n<iframe width="640" height="480" src="https://www.youtube.com/embed/j4JJIajjrAc?rel=0&amp;showinfo=0" frameborder="0" allowfullscreen></iframe>\n')


# самые внимательные могли заметить что в youtube попала предыдущая версия- без коррекции чёрных полос по краям- заменять уже не буду, в gifках в посте уже дана нормальная версия

# In[ ]: