Introducción al Algoritmo SIFT

La Transformación de Características Invariante a Escala, o SIFT (Scale-invariant Feature Transform), es un algoritmo fundamental en visión por computadora. Fue introducido por David Lowe en 1999 y refinado en 2004, diseñado para identificar y describir características locales dentro de una imagen. Su fortaleza radica en la detección de puntos extremos en el espacio de escala, extrayendo su posición, escala y orientación de manera invariante.

Las aplicaciones de SIFT son diversas, abarcando el reconocimiento de objetos, mapeo y navegación robótica, unión de imágenes (image stitching), construcción de modelos 3D, reconocimiento de gestos, seguimiento de objetos y comparación de movimientos.

La detección y descripción de características locales son cruciales para el reconocimiento de objetos. Las características SIFT se basan en puntos de interés local de un objeto, siendo robustas ante variaciones en tamaño, rotación, iluminación, ruido y leves cambios de perspectiva. Gracias a estas propiedades, son altamente distintivas y relativamente fáciles de obtener, lo que facilita el reconocimiento de objetos incluso en grandes bases de datos, con una baja tasa de falsos positivos. Los descriptores SIFT también son eficaces en la detección de objetos parcialmente ocluidos, requiriendo solo unos pocos puntos SIFT para estimar la posición y orientación del objeto. Con el hardware actual y bases de datos de características moderadas, la velocidad de reconocimiento se aproxima al tiempo real. La gran cantidad de información contenida en las características SIFT las hace ideales para un emparejamiento rápido y preciso en grandes volúmenes de datos.

En esencia, el algoritmo SIFT busca puntos clave (características) a través de diferentes escalas de una imagen y calcula su orientación. Estos puntos clave son ubicaciones prominentes que son estables ante cambios de iluminación, transformaciones afines y ruido, como esquinas, puntos en bordes, o puntos claros en regiones oscuras y oscuros en regiones claras.

Etapas del Algoritmo SIFT

El proceso de extracción de características SIFT se puede desglosar en las siguientes fases:

1. Construcción de Pirámides de Imágenes

1.1. Pirámide Gaussiana

La pirámide gaussiana de una imagen se construye aplicando una función gaussiana para difuminar la imagen y luego submuestrear la imagen resultante. Este proceso genera múltiples octavas (grupos) y niveles (capas) dentro de cada octava, donde cada nivel representa la imagen con un diferente grado de suavizado gaussiano.

La creación de la pirámide implica aplicar un filtro gaussiano con un coeficiente de suavizado σ específico. Un principio común es que el tamaño de la ventana del kernel (N) para el filtro gaussiano se determina por N = ⌈6σ+1⌉, tomando el número impar más cercano.

Convolución Gaussiana Separable

La convolución directa de la imagen con un kernel gaussiano 2D puede ser computacionalmente costosa y puede introducir artefactos en los bordes. Para optimizar este proceso, se utiliza la convolución separable. Esto implica aplicar un filtro 1xN a lo largo de la dirección X y luego un filtro Nx1 a lo largo de la dirección Y. Este método no solo ahorra tiempo de cómputo, sino que también minimiza la pérdida de información en los bordes.

Análisis de Código para la Pirámide Gaussiana

A continuación, se presenta un pseudocódigo que ilustra la construcción de la pirámide gaussiana:

// Pseudocódigo para construir la pirámide gaussiana
function construirPiramideGaussiana(imagenBase, numOctavas, numNivelesPorOctava, sigmas)
    piramide_gaussiana = array[numOctavas][numNivelesPorOctava + 3]

    for octava = 0 to numOctavas - 1:
        for nivel = 0 to numNivelesPorOctava + 2:
            if octava == 0 and nivel == 0:
                piramide_gaussiana[octava][nivel] = copiarImagen(imagenBase)
            else if nivel == 0: // Base de una nueva octava es la imagen submuestreada de la anterior
                imagen_anterior_octava = piramide_gaussiana[octava - 1][numNivelesPorOctava]
                piramide_gaussiana[octava][nivel] = submuestrear(imagen_anterior_octava)
            else: // Difuminar la imagen actual para crear la siguiente en la misma octava
                imagen_previa = piramide_gaussiana[octava][nivel - 1]
                nueva_imagen = crearImagen(tamano(imagen_previa), tipo_flotante_32, 1)
                aplicarSuavizadoGaussiano(imagen_previa, nueva_imagen, sigmas[nivel], sigmas[nivel])
                piramide_gaussiana[octava][nivel] = nueva_imagen
    return piramide_gaussiana

1.2. Pirámide de Diferencias de Gaussianas (DoG)

En 2002, Mikolajczyk demostró experimentalmente que los máximos y mínimos de la función Laplaciana de Gaussiana (LoG) normalizada en escala producen las características de imagen más estables en comparación con otros extractores. Lindeberg (1994) ya había descubierto que la función Diferencia de Gaussianas (DoG) es una excelente aproximación de la función LoG normalizada en escala, donde:

DoG(x, y, σ) ≈ (k-1) * σ^2 * LoG(x, y, σ)

Donde k-1 es una constante que no afecta la ubicación de los puntos extremos.

Construcción de la Pirámide DoG

La pirámide DoG se construye a partir de la pirámide gaussiana. En cada octava de la pirámide gaussiana, se restan imágenes adyacentes para generar las capas de la pirámide DoG (por ejemplo, la capa i+1 menos la capa i).

Análisis de Código para la Pirámide DoG

El siguiente pseudocódigo ilustra la generación de la pirámide DoG:

// Pseudocódigo para construir la pirámide DoG
function construirPiramideDoG(piramideGaussiana, numOctavas, numNivelesPorOctava)
    piramide_dog = array[numOctavas][numNivelesPorOctava + 2]

    for octava = 0 to numOctavas - 1:
        for nivel = 0 to numNivelesPorOctava + 1: // Necesitamos numNivelesPorOctava + 2 capas en DoG
            imagen_superior = piramideGaussiana[octava][nivel + 1]
            imagen_inferior = piramideGaussiana[octava][nivel]
            nueva_imagen_dog = crearImagen(tamano(imagen_superior), tipo_flotante_32, 1)
            restarImagenes(imagen_superior, imagen_inferior, nueva_imagen_dog)
            piramide_dog[octava][nivel] = nueva_imagen_dog
    return piramide_dog

2. Detección de Extremos en el Espacio de Escala (Puntos Clave)

Los puntos clave se identifican como puntos extremos locales en el espacio DoG. Para ello, cada píxel en la pirámide DoG se compara con sus 26 vecinos: los 8 vecinos en la misma capa, y los 9 vecinos (incluyendo el propio píxel) en las capas superior e inferior adyacentes. Un píxel es considerado un extremo si es mayor o menor que todos sus vecinos, tanto en el dominio espacial 2D como en el dominio de la escala.

Análisis de Código para la Detección de Extremos

El siguiente pseudocódigo demuestra la lógica para la detección de extremos:

// Pseudocódigo para la detección de extremos en la pirámide DoG
function esPuntoExtremo(piramideDoG, octava_idx, nivel_idx, fila_idx, col_idx):
    valor_central = obtenerValorPixel(piramideDoG[octava_idx][nivel_idx], fila_idx, col_idx)

    // Si el valor es positivo, buscar un máximo local
    if valor_central > 0:
        for i = -1 to 1:
            for j = -1 to 1:
                for k = -1 to 1:
                    // Excluir el propio punto central
                    if i == 0 and j == 0 and k == 0:
                        continue
                    
                    valor_vecino = obtenerValorPixel(piramideDoG[octava_idx][nivel_idx + i], fila_idx + j, col_idx + k)
                    if valor_central < valor_vecino:
                        return false // No es un máximo
        return true // Es un máximo
    else: // Si el valor es negativo o cero, buscar un mínimo local
        for i = -1 to 1:
            for j = -1 to 1:
                for k = -1 to 1:
                    // Excluir el propio punto central
                    if i == 0 and j == 0 and k == 0:
                        continue

                    valor_vecino = obtenerValorPixel(piramideDoG[octava_idx][nivel_idx + i], fila_idx + j, col_idx + k)
                    if valor_central > valor_vecino:
                        return false // No es un mínimo
        return true // Es un mínimo

3. Localización Precisa de Puntos Clave

Los puntos extremos detectados inicialmente son discretos. Para mejorar la estabilidad y precisión, se utiliza un ajuste tridimensional de una función cuadrática (expansión de Taylor) para localizar con exactitud la posición y escala subpíxel del punto clave. Este paso también filtra puntos clave con bajo contraste y elimina respuestas inestables en bordes.

3.1. Refinamiento de la Posición del Punto Clave

La expansión de Taylor de segundo orden de la función DoG D(x) alrededor del punto clave (x, y, σ) se utiliza para encontrar el extremo continuo. La derivada de esta expansión con respecto a x, y, σ, igualada a cero, permite calcular el desplazamiento (Δx, Δy, Δσ) desde el punto discreto. Si cualquier componente del desplazamiento es mayor que 0.5, el punto de interpolación se ha movido a un píxel vecino, y el proceso se repite en la nueva ubicación. Lowe sugiere un máximo de 5 iteraciones. Los puntos con un valor de función DoG muy bajo (menor a un umbral empírico, como 0.03 o 0.04/S) se descartan, ya que son susceptibles al ruido.

3.2. Eliminación de Respuestas de Borde

Los operadores DoG son propensos a generar fuertes respeustas a lo largo de los bordes, que son inestables y deben ser eliminadas. Para ello, se evalúa la matriz Hessiana 2x2 (H) en el punto clave. Las curvaturas principales de D (que son proporcionales a los valores propios de H) se utilizan para determinar si el punto clave está en un borde. Si los valores propios (α y β, donde α > β) están en una relación muy grande (α = rβ), el punto se considera una respuesta de borde. La condición para eliminar un punto de borde es:

Tr(H)^2 / Det(H) < (r+1)^2 / r

Donde Tr(H) es la traza de H y Det(H) es el determinante de H. Un valor común para r es 10.

Análisis de Código para la Interpolación de Taylor

El siguiente pseudocódigo muestra el bucle de refinamiento de puntos clave y la función de interpolación:

// Pseudocódigo para el bucle de refinamiento de puntos clave
function refinarPuntoClave(piramideDoG, octava_idx, nivel_idx, fila_idx, col_idx):
    max_iteraciones = 5 // SIFT_MAX_INTERP_STEPS
    umbral_borde = 5    // SIFT_IMG_BORDER

    for iter = 0 to max_iteraciones - 1:
        // Calcular los desplazamientos (delta_nivel, delta_fila, delta_columna) usando interpolación de Taylor
        (delta_nivel, delta_fila, delta_columna) = interp_paso(piramideDoG, octava_idx, nivel_idx, fila_idx, col_idx)

        // Si los desplazamientos son pequeños, se ha convergido
        if abs(delta_nivel) < 0.5 and abs(delta_fila) < 0.5 and abs(delta_columna) < 0.5:
            break

        // Actualizar la posición y el nivel del punto clave
        col_idx += redondear(delta_columna)
        fila_idx += redondear(delta_fila)
        nivel_idx += redondear(delta_nivel)

        // Verificar si el punto clave se ha movido fuera de los límites aceptables
        if nivel_idx < 1 or nivel_idx > numNivelesPorOctava or \
           col_idx < umbral_borde or fila_idx < umbral_borde or \
           col_idx > (anchoImagen(piramideDoG[octava_idx][0]) - umbral_borde) or \
           fila_idx > (altoImagen(piramideDoG[octava_idx][0]) - umbral_borde):
            return NULL // Punto no válido, descartar

    // Si no converge en max_iteraciones, también puede ser descartado
    // Realizar comprobación final de contraste y eliminación de bordes
    // ...

    return {fila_idx, col_idx, nivel_idx} // Devolver el punto clave refinado

// Función para un paso de interpolación
function interp_paso(piramideDoG, octava_idx, nivel_idx, fila_idx, col_idx):
    derivada_primera = calcularDerivada3D(piramideDoG, octava_idx, nivel_idx, fila_idx, col_idx)
    matriz_hessiana = calcularHessiana3D(piramideDoG, octava_idx, nivel_idx, fila_idx, col_idx)
    inversa_hessiana = calcularInversaMatriz(matriz_hessiana)

    // Calcular el desplazamiento X = -H_inv * dD
    desplazamiento = multiplicarMatriz(inversa_hessiana, derivada_primera, -1.0)

    return (desplazamiento[2], desplazamiento[1], desplazamiento[0]) // (delta_nivel, delta_fila, delta_columna)

4. Asignación de Orientación

Para lograr invarianza a la rotación, a cada punto clave se le asigna una orientación dominante basada en las propiedades locales de la imagen.

4.1. Cálculo del Gradiente del Punto Clave

En el punto clave detectado en la pirámide DoG, se calculan el módulo y la dirección del gradiente de los píxeles dentro de una ventana de 3σ_oct alrededor del punto, en la imagen gaussiana correspondiente a su escala. Los módulos del gradinete se ponderan con una función gaussiana centrada en el punto clave.

m(x,y) = sqrt( (L(x+1,y) - L(x-1,y))^2 + (L(x,y+1) - L(x,y-1))^2 )

theta(x,y) = atan2( (L(x,y+1) - L(x,y-1)), (L(x+1,y) - L(x-1,y)) )

Donde L es el valor de la imagen en la escala del punto clave.

4.2. Histograma de Orientación del Gradiente

Los gradientes de los píxeles dentro de la ventana se utilizan para construir un histograma de orientación. El rango de 0 a 360 grados se divide en 36 bins (cada uno de 10 grados). La magnitud ponderada de cada gradiente se añade al bin correspondiente a su orientación. El pico más alto del histograma representa la dirección principal del punto clave.

4.3. Determinación de la Orientación Principal y Auxiliar

La dirección del bin con el mayor valor en el histograma se convierte en la orientación principal del punto clave. Para mejorar la robustez, cualquier otro pico en el histograma que tenga un valor superior al 80% del pico principal se considera una orientación auxiliar. Esto puede resultar en la creación de múltiples puntos clave en la misma posición y escala, pero con diferentes orientaciones. Para una mayor precisión, se utiliza la interpolación parabólica para refinar el ángulo de orientación exacto para cada pico del histograma.

Análisis de Código para la Interpolación Parabólica

El pseudocódigo para asignar orientaciones, incluyendo la interpolación parabólica, es el siguiente:

// Pseudocódigo para interpolar picos del histograma
#define INTERP_PICO_HIST(l, c, r) (0.5 * ((l)-(r)) / ((l) - 2.0*(c) + (r)))

// Función para añadir características con orientaciones válidas
function anadirCaracteristicasDeOrientacion(lista_caracteristicas, histograma_orientaciones, num_bins, umbral_magnitud, caracteristica_base):
    PI2 = 2.0 * PI // Constante para 2*pi

    for i = 0 to num_bins - 1: // Iterar sobre cada bin del histograma (ej. 36 bins)
        izq = (i == 0) ? num_bins - 1 : i - 1 // Bin izquierdo (circular)
        der = (i + 1) % num_bins            // Bin derecho (circular)

        // Considerar un bin como un pico si es mayor que sus vecinos
        // Y su valor es mayor o igual al 80% del pico principal (umbral_magnitud)
        if histograma_orientaciones[i] > histograma_orientaciones[izq] and \
           histograma_orientaciones[i] > histograma_orientaciones[der] and \
           histograma_orientaciones[i] >= umbral_magnitud:

            // Interpolar para encontrar la orientación sub-bin
            bin_interpolado = i + INTERP_PICO_HIST(histograma_orientaciones[izq], histograma_orientaciones[i], histograma_orientaciones[der])

            // Asegurar que el ángulo esté dentro del rango [0, num_bins)
            if bin_interpolado < 0:
                bin_interpolado += num_bins
            elif bin_interpolado >= num_bins:
                bin_interpolado -= num_bins
            
            nueva_caracteristica = clonarCaracteristica(caracteristica_base)
            nueva_caracteristica.orientacion = (PI2 * bin_interpolado / num_bins) - PI // Convertir a radianes en [-PI, PI)
            
            anadirALista(lista_caracteristicas, nueva_caracteristica)
            liberarCaracteristica(nueva_caracteristica) // Asumiendo que anadirALista hace una copia

En este punto, cada punto clave tiene una posición (fila, columna), una escala (octava, nivel) y una orientación. Esto define una región SIFT única.

5. Generación del Descriptor de Característica

El paso final es crear un descriptor numérico para cada punto clave que sea robusto a las variaciones y distintivo para el emparejamiento.

5.1. Definición del Área para el Descriptor

Se selecciona una región alrededor del punto clave, orientada según su dirección principal. Esta región se divide en una cuadrícula de 4x4 subregiones (Lowe sugirió d=4). Para cada subregión, se calculan los histogramas de orientación. El tamaño de la ventana de imagen para el cálculo se define con un radio que considera la escala del punto clave y el tamaño de la cuadrícula.

5.2. Rotación de Ejes a la Orientación Principal

Para asegurar la invarianza a la rotación, el sistema de coordenadas local se rota de manera que el eje X coincida con la orientación principal del punto clave. Esto simplifica el cálculo de los gradientes en el nuevo sistema.

5.3. Generación de Histogramas de Gradiente

Dentro de cada subregión de la cuadrícula 4x4, se calculan los gradientes de los píxeles. Cada gradiente contribuye a un histograma de 8 bins de orientación para esa subregión. La contribución de cada píxel se pondera con una función gaussiana centrada en el punto clave y también espacialmente dentro de la subregión.

5.4. Interpolación Trilineal

Para que los descriptores sean más robustos a pequeños desplazamientos, la contribución de cada muestra de gradiente se distribuye a los bins vecinos mediante interpolación trilineal. Una muestra de gradiente contribuye a los 8 bins más cercanos en un cubo 2x2x2 en el espacio (subregión x, subregión y, bin de orientación).

5.5. Vector de Característica

Los 16 histogramas (4x4 subregiones, cada una con 8 bins) se concatenan para formar un vector de 128 dimensiones (4x4x8=128). Este es el descriptor SIFT. Para la invarianza a cambios de iluminación, el vector se normaliza a una longitud unitaria.

5.6. Umbralización del Descriptor

Para mitigar el efecto de cambios no lineales de iluminación o saturación del sensor, los valores del descriptor se umbralizan. Cualquier componente del vector con un valor superior a un umbral (comúnmente 0.2 después de la normalización) se recorta a ese umbral. Después de la umbralización, el vector se normaliza nuevamente para mejorar la distintividad.

Análisis de Código para la Interpolación Trilineal del Descriptor

El siguiente pseudocódigo demuestra la interpolación trilineal utilizada para generar los descriptores:

// Pseudocódigo para la interpolación trilineal de una entrada del histograma del descriptor
function interpolarEntradaHist(descriptor_hist, fila_bin, col_bin, orient_bin, magnitud_gradiente, tamano_grid, num_bins_orient):
    fila0 = piso(fila_bin)    // Parte entera del bin de fila
    col0 = piso(col_bin)    // Parte entera del bin de columna
    orient0 = piso(orient_bin) // Parte entera del bin de orientación

    dr = fila_bin - fila0    // Parte fraccionaria de fila
    dc = col_bin - col0    // Parte fraccionaria de columna
    do = orient_bin - orient0 // Parte fraccionaria de orientación

    for r_offset = 0 to 1: // Recorrido en la dimensión de fila (0 o 1)
        rb = fila0 + r_offset
        if rb >= 0 and rb < tamano_grid: // Si está dentro del grid del descriptor
            ponderacion_r = magnitud_gradiente * ((r_offset == 0) ? (1.0 - dr) : dr)
            
            for c_offset = 0 to 1: // Recorrido en la dimensión de columna (0 o 1)
                cb = col0 + c_offset
                if cb >= 0 and cb < tamano_grid: // Si está dentro del grid del descriptor
                    ponderacion_c = ponderacion_r * ((c_offset == 0) ? (1.0 - dc) : dc)
                    
                    for o_offset = 0 to 1: // Recorrido en la dimensión de orientación (0 o 1)
                        ob = (orient0 + o_offset) % num_bins_orient // Bin de orientación (circular)
                        ponderacion_o = ponderacion_c * ((o_offset == 0) ? (1.0 - do) : do)
                        
                        // Añadir la ponderación al bin correspondiente del histograma
                        descriptor_hist[rb][cb][ob] += ponderacion_o

Con estos pasos, el algoritmo SIFT extrae características robustas de la imagen, que luego pueden utilizarse para tareas de emparejamiento y reconocimiento de imágenes.