Sub-pixel shift via Fourier geeft rare resultaten

[317070]

Ik probeer via het Fourier-domein een functie op andere plaatsen te sampelen, bijvoorbeeld om te zien wat er gebeurt met mijn afbeelding als ik het een halve pixel opschuif. Dat doe ik met de translatie-eigenschap in het Fourier domein: https://en.wikipedia.org/wiki/Fourier_transform#Basic_properties

Zie bijvoorbeeld de publicaties:
http://www.mathworks.com/matlabcentral/fileexchange/18401-efficient-subpixel-image-registration-by-cross-correlation/content/html/efficient_subpixel_registration.html
Het einddoel: http://spin.ecn.purdue.edu/fmri/PDFLibrary/TongR_MRM_1999_41_253_256.pdf
 
De resultaten zijn daar steeds erg bevredigend, maar bij mij loopt er nog iets mis. Als ik bijvoorbeeld de afbeelding van mijn koe een halve pixel opschuif, krijg ik het volgende:
Origineel:

Het resultaat:

 
 
Ik plaats er nog even een minimal working code example in python bij:

Code: Selecteer alles

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

N = 2048

img=mpimg.imread('koe.png') #http://www.clipartbest.com/cliparts/RiA/yXL/RiAyXLriL.png

f = np.zeros((N,N))
f[:img.shape[0],:img.shape[1]] = img[:,:,3]

def fft(x):
    return np.fft.fft2(x)

def ifft(x):
    return np.fft.ifft2(x)

exp = np.exp
pi = np.pi
j = 1j

def s(r):
    plt.figure()
    #plt.imshow(np.real(r), vmin=0.0, vmax=1.0, interpolation='none')
    plt.imshow(np.absolute(r), interpolation='none', vmin=0.0, vmax=1.0)
    plt.colorbar()

x = np.tile(np.arange(N),(N,1)).T
y = x.T
result = ifft( fft(f)* np.exp(1j*2*np.pi*(-0.5*x/N-0.5*y/N)) )
s(result)
s(f)
plt.show()

Weet er iemand raad? Hoe kan ik wel een redelijk resultaat bekomen?

[Xenion]

Je afbeelding is niet vierkant, waarom definieer je dan een NxN grid?
Dan verder: de x en y grids die je definieert gaan van 0 tot het aantal rijen/kolommen. Als je het MATLAB voorbeeld volgt dan zie je dat die maar opbouwen tot de helft en dan negatief worden en weer opbouwen tot -1. Dit om rekening te houden met de volgorde waarin de fft functie het spectrum weergeeft.
 
Volgende code werkt voor mij wel:

Code: Selecteer alles

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

img=mpimg.imread('D:\koe.png') #http://www.clipartbest.com/cliparts/RiA/yXL/RiAyXLriL.png

nr = len(img)
nc = len(img[0])

f = np.zeros((nr,nc))
f[:img.shape[0],:img.shape[1]] = img[:,:,3]

def fft(x):
    return np.fft.fft2(x)

def ifft(x):
    return np.fft.ifft2(x)

exp = np.exp
pi = np.pi
j = 1j

def s(r):
    plt.figure()
    #plt.imshow(np.real(r), vmin=0.0, vmax=1.0, interpolation='none')
    plt.imshow(np.absolute(r), interpolation='none', vmin=0.0, vmax=1.0)
    plt.colorbar()
    
Nc = np.tile(np.arange(nr),(nc,1)).T
Nc[nr/2:] -= nr
Nr = np.tile(np.arange(nc),(nr,1)).T
Nr[nc/2:] -= nc
Nr = Nr.T
result = ifft( fft(f)* np.exp(1j*2*np.pi*(-0.5*Nc/nc-0.5*Nc/nr)) )
s(f)
s(result)
plt.show()

[317070]

Je afbeelding is niet vierkant, waarom definieer je dan een NxN grid?
Dan verder: de x en y grids die je definieert gaan van 0 tot het aantal rijen/kolommen. Als je het MATLAB voorbeeld volgt dan zie je dat die maar opbouwen tot de helft en dan negatief worden en weer opbouwen tot -1. Dit om rekening te houden met de volgorde waarin de fft functie het spectrum weergeeft.

Dat vierkant is omdat de uiteindelijke implementatie toch op GPU komt, en die doet enkel power-2 fft's.
 
Je oplossing is inderdaad juist. Ongelofelijk bedankt! Maar nu zit ik met het volgende probleem bij rotatie via het Fourier-domein. Weet je waar zo'n effect vandaan kan komen? Het heeft waarschijnlijk ook iets te maken met het kiezen van de N'en, aangezien het voor die aanpassing wel goed werkte (op die interpolatie met te hoge frequenties na, zoals in #1) maar ik zie begot niet wat dit effect kan veroorzaken.
 

koe2 551 keer bekeken

Code: Selecteer alles

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import math

N = 2048

img=mpimg.imread('koe.png')

f = np.zeros((N,N))
f[:img.shape[0],:img.shape[1]] = img[:,:,3]
f = f[::8,::8]
N=N/8

def fft(x):
    return np.fft.fft2(x)

def ifft(x):
    return np.fft.ifft2(x)

exp = np.exp
pi = np.pi
j = 1j

def convertfreqs(x):
    return np.mod((x + N/2),N) - N/2

def s(r):
    plt.figure()
    #plt.imshow(np.real(r), vmin=0.0, vmax=1.0, interpolation='none')
    plt.imshow(np.absolute(r), interpolation='none', vmin=0.0, vmax=1.0)
    plt.colorbar()

x = np.tile(convertfreqs(np.arange(N)),(N,1)).T
y = x.T

freq = np.zeros((N,N),dtype='complex64')
S = 2
freq[:N/S:,:N/S:] = fft(f)[:N/S:,:N/S:]
freq[(S-1)*N/S::,(S-1)*N/S::] = fft(f)[(S-1)*N/S::,(S-1)*N/S::]
freq[(S-1)*N/S::,:N/S:] = fft(f)[(S-1)*N/S::,:N/S:]
freq[:N/S:,(S-1)*N/S::] = fft(f)[:N/S:,(S-1)*N/S::]

angle = pi/4
translation = (0.0,0.0) #( 0.5 * N , (-np.sqrt(2.)/2+0.5) * N)
a = np.cos(angle)
b = np.sin(angle)

x = np.tile(convertfreqs(np.arange(N)),(N,1)).T
y = x.T
z1  = exp(  pi*1j * ( (x**2 - y**2)*a + 2*x*y*b)/N)

#still a bug!

t = np.zeros((4*N,4*N),dtype='complex64')
t[:N:,:N:]=z1 * freq * 1./N**2 * exp( 2*pi*1j* (x*translation[0] + y*translation[1]) / N )

x = np.tile(convertfreqs(np.arange(N-1,-N-1,-1)),(2*N,1)).T
y = np.tile(convertfreqs(np.arange(2*N-1,-1,-1)),(2*N,1))

big_z2 = np.zeros((4*N,4*N),dtype='complex64')
big_z2[:2*N,:2*N] = exp(  pi*1j * ( (x**2 - y**2)*a + 2*x*y*b)/N)

c1 = ifft( fft(t) * fft(big_z2.conjugate()) )
c1 = c1[N:2*N, 2*N-1:N-1:-1]
c1 = c1 * z1.T.conjugate()
s(c1)
plt.show()

[Xenion]

Ben ni zeker maar ik vermoed dat je nu die freq matrix niet goed definieert. Kijk telkens goed na waar de negatieve frequenties zitten en zie dat je daar juist mee rekent. Die convertFreqs functie definieert het x y grid om daar rekening mee te houden maar mss dat je nu de matrix ineen puzzelt a la fftshift in MATLAB?

Ik kan er zelf morgen pas aan prutsen maar tegen dan heb je het zelf mss al gefixed. Plot ook mss eens de spectra en exp functies om te debuggen.

[317070]

Ik kan er zelf morgen pas aan prutsen maar tegen dan heb je het zelf mss al gefixed. Plot ook mss eens de spectra en exp functies om te debuggen.

Er ontbrak een fftshift bij het definiëren van de matrix t. Ik zit nog met off-by-one errors e.d. Maar de essentie is er nu  Nog eens bedankt voor je hulp!

[Xenion]

Cool Die off by one errors komen mogelijk van een kleine fout in het opsplitsen van het spectrum. Probeer eens de splits grens met 1 naar links of rechts op te schuiven.

Waar ben je eigenlijk mee bezig? Ik dacht dat jij met robots werkte?

Bij welk departement werk jij trouwens? Ik was vandaag op labo bezoek bij UGent TELIN-IPI.

[317070]

Waar ben je eigenlijk mee bezig? Ik dacht dat jij met robots werkte?

Bij welk departement werk jij trouwens? Ik was vandaag op labo bezoek bij UGent TELIN-IPI.

Ik ben bezig met de grens tussen machine learning en robots, en dit valt wat meer in de eerste categorie. En vandaar mijn rare benadering, maar mijn transformatie moet in alle mogelijke richting differentieerbaar zijn op GPU om de gradient er door te kunnen berekenen voor een gradient descent.
Vakgroep ELIS trouwens 
 
Moest je trouwens een (vaag) idee hebben waar de spookbeelden die ik nog heb (in het lichtblauw tussen de gezichten) vandaan kunnen komen, mag je altijd suggesties doen.   De spookbeelden zien eruit. de afbeelding die ik heb, maar verschoven over (N/2;N/2) en met enkel de randen zichtbaar.

result 553 keer bekeken

[Xenion]

Ik weet niet wàt precies het probleem is, want echt praktische ervaring heb ik hier niet mee maar op het eerste zicht denk ik aan aliasing-achtige dingen en omdat je je afbeelding zero-padd introduceer je mogelijk niet-bestaande hoge frequenties. Probeer eens een symmetric boundary extension ofzo.