Optimale Implementierung von Transport Warping in Matlab

11

Ich implementiere das Papier " Optimaler Massentransport für Registrierung und Warping ". Mein Ziel ist es, es online zu stellen, da ich online keinen eulerschen Massentransportcode finden kann. Dies wäre zumindest für die Forschungsgemeinschaft in der Bildverarbeitung interessant.

Das Papier kann wie folgt zusammengefasst werden:
- Finden einer Anfangskarte u Verwendung von 1D-Histogramm-Übereinstimmungen entlang der x- und y-Koordinaten
- Auflösen nach dem festen Punkt von ut=1μ0Du1div(u), wobeiufür eine Drehung um 90 Grad gegen den Uhrzeigersinn steht,1für die Lösung der Poisson-Gleichung mit Dirichlet-Randbedingungen (= 0) undDuist die Determinante der Jacobi-Matrix.
- Stabilität ist für einen Zeitschrittdt<mingarantiertdt<min|1μ01div(u)|

Für numerische Simulationen (die auf einem regulären Gitter durchgeführt werden) geben sie an, dass zur Lösung der Poisson-Gleichung matlabs Poicalc verwendet wird. Für räumliche Ableitungen verwenden sie zentrierte endliche Differenzen, mit Ausnahme von Du das unter Verwendung eines Aufwindschemas berechnet wird.

Bei Verwendung meines Codes nehmen die Energiefunktion und die Kräuselung des Mappings für einige Iterationen ordnungsgemäß ab (von einigen zehn auf einige tausend, abhängig vom Zeitschritt). Danach explodiert die Simulation: Die Energie steigt, um in sehr wenigen Iterationen ein NAN zu erreichen. Ich habe mehrere Ordnungen für die Differenzierungen und Integrationen (ein Ersatz für Cumptrapz höherer Ordnung finden Sie hier ) und verschiedene Interpolationsschemata ausprobiert, aber ich habe immer das gleiche Problem (auch bei sehr glatten Bildern, überall ungleich Null usw.).
Möchte sich jemand den Code und / oder das theoretische Problem ansehen, mit dem ich konfrontiert bin? Der Code ist ziemlich kurz.

Bitte ersetzen Sie gradient2 () am Ende durch gradient (). Dies war ein Gradient höherer Ordnung, löst aber auch keine Probleme.

Ich interessiere mich vorerst nur für den optimalen Transportteil des Papiers, nicht für den zusätzlichen Regularisierungsbegriff.

Vielen Dank !

WhitAngl
quelle

Antworten:

5

Mein guter Freund Pascal hat das vor ein paar Jahren gemacht (es ist fast in 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)

Testlauf dauert ca. 2 Sekunden.

Der Gradientenabstiegsansatz wird hier erläutert: people.clarkson.edu/~ebollt/Papers/quadcost.pdf

dranxo
quelle
ausgezeichnet .. vielen Dank! Ich werde diesen Code ausprobieren und mit meinem vergleichen, um nach meinen Fehlern zu suchen. Dieser Ansatz scheint tatsächlich die lokale Version des Papiers von Haker et al. das habe ich erwähnt - dh ohne nach einem Laplace zu suchen. Danke noch einmal !
WhitAngl
Ich habe endlich ein paar Probleme mit diesem Code ...: Wenn ich berechne, bin ich ziemlich weit von (mit dem Hessischen) - selbst wenn ich den Gaußschen entferne verwischen. Wenn ich nur die Anzahl der Iterationen auf ein paar Tausend erhöhe, explodiert der Code und gibt NaN-Werte aus (und stürzt ab). Irgendeine Idee ? Vielen Dank ! f 1 H.f2(ϕ)detHϕf1H
WhitAngl
Hilft mehr Unschärfe beim NaN-Problem?
Dranxo
Ja, nach mehr Tests hilft es tatsächlich bei mehr Unschärfe - danke!. Ich versuche jetzt 8 Schritte mit einer Startunschärfe von 140 Pixel Standardabweichung bis 1 Pixel stdev (noch zu berechnen). Ich habe jedoch immer noch eine erhebliche Menge des Originalbilds in meinem letzten Ergebnis (mit 64px Unschärfe). Ich werde auch nach verbleibenden Locken in . ϕ
WhitAngl
Oh ok gut. Ich denke. Die Unschärfe ist da, da Bilder natürlich nicht kontinuierlich sind (Kanten) und der Farbverlauf problematisch ist. Hoffentlich erhalten Sie immer noch gute Antworten, ohne zu viel zu verwischen.
Dranxo