Optymalne wdrożenie wypaczania transportu w Matlab

Wdrażam artykuł „ Optymalny transport masowy do rejestracji i wypaczania ”, moim celem jest umieszczenie go w Internecie, ponieważ po prostu nie mogę znaleźć żadnego kodu eulerowskiego transportu masowego w Internecie, co byłoby interesujące, przynajmniej dla społeczności naukowej zajmującej się przetwarzaniem obrazu.

Artykuł można podsumować w następujący sposób:
- znajdź początkową mapę $u$ przy użyciu dopasowań histogramu 1D wzdłuż współrzędnych x i y
- rozwiąż dla stałego punktu $u_t = \frac{1}{\mu_0} Du \nabla^\perp\triangle^{-1}div(u^\perp)$, gdzie$u^\perp$oznacza obrót o 90 stopni przeciwnie do ruchu wskazówek zegara,$\triangle^{-1}$dla rozwiązania równania Poissona z warunkami brzegowymi Dirichleta (= 0), i$Du$jest wyznacznikiem macierzy Jakubowej.
- stabilność jest gwarantowana przez czas$dt<\min|\frac{1}{\mu_0}\nabla^\perp\triangle^{-1}div(u^\perp)|$

W przypadku symulacji numerycznych (wykonywanych na regularnej siatce), wskazują one użycie poicalc Matlaba do rozwiązania równania Poissona , używają wyśrodkowanych różnic skończonych dla pochodnych przestrzennych, z wyjątkiem $Du$ który jest obliczany przy użyciu schematu wiatru.

Używając mojego kodu, funkcjonalność energetyczna i zawijanie mapowania odpowiednio zmniejszają się przez kilka iteracji (od kilkudziesięciu do kilku tysięcy w zależności od kroku czasowego). Ale potem symulacja eksploduje: energia wzrasta, aby osiągnąć NAN w bardzo niewielu iteracjach. Próbowałem kilku rozróżnień i integracji (można znaleźć tutaj zamiennik wyższego rzędu dla cumptrapz ) i różnych schematów interpolacji, ale zawsze mam ten sam problem (nawet na bardzo płynnych obrazach, wszędzie niezerowe itp.).
Czy ktoś byłby zainteresowany zapoznaniem się z kodem i / lub teoretycznym problemem, przed którym stoję? Kod jest raczej krótki.

• Kod z funkcjami debugowania

  • funkcja rejestracji
  • kod testowy, pod warunkiem, że masz do zarejestrowania dwa obrazy tego samego rozmiaru

• Tylko niezbędna funkcja bez testów (<100 linii)

Zastąp gradient2 () na końcu gradientem (). To był gradient wyższego rzędu, ale też nie rozwiązuje problemu.

Interesuje mnie tylko optymalna część transportowa papieru, a nie dodatkowy termin regularyzacji.

Dzięki !

Mój dobry przyjaciel Pascal zrobił to kilka lat temu ( prawie w Matlabie):

#! /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)

Metoda opadania gradientu jest wyjaśniona tutaj: people.clarkson.edu/~ebollt/Papers/quadcost.pdf