Implementasi warping transportasi optimal di Matlab

11

Saya mengimplementasikan makalah " Transportasi Massal Optimal untuk Registrasi dan Warping ", tujuan saya adalah untuk membuatnya online karena saya tidak dapat menemukan kode transportasi massal euler online dan ini akan menarik setidaknya untuk komunitas riset dalam pemrosesan gambar.

Makalah ini dapat diringkas sebagai berikut:
- temukan peta awal u menggunakan pencocokan histogram 1D di sepanjang koordinat x dan y
- pecahkan titik tetap ut=1μ0Du1div(u)u1Du
dt<min|1μ01div(u)|

Untuk simulasi numerik (dilakukan pada grid biasa), mereka mengindikasikan menggunakan matlab's poicalc untuk menyelesaikan persamaan poisson, mereka menggunakan perbedaan hingga terpusat untuk turunan spasial, kecuali untuk Du yang dihitung menggunakan skema arah angin.

Menggunakan kode saya, energi fungsional dan ikal pemetaan berkurang dengan benar untuk beberapa iterasi (dari beberapa puluh hingga beberapa ribu tergantung pada langkah waktu). Tetapi setelah itu, simulasi meledak: energi meningkat untuk mencapai NAN dalam iterasi yang sangat sedikit. Saya mencoba beberapa pesanan untuk diferensiasi dan integrasi (pengganti orde yang lebih tinggi untuk cumptrapz dapat ditemukan di sini ), dan skema interpolasi yang berbeda, tetapi saya selalu mendapatkan masalah yang sama (bahkan pada gambar yang sangat halus, bukan nol di mana-mana dll).
Adakah yang akan tertarik melihat kode dan / atau masalah teoretis yang saya hadapi? Kode ini agak pendek.

Silakan ganti gradient2 () di bagian akhir dengan gradien (). Ini adalah gradien urutan yang lebih tinggi tetapi juga tidak menyelesaikan masalah.

Saya hanya tertarik pada bagian pengangkutan kertas yang optimal untuk saat ini, bukan istilah regularisasi tambahan.

Terima kasih!

WhitAngl
sumber

Jawaban:

5

Teman baik saya Pascal membuat ini beberapa tahun yang lalu ( hampir di 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)

Uji coba, membutuhkan waktu sekitar 2 detik.

Pendekatan gradient descent dijelaskan di sini: people.clarkson.edu/~ebollt/Papers/quadcost.pdf

dranxo
sumber
luar biasa .. terima kasih banyak! Saya akan mencoba kode ini, dan membandingkannya dengan saya untuk memeriksa kesalahan saya. Pendekatan ini sebenarnya tampaknya merupakan versi lokal dari makalah ini oleh Haker et al. yang saya sebutkan - yaitu., tanpa memecahkan laplacian. Terima kasih lagi !
WhitAngl
Saya akhirnya menemukan beberapa masalah dengan kode ini ...: jika saya menghitung Saya cukup jauh dari (dengan the hessian) - bahkan ketika melepas gaussian mengaburkan. Juga, jika saya hanya menambah jumlah iterasi ke beberapa ribu, kode meledak dan memberikan nilai NaN (dan crash). Ada ide ? Terima kasih! f2(ϕ)detHϕf1H
WhitAngl
Apakah lebih banyak blur membantu masalah NaN?
dranxo
Ya, memang, setelah lebih banyak tes itu membantu dengan lebih kabur - terima kasih !. Saya sekarang mencoba 8 langkah dengan blur awal 140 piksel standar deviasi hingga 1 piksel stdev (masih komputasi). Saya masih memiliki sejumlah besar gambar asli dalam hasil terakhir saya (dengan 64px blur). Saya juga akan memeriksa semua ikal yang tersisa di . ϕ
WhitAngl
Oh, bagus. Kupikir. Keburaman ada di sana karena gambar secara alami tidak kontinu (tepi) dan gradien akan bermasalah. Semoga Anda masih mendapatkan jawaban yang bagus tanpa mengaburkan terlalu banyak.
dranxo