Optymalne wdrożenie wypaczania transportu w Matlab

11

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 ut=1μ0Du1div(u), gdzieuoznacza obrót o 90 stopni przeciwnie do ruchu wskazówek zegara,1dla rozwiązania równania Poissona z warunkami brzegowymi Dirichleta (= 0), iDujest wyznacznikiem macierzy Jakubowej.
- stabilność jest gwarantowana przez czasdt<min|1μ01div(u)|

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.

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 !

WhitAngl
źródło

Odpowiedzi:

5

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)

Uruchomienie testowe zajmuje około 2 sekund.

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

dranxo
źródło
doskonale .. wielkie dzięki! Wypróbuję ten kod i porównuję z moim, aby sprawdzić moje błędy. Podejście to wydaje się być lokalną wersją artykułu Hakera i in. o którym wspomniałem - tj. bez rozwiązania dla laplaciana. Dzięki jeszcze raz !
WhitAngl
W końcu napotykam kilka problemów z tym kodem ...: jeśli obliczę , jestem całkiem daleko od (z ) - nawet przy usuwaniu gaussa plama. Ponadto, jeśli zwiększę liczbę iteracji do kilku tysięcy, kod eksploduje i podaje wartości NaN (i ulega awarii). Dowolny pomysł ? Dzięki ! f2(ϕ)detHϕf1H
WhitAngl
Czy więcej rozmycia pomaga rozwiązać problem NaN?
dranxo
Tak, rzeczywiście, po większej liczbie testów pomaga to uzyskać więcej rozmycia - dzięki! Próbuję teraz 8 kroków z rozmyciem początkowym odchylenia standardowego 140 pikseli, aż do 1 piksela standardowego (wciąż przetwarzam). Nadal mam jednak znaczną ilość oryginalnego obrazu w moim ostatnim wyniku (z rozmyciem 64px). Sprawdzę również, czy nie pozostały zwijające się w . ϕ
WhitAngl
Ok, w porządku. Myślę. Rozmycie występuje, ponieważ obrazy naturalnie nie są ciągłe (krawędzie), a gradient będzie problematyczny. Mam nadzieję, że nadal otrzymujesz dobre odpowiedzi bez zbytniego rozmycia.
dranxo