[PYTHON] Implémentation de Scale-Space pour SIFT

Objectif de cet article

Vous serez en mesure de comprendre et de mettre en œuvre l'Echelle-espace qui prend en charge SIFT (Scale Invariant Feature Transform), qui décrit les caractéristiques des images qui résistent aux changements d'échelle et à la rotation.

Connaissances connexes

• Lissage d'image avec convolution gaussienne
• Transformation de Fourier et théorème d'échantillonnage

Environnement à utiliser

Remarque

Cet article a été rédigé en tant que matériel pour la 41e réunion de Morning Project Samurai.

S'il vous plaît

Ceci est une version provisoire. Si vous constatez des erreurs, nous vous serions reconnaissants de bien vouloir nous contacter.

Scale-space

Définition

• Générez plusieurs images avec différentes échelles à partir d'une image donnée.
• L'espace composé de ces images est appelé Scale-space.

motivation

Considérez la détection d'un objet particulier à partir d'une image donnée. Par exemple, envisagez de détecter une personne à partir d'une scène (image) d'une vidéo d'une route urbaine. À ce stade, il n'est pas nécessaire d'avoir une image à l'échelle qui montre les détails du motif des vêtements qu'une personne porte. Il est préférable d'avoir une image à une échelle qui vous permette de comprendre la forme d'une personne mais pas les vêtements que vous portez. Les informations trop détaillées sont plutôt bruyantes.

Comme mentionné ci-dessus, il existe une échelle d'image appropriée pour détecter l'objet d'intérêt. Cependant, en réalité, une image donnée n'est pas toujours à la bonne échelle. Par conséquent, un espace d'échelle est construit en générant une pluralité d'images ayant des échelles différentes à partir d'une image donnée, et une image d'une échelle appropriée est trouvée à partir de l'espace d'échelle pour détecter un objet.

Échelle et mise à l'échelle

Dans Scale-space, la mise à l'échelle est la suppression d'informations détaillées de l'image d'origine à l'aide d'une technique appelée convolution gaussienne. Le degré de mise à l'échelle est déterminé par le paramètre de convolution gaussienne $ \ sigma $. Plus l'échelle $ \ sigma $ est grande, plus les informations sont supprimées de l'image d'origine.

Par exemple, dans la figure ci-dessous, l'image originale à l'extrême gauche permet de voir clairement les détails du visage et même les poils en fourrure du chapeau. Dans l'image à l'extrême droite de l'échelle $ \ sigma = 3,2 $, vous pouvez voir qu'il y a une femme vague et les yeux sont clairs, mais rien de plus.

Comment configurer

1. Préparez l'image originale (échelle de gris).
2. Mettez à l'échelle $ \ sigma = k ^ {n} \ sigma_ {0} $ et définissez la valeur initiale de $ n $ sur $ 0 $.
3. Agrandissez l'image d'origine en utilisant la convolution gaussienne.
4. Répétez l'étape 3 en augmentant $ n $ à $ 1, 2, 3, ... $ jusqu'à $ k ^ {n} = 2,0 $. $ k $ est déterminé par le nombre d'images que vous souhaitez créer avant $ k ^ {n} = 2,0 $.
5. Sous-échantillonner l'image d'échelle $ \ sigma = 2.0 \ sigma_ {0} $ pour qu'elle fasse la moitié de la taille verticalement et horizontalement.
6. En utilisant l'image sous-échantillonnée comme image d'origine, définissez $ n = 1 $ et revenez à l'étape 3.

La méthode de construction détaillée et son programme seront décrits après avoir appris la convolution gaussienne, la transformation de Fourier, le théorème d'échantillonnage, etc.

Gaussian convolution

• Scale-space est utilisé pour générer une image d'échelle $ \ sigma $ à partir de l'image originale.
• Utilisé pour la suppression du bruit d'image et le flou dans le contexte général du traitement d'image.
• Mathématiquement, le noyau gaussien est utilisé pour trouver une série de moyennes pondérées pour n'importe quelle fonction.

Convolution gaussienne 1D

La convolution gaussienne unidimensionnelle est exprimée par l'équation suivante.

F(x, \sigma)=\int_{-\infty}^{\infty} f(u) \frac{1}{\sqrt{2 \pi \sigma^{2}}}e^{-\frac{(u-x)^2}{2\sigma^{2}}} du

$ f $ est le signal d'origine. $ \ frac {1} {\ sqrt {2 \ pi \ sigma ^ {2}}} e ^ {- \ frac {(ux) ^ 2} {2 \ sigma ^ {2}}} $ est appelé noyau gaussien .. $ \ Sigma $ est un paramètre qui détermine la forme du noyau gaussien. Le noyau gaussien a la forme d'une cloche avec une valeur maximale au centre $ x $, et sa hauteur et son étendue sont déterminées par $ \ sigma $.

Gaussian kernel Essayons de dessiner un noyau gaussien en utilisant numpy, scipy et matplotlib.

import numpy as np
from scipy.signal import gaussian
from matplotlib import pyplot as plt


if __name__ == '__main__':
    npoints = 101
    sigmas = np.array([6.0, 12.0])
    for i, sigma in enumerate(sigmas):
        y = gaussian(npoints, sigma) / (np.sqrt(2.0 * np.pi) * sigma)
        plt.subplot(len(sigmas), 1, i + 1)
        plt.title('sigma = %s' % sigma)
        plt.ylim(ymax=0.08)
        plt.plot(np.arange(npoints/2 - npoints, npoints/2, dtype=np.int), y)
    plt.tight_layout()
    plt.show()

Si vous exécutez le programme ci-dessus, vous obtiendrez le résultat ci-dessous.

[Ajout de la nature du noyau gaussien plus tard]

Exemple d'interprétation et d'application

La formule de convolution gaussienne est "le signal original $ f $ multiplié par le poids du centre $ x $ spread $ \ sigma $ (noyau gaussien) et la somme (moyenne pondérée) de $ F (x, \ sigma) $. Il peut être interprété comme "représenté par". A partir de là, $ F $ a l'information de la totalité de $ f $ tout en reflétant fortement l'information au voisinage de $ u = x $ du signal original $ f (u) $ avec la force déterminée par $ \ sigma $. Il peut être interprété comme un "nouveau signal constitué de $ F (x, \ sigma) ". Cela signifie que " F $ est un signal obtenu par lissage (flou, suppression du bruit) du signal original $ f $ avec une intensité de $ \ sigma ". Dans le contexte de Scale-space, on peut dire que " F $ est un signal obtenu en mettant à l'échelle $ f $ à l'échelle $ \ sigma $ (sans informations détaillées)".

Comme nous l'avons vu dans la section précédente, le noyau gaussien devient plus plat à mesure que $ \ sigma $ augmente. Cela signifie que "plus $ \ sigma $ est grand, plus la couleur de l'information est claire au voisinage de $ u = x $ du signal original $ f $ contenu dans $ F (x, \ sigma) $". Au fur et à mesure que sigma $ devient de plus en plus gros, $ F $ finira par avoir la même valeur pour tous les $ x $. Dans le contexte de Scale-space, "plus l'échelle $ \ sigma $ est grande, plus le signal est visualisé macroscopiquement (les informations détaillées sont perdues) et finalement toutes les informations contenues dans le signal sont identifiées." Peut être interprété comme.

Visualisons l'histoire jusqu'à présent avec un programme.

import numpy as np
from scipy.ndimage.filters import gaussian_filter1d
from matplotlib import pyplot as plt


if __name__ == '__main__':
    x = np.arange(0, 100)
    f = np.sin(0.5 * x) + np.random.normal(0, 6.0, x.shape)

    sigmas = [1.6, 3.2, 6.4]

    plt.subplot(len(sigmas) + 1, 1, 1)
    plt.ylim(ymin=np.min(f), ymax=np.max(f))
    plt.title('Original')
    plt.plot(f)

    for i, sigma in enumerate(sigmas):
        plt.subplot(len(sigmas) + 1, 1, 2 + i)
        plt.title('Sigma = %s' % sigma)
        plt.plot(gaussian_filter1d(f, sigma))
    plt.tight_layout()
    plt.show()

Lorsque le programme ci-dessus est exécuté, la figure suivante est obtenue.

Au fur et à mesure que $ \ sigma $ grandit (augmente), vous pouvez voir que les informations détaillées contenues dans le signal $ f $ sont perdues de $ F $.

Convolution gaussienne 2D

L'image est un signal bidimensionnel. Par conséquent, afin de mettre à l'échelle une image en utilisant la convolution gaussienne, une convolution gaussienne pour un signal bidimensionnel est nécessaire. Par la suite, lorsque nous écrivons simplement la convolution gaussienne, nous nous référerons à la convolution gaussienne bidimensionnelle.

La convolution gaussienne pour un signal bidimensionnel est exprimée par l'équation suivante.

F(x, y, \sigma) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(u, v) \frac{1}{2\sigma^{2}} e^{-\frac{(u - x)^2 + (v-y)^{2}}{\sigma^{2}}} dudv

$ \ frac {1} {2 \ sigma ^ {2}} e ^ {- \ frac {(u-x) ^ 2 + (vy) ^ {2}} {\ sigma ^ {2}}} $ est gaussien C'est un noyau. Ce noyau gaussien est un type de cloche tridimensionnel, avec $ (x, y) $ prenant la valeur maximale et $ \ sigma $ représentant sa propagation.

Gaussian kernel Dessinons le noyau gaussien par programmation.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def gaussian_kernel_2d(x, y, sigma):
    return np.exp(-(np.power(x/sigma, 2) + np.power(y/sigma, 2)))/(2 * np.power(sigma, 2))


if __name__ == '__main__':
    xrange = np.arange(-10.0, 10.5, 0.5)
    yrange = np.arange(-10.0, 10.5, 0.5)
    kernel_values = np.zeros(shape=(len(yrange), len(xrange)))

    sigmas = [1.6, 2.4]

    fig = plt.figure()
    X, Y = np.meshgrid(xrange, yrange)

    for i, sigma in enumerate(sigmas):
        ax = fig.add_subplot(len(sigmas), 1, 1 + i,  projection='3d')
        plt.title('Sigma = %s' % sigma)
        ax.set_zlim(zmax=0.2)
        ax.plot_wireframe(X, Y, gaussian_kernel_2d(X, Y, sigma))
    plt.show()

La figure suivante peut être obtenue en exécutant le programme ci-dessus. Il est possible de changer le point de vue avec la souris.

Convolution gaussienne discrète

Gaussian filter Les ordinateurs ne peuvent gérer que des valeurs discrètes et leur plage est finie. Le noyau gaussien est une fonction continue avec une zone de définition de $ \ pm \ infty $. Par conséquent, il est nécessaire d'approcher le noyau gaussien avec un nombre fini de chaînes numériques et de le disperser pour qu'il puisse être manipulé par un ordinateur.

Prenons une grille de $ s $ lignes et $ t $ colonnes. Par la suite, cette grille est appelée le filtre de $ s \ times t $. Ceci est représenté par programme dans le tableau à deux dimensions suivant.

filter = np.zeros(shape=(s, t))

Si vous pouvez définir tous les éléments du filtre variable pour représenter le noyau gaussien lorsque $ u = 0 $, vous pouvez créer un filtre gaussien $ g $ qui est une version discrète du noyau gaussien avec $ u = 0 $. Ici, l'élément $ k $ row $ l $ column de $ g $ $ g (k, l) $ = $ \ frac {1} {\ alpha} e ^ {- \ frac {l ^ {2} + k Définir comme ^ {2}} {\ sigma ^ {2}}} $. Ici, $ \ alpha $ est une constante qui se normalise de sorte que lorsque tous les éléments de $ g $ sont ajoutés, la somme est de 1 $.

Convolution gaussienne utilisant un filtre gaussien

L'équation suivante est une convolution gaussienne utilisant un filtre gaussien. C'est ce qu'on appelle la convolution gaussienne discrète. Par la suite, lorsque nous nous référons simplement à la convolution gaussienne, nous nous référerons à cette convolution gaussienne discrète.

F(i, j) = \sum_{k=0}^{s}\sum_{l=0}^{t} f(i+k-[\frac{s}{2}], j+l-[\frac{t}{2}]) g(k, l)

• Image originale de $ f $: $ m \ fois n $
• $ g $: $ s \ times t $ Filtre gaussien
• Image de sortie de $ F $: $ m \ times n $

Si vous créez un programme docilement selon la formule ci-dessus, ce sera comme suit.

for k in range(s):
    for l in range(t):
        F[i, j] = f[i + k - int(s/2), j + l - int(t/2)] * g(k, l)

Dans ce cas, $ s \ fois t $ opérations sont nécessaires pour obtenir un pixel de l'image de sortie $ F $, et $ m \ fois n \ fois s \ fois t $ pour obtenir tous les pixels de l'image de sortie. Cela nécessite un certain nombre d'opérations. La convolution avec une image à 320 $ \ fois 240 $ et un filtre gaussien à 5 $ \ fois 5 $ nécessite un total de 1 920 000 $ d'opérations. Le fait que ce nombre de calculs soit critique ou non dépend des spécifications de la machine et de l'application, mais c'est essentiel sur mon MacBook Air (13 pouces, début 2015).

Accélérant

La convolution gaussienne a la propriété que l'image de sortie $ F $ de convolution gaussienne utilisant le filtre gaussien de $ s \ fois t $ est représentée par l'image $ F_ {1} $ sortie dans les deux étapes suivantes. Il y a.

1. Convolution du filtre gaussien unidimensionnel de $ 1 \ fois t $ avec chaque ligne de l'image originale $ f $ et sortie de l'image $ F_ {0} $.
2. Convolution du filtre gaussien unidimensionnel de $ s \ fois 1 $ avec chaque colonne de $ F_ {0} $ et sortie de l'image $ F_ {1} $.

Ceci est prouvé comme suit.

\begin{align}
F(i, j) &=\sum_{k=0}^{s}\sum_{l=0}^{t} f(i+k-[\frac{s}{2}], j+l-[\frac{t}{2}]) g(k, l)\\

&= \sum_{k=0}^{s}\sum_{l=0}^{t} f(i+k-[\frac{s}{2}], j+l-[\frac{t}{2}]) \frac{1}{\alpha}e^{-\frac{l^{2} + k^{2}}{\sigma^{2}}}\\

&= \sum_{k=0}^{s}\sum_{l=0}^{t} f(i+k-[\frac{s}{2}], j+l-[\frac{t}{2}]) \frac{1}{\alpha}e^{-\frac{l^{2}}{\sigma^{2}}}e^{-\frac{k^{2}}{\sigma^{2}}}\\

&=\frac{\alpha_{0}\alpha_{1}}{\alpha} \sum_{k=0}^{s}\frac{1}{\alpha_{1}}e^{\frac{-k^{2}}{\sigma^{2}}}\sum_{l=0}^{t} f(i+k-[\frac{s}{2}], j+l-[\frac{t}{2}]) \frac{1}{\alpha_{0}}e^{-\frac{l^{2}}{\sigma^{2}}}\\

&=\frac{\alpha_{0}\alpha_{1}}{\alpha} \sum_{k=0}^{s}\frac{1}{\alpha_{1}}e^{\frac{-k^{2}}{\sigma^{2}}}\sum_{l=0}^{t} f(i+k-[\frac{s}{2}], j+l-[\frac{t}{2}]) g_{0}(l)\\


&=\frac{\alpha_{0}\alpha_{1}}{\alpha} \sum_{k=0}^{s}\frac{1}{\alpha_{1}}e^{-\frac{k^{2}}{\sigma^{2}}} F_{0}(i, j)\\


&=\frac{\alpha_{0}\alpha_{1}}{\alpha} \sum_{k=0}^{s} F_{0}(i, j) g_{1}(k)\\


&=\frac{\alpha_{0}\alpha_{1}}{\alpha} F_{1}(i, j)
\end{align}

• $ g_ {0} $: $ 1 \ times t $ Filtre gaussien
• $ g_ {1} $: $ s \ fois 1 $ Filtre gaussien
• $ \ alpha $: Constante de normalisation de $ g $
• $ \ alpha_ {0} $: Constante de normalisation de $ g_ {0} $
• $ \ alpha $: Constante de normalisation de $ g_ {1} $

Ce faisant, le nombre d'opérations est de $ s \ fois m + t \ fois n $. Pour une image à 320 $ \ fois 240 $ et une convolution de filtre gaussien à 5 $ \ fois 5 $, le nombre d'opérations est de 768 000 $. En accélérant, le nombre d'opérations peut être réduit de 60 $ % $.

Appliquer à l'image

Appliquons en fait la convolution gaussienne à l'image. Cette fois, Pillow, Numpy et Scipy sont utilisés. Le programme suivant génère une image d'échelle $ \ sigma = 1,6 $ et une image d'échelle $ \ sigma = 3,2 $.

from PIL import Image
import numpy as np
from scipy.ndimage.filters import gaussian_filter
from matplotlib import pyplot as plt

if __name__ == '__main__':
    orig_image = Image.open('lena.jpg').convert('L')
    orig_image = np.array(orig_image, dtype=np.uint8)

    sigmas = [1.6, 3.2]

    plt.subplot(1, len(sigmas) + 1, 1)
    plt.title('Orig')
    plt.imshow(orig_image, cmap='Greys_r')
    for i, sigma in enumerate(sigmas):
        plt.subplot(1, len(sigmas) + 1, 2 + i)
        plt.title('Sigma=%s' % sigma)
        plt.imshow(gaussian_filter(orig_image, sigma), cmap='Greys_r')
    plt.tight_layout()
    plt.savefig('gausian_convolution_2d_examples.png')
    plt.show()

Lorsque le programme ci-dessus est exécuté, la figure suivante est obtenue.

OpenCV est une bibliothèque bien connue pour travailler avec des images. Si vous l'avez installé, vous pouvez l'utiliser. Vous pouvez également essayer d'implémenter vous-même la convolution gaussienne.

À ce stade, vous pouvez générer une image de n'importe quelle échelle $ \ sigma $ par rapport à l'image d'origine.

Bâtiment échelle-espace 1ère octave

Enfin, nous construirons Scale-space. La 1ère octave du titre du chapitre sera décrite plus tard.

Flux de construction

1. Déterminez l'échelle de base $ \ sigma_ {0} $.
2. Générez une image avec l'échelle $ k ^ {n} \ sigma $ pour $ n = 0, ..., s $. $ s $ est le nombre de divisions de la 1ère octave de l'espace d'échelle. Il est représenté par $ k = 2 ^ {1 / s} $.

Programme pour construire Scale-space 1ère octave

from PIL import Image
import numpy as np
from scipy.ndimage.filters import gaussian_filter
from matplotlib import pyplot as plt


if __name__ == '__main__':
    orig_image = Image.open('lena.jpg').convert('L')
    orig_image = np.array(orig_image, dtype=np.uint8)

    sigma = 1.6
    s = 3
    k = np.power(2, 1/s)
    scale_space = []

    for n in range(s + 1):
        scale_space.append(gaussian_filter(orig_image, np.power(k, n) * sigma))

    for n, img in enumerate(scale_space):
        plt.subplot(1, len(scale_space), 1 + n)
        plt.title('Sigma=%s' % np.round(np.power(k, n) * sigma, 2))
        plt.imshow(img, cmap='Greys_r')
    plt.tight_layout()
    plt.show()

Lorsque le programme ci-dessus est exécuté, quatre images avec des échelles de 1,6 $ à 3,2 $ sont stockées dans la variable scale_space, et la figure ci-dessous est affichée.

L'espace constitué des images de $ k \ sigma $ à $ 2k \ sigma $ est appelé 1ère octave.

Transformation de Fourier et théorème d'échantillonnage

Lors de la création de la 2e octave et des suivantes, il est nécessaire de mentionner la transformation de Fourier et le théorème d'échantillonnage. Cependant, cette fois, je ne vais pas aller plus loin et simplement laisser couler doucement.

Conversion de Fourier

La transformation de Fourier transforme une image d'une région de pixels en une région de fréquence.

Formule de transformation de Fourier 1D

F(j\omega) = \int_{-\infty}^{\infty} f(x) e^{-j\omega x} dx

• $ f $: Signal de la zone de localisation (par exemple, luminosité de la lumière passant à travers une fente allongée)
• $ F $: Représentation de $ f $ dans la plage de fréquences
• $ \ omega : Fréquence carrée ( 2 \ pi \ omega_ {0} $, où $ \ omega_ {0} $ est la fréquence)
• $e^{-j\omega x} = \cos(\omega x) - j \sin(\omega x) $
• $ j $: Unité imaginaire

Exemple d'application de la transformation de Fourier unidimensionnelle

Essayons la transformation de Fourier par FFT (Fast Fourier Transform) en utilisant numpy, scipy, matplotlib. Dans le programme ci-dessous, la fréquence d'échantillonnage est de 20 et la fréquence angulaire du signal d'origine est de 2 $ \ pi $.

import numpy as np
from scipy.fftpack import fft
from matplotlib import pyplot as plt


if __name__ == '__main__':
    sampling_rate = 20
    sampling_interval = 1.0/sampling_rate
    x = np.arange(0, 1, sampling_interval)

    omega0 = 1.0
    omega = 2.0 * np.pi * omega0
    f = np.sin(omega * x)

    F = fft(f)

    plt.plot(np.arange(-len(f)/2, len(f)/2), np.abs(np.concatenate((F[len(f)/2:], F[:len(f)/2]))), 'bo-')
    plt.show()

Lorsque le programme ci-dessus est exécuté, la figure suivante est obtenue.

A partir de cette figure, on peut voir que le signal original $ f $ contient un signal avec une fréquence de $ 1 $. C'est une bonne idée d'effectuer une conversion de Fourier sur la superposition de plusieurs fonctions périodiques.

[J'écrirai un peu plus en détail plus tard]

Théorème d'échantillonnage de signal unidimensionnel

Le théorème d'échantillonnage déclare: «Lorsque la période maximale du signal contenu dans le signal original $ f $ est $ W $, si la période d'échantillonnage $ T $ du signal original est $ T \ leq \ frac {1} {W} $ , Le signal original $ f $ est complètement récupérable à partir du signal échantillonné $ f_ {d} $. "

[Une preuve sera ajoutée à une date ultérieure]

Expérience avec un théorème d'échantillonnage pour un signal unidimensionnel

Soit la période d'échantillonnage de 1/20 $. Ensuite, selon le théorème d'échantillonnage, si le signal a une fréquence maximale de 10 $, toutes les informations de fréquence doivent être extraites avec succès par conversion de Fourier.

Vérifions avec le programme ci-dessous.

import numpy as np
from scipy.fftpack import fft
from matplotlib import pyplot as plt

if __name__ == '__main__':
    sampling_rate = 20
    sampling_interval = 1.0/sampling_rate
    x = np.arange(0, 1, sampling_interval)

    omega0s = [9, 10, 11, 30]

    for i, omega0 in enumerate(omega0s):
        omega = 2.0 * np.pi * omega0
        f = np.sin(omega * x)

        F = fft(f)

        plt.subplot(len(omega0s),  1,  1 + i)
        plt.plot(np.arange(-len(f)/2, len(f)/2), np.abs(np.concatenate((F[len(f)/2:], F[:len(f)/2]))), 'bo-')
        plt.title('Max frequency = %s' % omega0)
        plt.xlabel('Frequency')
        plt.ylabel('Amplitude spectrum')
    plt.tight_layout()
    plt.show()

Une fois exécuté, la figure suivante est obtenue. À partir de ce résultat, on peut voir que les informations de fréquence peuvent être extraites avec succès avec le signal ayant la fréquence maximale de 9 $ (les détails seront ajoutés plus tard).

Transformation de Fourier bidimensionnelle

[Postscript]

Construire après la 2e octave

Sous-échantillonnage

Aperçu

[Postscript]

Preuve de sous-échantillonnage

[Postscript]

Programme de construction de 2e octave

[Postscript]

Les références

• Distinctive Image Features from Scale-Invariant Keypoints (David G. Lowe)
• Traitement du signal (Iwao Morishita, Hidefumi Obata)