Implementação ideal de distorção de transporte no Matlab

11

Estou implementando o artigo " Transporte de massa ideal para registro e entortamento ", meu objetivo é colocá-lo online, pois não consigo encontrar nenhum código de transporte de massa euleriano on-line e isso seria interessante pelo menos para a comunidade de pesquisa em processamento de imagens.

O artigo pode ser resumido da seguinte forma:
- encontre um mapa inicial usando as correspondências do histograma 1D ao longo das coordenadas xey - resolva o ponto fixo de , em que representa uma rotação de 90 graus no sentido anti-horário, para a solução da equação de poisson com condições de contorno de Dirichlet (= 0), e é o determinante da matriz jacobiana. - a estabilidade é garantida para um intervalo de tempou
ut=1μ0Du1div(u)u1Du
dt<min|1μ01div(u)|

Para simulações numéricas (realizadas em uma grade regular), elas indicam o uso de poicalc do matlab para resolver a equação de poisson, usam diferenças finitas centralizadas para derivadas espaciais, exceto Du que é calculado usando um esquema a favor do vento.

Usando meu código, a energia funcional e a curvatura do mapeamento estão diminuindo adequadamente em algumas iterações (de algumas dezenas a alguns milhares, dependendo da etapa do tempo). Mas depois disso, a simulação explode: a energia aumenta para alcançar um NAN em muito poucas iterações. Tentei várias ordens para diferenciações e integrações (uma substituição de ordem superior ao cumptrapz pode ser encontrada aqui ) e esquemas de interpolação diferentes, mas sempre recebo o mesmo problema (mesmo em imagens muito suaves, diferentes de zero em todos os lugares etc.).
Alguém estaria interessado em olhar para o código e / ou o problema teórico que estou enfrentando? O código é bastante curto.

Substitua gradiente2 () no final por gradiente (). Este era um gradiente de ordem mais alta, mas também não resolve as coisas.

Por enquanto, estou interessado apenas na parte ótima do transporte, não no termo adicional de regularização.

Obrigado !

WhitAngl
fonte

Respostas:

5

Meu bom amigo Pascal fez isso há alguns anos (está quase no 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)

A execução do teste leva cerca de 2 segundos.

A abordagem de descida em gradiente é explicada aqui: people.clarkson.edu/~ebollt/Papers/quadcost.pdf

dranxo
fonte
excelente .. muito obrigado! Vou tentar esse código e comparar com o meu para verificar meus erros. Essa abordagem, na verdade, parece ser a versão local do artigo de Haker et al. que eu mencionei - isto é, sem resolver para um laplaciano. Obrigado novamente !
WhitAngl
Finalmente estou encontrando alguns problemas com esse código ...: se eu calcular , estou bem longe de (com o hessian) - mesmo ao remover o gaussian borrão. Além disso, se eu apenas aumentar o número de iterações para alguns milhares, o código explode e fornece valores de NaN (e falha). Qualquer ideia ? Obrigado ! f2(ϕ)detHϕf1H
WhitAngl
Mais desfoque ajuda no problema do NaN?
Dranxo 22/05/12
Sim, de fato, depois de mais testes, ajuda com mais desfoque - obrigado !. Agora estou tentando 8 etapas com um desfoque inicial de desvio padrão de 140 pixels até 1 pixel stdev (ainda computando). Ainda tenho uma quantidade significativa da imagem original no meu último resultado (com desfoque de 64px). Também vou verificar se há alguma curvatura restante no . ϕ
WhitAngl
Ah, tudo bem. Eu acho que. O desfoque existe porque as imagens naturalmente não são contínuas (bordas) e o gradiente será problemático. Espero que você ainda esteja obtendo boas respostas sem desfocar muito.
Dranxo 22/05