Implementación de deformación de transporte óptima en Matlab

11

Estoy implementando el documento " Transporte masivo óptimo para el registro y la deformación ", mi objetivo es ponerlo en línea ya que simplemente no puedo encontrar ningún código de transporte masivo euleriano en línea y esto sería interesante al menos para la comunidad de investigadores en el procesamiento de imágenes.

El documento se puede resumir de la siguiente manera:
- encuentre un mapa inicial utilizando las coincidencias de histograma 1D a lo largo de las coordenadas x e y - resuelva el punto fijo de u t = 1u
ut=1μ0Du1div(u)u1Du
dt<min|1μ01div(u)|

Para las simulaciones numéricas (realizadas en una cuadrícula regular), indican el uso de poicalc de matlab para resolver la ecuación de poisson, usan diferencias finitas centradas para derivadas espaciales, excepto que se calcula utilizando un esquema de viento ascendente.Du

Usando mi código, la energía funcional y la curvatura del mapeo disminuyen adecuadamente durante un par de iteraciones (de unas pocas decenas a unos miles, dependiendo del paso de tiempo). Pero después de eso, la simulación explota: la energía aumenta para alcanzar una NAN en muy pocas iteraciones. Intenté varias órdenes para las diferenciaciones e integraciones ( aquí se puede encontrar un reemplazo de orden superior a cumptrapz ), y diferentes esquemas de interpolación, pero siempre obtengo el mismo problema (incluso en imágenes muy suaves, distintas de cero en todas partes, etc.).
¿Alguien estaría interesado en mirar el código y / o el problema teórico que estoy enfrentando? El código es bastante corto.

Reemplace gradiente2 () al final por gradiente (). Este era un gradiente de orden superior pero tampoco resuelve las cosas.

Por ahora solo me interesa la parte de transporte óptima del documento, no el término de regularización adicional.

Gracias !

WhitAngl
fuente

Respuestas:

5

Mi buen amigo Pascal hizo esto hace unos años (está casi en 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)

Prueba de funcionamiento, tarda unos 2 segundos.

El enfoque de descenso de gradiente se explica aquí: people.clarkson.edu/~ebollt/Papers/quadcost.pdf

dranxo
fuente
excelente .. muchas gracias! Probaré este código y lo compararé con el mío para verificar mis errores. Este enfoque en realidad parece ser la versión local del documento de Haker et al. que mencioné, es decir, sin resolver por un laplaciano. Gracias de nuevo !
WhitAngl
Finalmente me encuentro con un par de problemas con este código ...: si calculo Estoy bastante lejos de (con el hessian), incluso cuando elimino el gaussiano difuminar. Además, si solo aumento el número de iteraciones a un par de miles, el código explota y da valores NaN (y se bloquea). Alguna idea ? Gracias ! f2(ϕ)detHϕf1H
WhitAngl
¿Más desenfoque ayuda con el problema de NaN?
dranxo
Sí, de hecho, después de más pruebas, ayuda con más desenfoque, ¡gracias! Ahora estoy intentando 8 pasos con un desenfoque inicial de 140 píxeles de desviación estándar hasta 1 píxel stdev (todavía computación). Sin embargo, todavía tengo una cantidad significativa de la imagen original en mi último resultado (con un desenfoque de 64 píxeles). También comprobaré si hay rizos restantes en . ϕ
WhitAngl
Oh ok bien. Yo creo que. El desenfoque está ahí ya que las imágenes, naturalmente, no son continuas (bordes) y el gradiente será problemático. Espero que sigas obteniendo buenas respuestas sin desdibujar demasiado.
dranxo