Optimale Implementierung von Transport Warping in Matlab

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 $u_t = \frac{1}{\mu_0} Du \nabla^\perp\triangle^{-1}div(u^\perp)$, wobei$u^\perp$für eine Drehung um 90 Grad gegen den Uhrzeigersinn steht,$\triangle^{-1}$für die Lösung der Poisson-Gleichung mit Dirichlet-Randbedingungen (= 0) und$Du$ist die Determinante der Jacobi-Matrix.
- Stabilität ist für einen Zeitschrittdt<mingarantiert$dt<\min|\frac{1}{\mu_0}\nabla^\perp\triangle^{-1}div(u^\perp)|$

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.

• Code mit Debugging-Funktionen

  • Registrierungsfunktion
  • Testcode vorausgesetzt, Sie müssen zwei Bilder derselben Größe registrieren

• Nur die notwendige Funktion ohne Testmaterial (<100 Zeilen)

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 !

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)

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