Оптимальная реализация транспортной деформации в Matlab

Я внедряю документ « Оптимальный массовый транспорт для регистрации и деформации », моя цель - выложить его в Интернет, поскольку я просто не могу найти в Интернете код эйлерова массового транспорта, и это было бы интересно, по крайней мере, для исследовательского сообщества в области обработки изображений.

Работа может быть кратко изложена следующим образом:
- найти исходную карту $u$ используя одномерные гистограммные совпадения по координатам x и y
- решить для фиксированной точки , где обозначает поворот на 90 градусов против часовой стрелки, - решение уравнения Пуассона с граничными условиями Дирихле (= 0), и является определителем матрицы Якоби. - стабильность гарантируется для временного шага dt <\ min | \ frac {1} {\ mu_0} \ nabla ^ \ perp \ triangle ^ {- 1} div (u ^ \ perp) |$u_t = \frac{1}{\mu_0} Du \nabla^\perp\triangle^{-1}div(u^\perp)$$u^\perp$$\triangle^{-1}$$Du$
$dt<\min|\frac{1}{\mu_0}\nabla^\perp\triangle^{-1}div(u^\perp)|$

При численном моделировании (выполняемом на регулярной сетке) они указывают, что для решения уравнения Пуассона используют poicalc от matlab , они используют центрированные конечные разности для пространственных производных, за исключением $Du$ который вычисляется с использованием схемы против ветра.

Используя мой код, функционал энергии и завиток отображения должным образом уменьшаются за пару итераций (от нескольких десятков до нескольких тысяч в зависимости от временного шага). Но после этого симуляция взрывается: энергия увеличивается, чтобы достичь NAN за несколько итераций. Я попробовал несколько заказов на дифференцировки и интеграции (более высокая замену для того , чтобы cumptrapz можно найти здесь ), а также различные схемы интерполяции, но я всегда получаю один и тот же вопрос (даже на очень гладкие изображения, ненулевой везде и т.д.).
Кому-нибудь было бы интересно посмотреть на код и / или на теоретическую проблему, с которой я столкнулся? Код довольно короткий.

• Код с функциями отладки

  • функция регистрации
  • тестовый код, если у вас есть два изображения одинакового размера для регистрации

• Просто необходимая функция без тестовых материалов (<100 строк)

Пожалуйста, замените градиент2 () в конце градиентом (). Это был градиент более высокого порядка, но он тоже ничего не решает.

Сейчас меня интересует только оптимальная транспортная часть статьи, а не дополнительный термин регуляризации.

Спасибо !

Мой хороший друг Паскаль сделал это несколько лет назад (это почти в Matlab):

#! /usr/bin/env python

#from scipy.interpolate import interpolate
from pylab import *
from numpy import *


def GaussianFilter(sigma,f):
    """Apply Gaussian filter to an image"""
    if sigma > 0:
        n = ceil(4*sigma)
        g = exp(-arange(-n,n+1)**2/(2*sigma**2))
        g = g/g.sum()

        fg = zeros(f.shape)

        for i in range(f.shape[0]):
            fg[i,:] = convolve(f[i,:],g,'same')
        for i in range(f.shape[1]):
            fg[:,i] = convolve(fg[:,i],g,'same')
    else:
        fg = f

    return fg


def clamp(x,xmin,xmax):
    """Clamp values between xmin and xmax"""
    return minimum(maximum(x,xmin),xmax)


def myinterp(f,xi,yi):
    """My bilinear interpolator (scipy's has a segfault)"""
    M,N = f.shape
    ix0 = clamp(floor(xi),0,N-2).astype(int)
    iy0 = clamp(floor(yi),0,M-2).astype(int)
    wx = xi - ix0
    wy = yi - iy0
    return ( (1-wy)*((1-wx)*f[iy0,ix0] + wx*f[iy0,ix0+1]) +
        wy*((1-wx)*f[iy0+1,ix0] + wx*f[iy0+1,ix0+1]) )


def mkwarp(f1,f2,sigma,phi,showplot=0):
    """Image warping by solving the Monge-Kantorovich problem"""
    M,N = f1.shape[:2]

    alpha = 1
    f1 = GaussianFilter(sigma,f1)
    f2 = GaussianFilter(sigma,f2)

    # Shift indices for going from vertices to cell centers
    iUv = arange(M)             # Up
    iDv = arange(1,M+1)         # Down
    iLv = arange(N)             # Left
    iRv = arange(1,N+1)         # Right
    # Shift indices for cell centers (to cell centers)
    iUc = r_[0,arange(M-1)]
    iDc = r_[arange(1,M),M-1]
    iLc = r_[0,arange(N-1)]
    iRc = r_[arange(1,N),N-1]
    # Shifts for going from centers to vertices
    iUi = r_[0,arange(M)]
    iDi = r_[arange(M),M-1]
    iLi = r_[0,arange(N)]
    iRi = r_[arange(N),N-1]


    ### The main gradient descent loop ###      
    for iter in range(0,30):
        ### Approximate derivatives ###
        # Compute gradient phix and phiy at pixel centers.  Array phi has values
        # at the pixel vertices.
        phix = (phi[iUv,:][:,iRv] - phi[iUv,:][:,iLv] + 
            phi[iDv,:][:,iRv] - phi[iDv,:][:,iLv])/2
        phiy = (phi[iDv,:][:,iLv] - phi[iUv,:][:,iLv] + 
            phi[iDv,:][:,iRv] - phi[iUv,:][:,iRv])/2
        # Compute second derivatives at pixel centers using central differences.
        phixx = (phix[:,iRc] - phix[:,iLc])/2
        phixy = (phix[iDc,:] - phix[iUc,:])/2
        phiyy = (phiy[iDc,:] - phiy[iUc,:])/2
        # Hessian determinant
        detD2 = phixx*phiyy - phixy*phixy

        # Interpolate f2 at (phix,phiy) with bilinear interpolation
        f2gphi = myinterp(f2,phix,phiy)

        ### Update phi ###
        # Compute M'(phi) at pixel centers
        dM = alpha*(f1 - f2gphi*detD2)
        # Interpolate to pixel vertices
        phi = phi - (dM[iUi,:][:,iLi] + 
            dM[iDi,:][:,iLi] + 
            dM[iUi,:][:,iRi] + 
            dM[iDi,:][:,iRi])/4


    ### Plot stuff ###      
    if showplot:
        pad = 2
        x,y = meshgrid(arange(N),arange(M))
        x = x[pad:-pad,:][:,pad:-pad]
        y = y[pad:-pad,:][:,pad:-pad]
        phix = phix[pad:-pad,:][:,pad:-pad]
        phiy = phiy[pad:-pad,:][:,pad:-pad]

        # Vector plot of the mapping
        subplot(1,2,1)
        quiver(x,y,flipud(phix-x),-flipud(phiy-y))
        axis('image')
        axis('off')
        title('Mapping')

        # Grayscale plot of mapping divergence
        subplot(1,2,2)  
        divs = phixx + phiyy # Divergence of mapping s(x)
        imshow(divs[pad:-pad,pad:-pad],cmap=cm.gray)
        axis('off')
        title('Divergence of Mapping')
        show()

    return phi


if __name__ == "__main__":  # Demo
    from pylab import *
    from numpy import * 

    f1 = imread('brain-tumor.png')
    f2 = imread('brain-healthy.png')
    f1 = f1[:,:,1]
    f2 = f2[:,:,1]

    # Initialize phi as the identity map
    M,N = f1.shape
    n,m = meshgrid(arange(N+1),arange(M+1))
    phi = ((m-0.5)**2 + (n-0.5)**2)/2

    sigma = 3
    phi = mkwarp(f1,f2,sigma,phi)
    phi = mkwarp(f1,f2,sigma/2,phi,1)
#   phi = mkwarp(f1,f2,sigma/4,phi,1)

Метод градиентного спуска объясняется здесь: people.clarkson.edu/~ebollt/Papers/quadcost.pdf