Este artículo explora la implementación del algoritmo SIFT (Scale-Invariant Feature Transform) en Python, detallando sus componentes clave y proporcionando ejemplos de código. Se aborda la construcción del espacio de escala, la detección de puntos clave, la asignación de orientaciones y la generación de descriptores de características.

El algoritmo SIFT es fundamental para la extracción de características de imágenes que son invariantes a la escala, rotación y cambios de iluminación. Su objetivo principal es identificar puntos de interés únicos en una imagen y generar un descriptor para cada uno, facilitando la comparación y el emparejamiento entre diferentes imágenes.

La implementación del SIFT abarca varios pasos:

1. Construcción del Espacio de Escala

Este paso implica la creación de una pirámide gaussiana y una pirámide de diferencias de gaussianas (DoG). La pirámide gaussiana se genera aplicando filtros gaussianos con diferentes desviaciones estándar (σ) a la imagen original, y luego submuestreando la imagen para crear octavas. La pirámide DoG se obtiene calculando la diferencia entre imágenes gaussianas consecutivas dentro de cada octava. Esto ayuda a detectar características a diferentes escalas.

Ejemplo de Estructura de Clase para Parámetros SIFT:

class CSiftConfig:
    def __init__(self, num_octave_layers, initial_sigma, contrast_threshold, edge_threshold, scale_factor):
        self.num_scale = num_octave_layers  # Número de capas por octava
        self.sigma = initial_sigma          # Sigma inicial para el filtro Gaussiano
        self.contrast_t = contrast_threshold # Umbral de contraste para eliminar puntos débiles
        self.eigenvalue_r = edge_threshold   # Ratio de valores propios para eliminar bordes
        self.scale_factor = scale_factor     # Factor de escala para la orientación
        self.num_bins = 36                   # Número de bins para el histograma de orientación
        self.peak_ratio = 0.8                # Ratio para considerar direcciones secundarias

Función para Preprocesar la Imagen Inicial:

def preprocess_image(image, target_sigma, camera_sigma=0.5):
    # Calcula el sigma intermedio necesario considerando la posibleGaussianBlur inicial de la cámara
    sigma_mid = np.sqrt(target_sigma**2 - (2 * camera_sigma)**2)
    # Duplica el tamaño de la imagen y aplica el filtro Gaussiano
    resized_image = cv2.resize(image, (image.shape[1] * 2, image.shape[0] * 2), interpolation=cv2.INTER_LINEAR)
    blurred_image = cv2.GaussianBlur(resized_image, (0, 0), sigmaX=sigma_mid, sigmaY=sigma_mid)
    return blurred_image

2. Detección de Puntos Clave (Extremas)

Los puntos clave se detectan identificando los extremos (máximos y mínimos locales) en la pirámide DoG. Para cada píxel en una capa DoG, se compara con sus vecinos en la misma capa y en las capas adyacentes (superior e inferior). Esto implica comparar cada punto con 26 vecinos en total. Se aplica un umbral de contraste para descartar puntos con respuestas débiles.

Para refinar la ubicación de los puntos clave, se utiliza la interpolación usando expansiones de Taylor de segundo orden de la función DoG. Si la interpolación resulta en un desplazamiento significativo o si el valor DoG interpolado es demasiado bajo, el punto se descarta. Además, se eliminan los puntos que corresponden a bordes, utilizando la matriz Hessiana. Se calcula la relación entre el determinante y la traza de la matriz Hessiana proyectada sobre la dirección de máxima curvatura.

Función para Construir la Pirámide Gaussiana:

def build_gaussian_pyramid(base_image, config: CSiftConfig):
    pyramid = []
    current_image = base_image.copy()
    num_octaves = int(np.round(np.log(min(current_image.shape)) / np.log(2))) - 1 # Calcula el número de octavas

    for i in range(num_octaves):
        octave_layers = build_octave_layers(current_image, config.num_scale, config.sigma)
        pyramid.append(octave_layers)
        # La siguiente octava se obtiene submuestreando la tercera capa desde el final
        current_image = octave_layers[-3]
        current_image = cv2.resize(current_image, (int(current_image.shape[1] / 2), int(current_image.shape[0] / 2)), interpolation=cv2.INTER_NEAREST)
    return pyramid

def build_octave_layers(img_src, num_layers, base_sigma):
    layers = []
    layers.append(img_src) # La imagen de entrada ya tiene el sigma inicial aplicado
    scale_k = 2**(1.0/num_layers)

    for i in range(1, num_layers + 3): # Necesitamos num_layers + 2 capas para obtener num_layers DoG layers
        prev_sigma_total = base_sigma * (scale_k**(i-1))
        current_sigma_total = base_sigma * (scale_k**i)
        # Sigma necesario para el filtro Gaussiano actual
        sigma_needed = np.sqrt(current_sigma_total**2 - prev_sigma_total**2)
        blurred_layer = cv2.GaussianBlur(layers[-1], (0,0), sigmaX=sigma_needed, sigmaY=sigma_needed)
        layers.append(blurred_layer)
    return layers

Función para Construir la Pirámide DoG y Detectar Puntos Clave:

def detect_keypoints(gaussian_pyramid, config: CSiftConfig):
    keypoints = []
    contrast_threshold_scaled = np.floor(0.5 * config.contrast_t / config.num_scale * 255)

    for octave_idx, octave in enumerate(gaussian_pyramid):
        dog_octave = []
        for layer_idx in range(len(octave) - 1):
            diff = octave[layer_idx+1] - octave[layer_idx]
            dog_octave.append(diff)

        for s in range(1, len(dog_octave) - 1):
            bottom_dog, mid_dog, top_dog = dog_octave[s-1], dog_octave[s], dog_octave[s+1]
            # Iterar sobre los píxeles de la capa intermedia, excluyendo bordes
            border_width = 5
            for r in range(border_width, bottom_dog.shape[0] - border_width):
                for c in range(border_width, bottom_dog.shape[1] - border_width):
                    pixel_value = mid_dog[r, c]
                    # Verificar si es un extremo local (mayor que vecinos o menor que vecinos)
                    is_local_extremum = is_extremum(bottom_dog, mid_dog, top_dog, r, c, pixel_value, contrast_threshold_scaled)

                    if is_local_extremum:
                        # Intentar ajustar la posición del extremo para obtener coordenadas más precisas
                        fitted_keypoint = refine_extremum_location(dog_octave, s, r, c, border_width, octave_idx, config, gaussian_pyramid[octave_idx][s])
                        if fitted_keypoint:
                            keypoints.append(fitted_keypoint)
    return keypoints

def is_extremum(bottom_layer, mid_layer, top_layer, r, c, value, threshold):
    # Compara el valor central con sus 8 vecinos en las 3 capas (total 26 vecinos)
    is_max = value > threshold and np.all(value > mid_layer[r-1:r+2, c-1:c+2]) and np.all(value > bottom_layer[r-1:r+2, c-1:c+2]) and np.all(value > top_layer[r-1:r+2, c-1:c+2])
    is_min = value < -threshold and np.all(value < mid_layer[r-1:r+2, c-1:c+2]) and np.all(value < bottom_layer[r-1:r+2, c-1:c+2]) and np.all(value < top_layer[r-1:r+2, c-1:c+2])
    return is_max or is_min

# Funciones para refine_extremum_location (fit_extremum, get_gradient, get_hessian) y remove_edge_response deben ser implementadas detalladamente.
# Estas funciones involucran cálculos matriciales y derivadas parciales complejas.

3. Asignación de Orientación

Para lograr la invarianza a la rotación, a cada punto clave se le asigna una o más orientaciones. Esto se hace calculando el gradiente de magnitud y dirección en una vecindad alrededor del punto clave, en la imagen gaussiana correspondiente a su escala. Se construye un histograma de orientaciones (generalmente con 36 bins) ponderado por la magnitud del gradiente y un peso gaussiano. La orientación principal se toma del pico más alto del histograma suavizado. Si hay otros picos que superan el 80% del pico principal, se crean puntos clave adicionales para esas orientaciones.

Función Esquemática para Calcular Orientación:

def compute_orientation(keypoint, octave_image, config: CSiftConfig):
    orientations = []
    scale = keypoint.size # El tamaño del KeyPoint ya contiene información de escala
    radius = int(round(config.radius_factor * config.scale_factor * scale))
    weight_sigma_sq = -0.5 / ((config.scale_factor * scale) ** 2)
    histogram = np.zeros(config.num_bins)

    center_x = int(round(keypoint.pt[0] * octave_image.shape[1]))
    center_y = int(round(keypoint.pt[1] * octave_image.shape[0]))

    for y in range(-radius, radius + 1):
        for x in range(-radius, radius + 1):
            img_y, img_x = center_y + y, center_x + x
            if 0 < img_y < octave_image.shape[0] - 1 and 0 < img_x < octave_image.shape[1] - 1:
                dx = float(octave_image[img_y, img_x + 1]) - float(octave_image[img_y, img_x - 1])
                dy = float(octave_image[img_y - 1, img_x]) - float(octave_image[img_y + 1, img_x])
                magnitude = np.sqrt(dx**2 + dy**2)
                orientation_rad = np.arctan2(dy, dx)
                orientation_deg = np.rad2deg(orientation_rad)
                if orientation_deg < 0: orientation_deg += 360

                orientation_bin = int(round(orientation_deg / (360.0 / config.num_bins))) % config.num_bins
                weight = np.exp(weight_sigma_sq * (x**2 + y**2))
                histogram[orientation_bin] += magnitude * weight

    # Suavizado del histograma y búsqueda de picos
    # ... (implementación del suavizado y la interpolación parabólica para encontrar orientaciones principales y secundarias)
    # Crear nuevos KeyPoints para cada orientación principal/secundaria
    # ...
    return orientations # Lista de KeyPoints con orientación

4. Generación de Descriptores de Características

Para cada punto clave con su orientación asignada, se genera un descriptor de 128 dimensiones. Se define una ventana de vecindad alrededor del punto clave, cuya escala está determinada por el tamaño del punto clave. Esta ventana se divide en una cuadrícula de subregiones (comúnmente 4x4). Se rota el sistema de coordenadas para que el eje principal del punto clave se alinee con el eje X. Luego, se calculan los gradientes de magnitud y dirección para cada píxel dentro de la vecindad rotada. Estos gradientes se ponderan y se distribuyen en un histograma de 8 orientaciones por cada subregión de la cuadrícula (4x4x8 = 128 dimensiones). Finalmente, el vector descriptor resultante se normaliza para lograr la invarianza a la iluminación y se trunca para reducir la sensibilidad a cambios no lineales de intensidad.

Función Esquemática para Calcular Descriptores:

def compute_descriptors(keypoints, gaussian_pyramid, config: CSiftConfig):
    descriptors = []
    for kp in keypoints:
        # Obtener imagen gaussiana y escala
        octave = kp.octave & 255 # Información del octavo
        layer = (kp.octave >> 8) & 255 # Información de la capa
        image = gaussian_pyramid[octave][layer]
        scale = kp.size
        orientation = kp.angle
        
        # Definir parámetros del descriptor (ej. 4x4 subregiones, 8 bins de orientación)
        num_subregions = 4
        num_orientation_bins = 8
        descriptor_dims = num_subregions * num_subregions * num_orientation_bins
        descriptor = np.zeros(descriptor_dims)

        # Rotar el sistema de coordenadas y calcular gradientes
        # ... (implementación detallada del cálculo de gradientes y la rotación)

        # Distribuir gradientes en el histograma 3D (4x4x8) usando interpolación
        # ...

        # Aplanar el histograma y normalizar el vector descriptor
        flat_descriptor = descriptor.flatten()
        norm_val = np.linalg.norm(flat_descriptor)
        if norm_val > 1e-6: # Evitar división por cero
             flat_descriptor /= norm_val

        # Truncar valores y volver a escalar (ej. a 0-255)
        flat_descriptor[flat_descriptor > config.des_max_value] = config.des_max_value
        final_descriptor = np.round(512 * flat_descriptor)
        # Asegurar que los valores estén en el rango [0, 255]
        final_descriptor[final_descriptor < 0] = 0
        final_descriptor[final_descriptor > 255] = 255
        
        descriptors.append(final_descriptor)
    return descriptors

5. Emparejamiento de Descriptores

Una vez que se tienen los descriptores de dos imágenes, se pueden emparejar utilizando métricas de distancia como la distancia Euclidiana. Se utiliza comúnmente el algoritmo k-NN (k-Nearest Neighbors) para encontrar los mejores emparejamientos. El ratio de distancia entre el mejor y el segundo mejor vecino se usa a menudo como un criterio para aceptar un emparejamiento, descartando aquellos que son ambiguos.

Ejemplo de Emparejamiento con FLANN:

def match_descriptors(des1, des2, kpts1, kpts2, img1, img2, min_match_count=10):
    FLANN_INDEX_KDTREE = 1
    index_params = dict(algorithm=FLANN_INDEX_KDTREE, trees=5)
    search_params = dict(checks=50)
    flann = cv2.FlannBasedMatcher(index_params, search_params)
    
    matches = flann.knnMatch(np.array(des1).astype(np.float32), np.array(des2).astype(np.float32), k=2)
    
    good_matches = []
    for m, n in matches:
        if m.distance < 0.75 * n.distance: # Ratio test
            good_matches.append(m)

    # Dibujar las buenas coincidencias si hay suficientes
    if len(good_matches) > min_match_count:
        # ... (código para dibujar líneas entre puntos correspondientes)
        # Usar cv2.drawMatches o una implementación manual
        pass # Placeholder
    else:
        print(f"Not enough matches are found - {len(good_matches)}/{min_match_count}")
        return None # O una imagen vacía, o un indicador de fallo

    # Construir la imagen resultante con las líneas de emparejamiento
    # ...
    return result_image # Imagen con las correspondencias dibujadas

La implementación completa de SIFT es compleja y requiere un manejo cuidadoso de las matemáticas subyacentes, especialmente en el cálculo de derivadas, interpolación y transformaciones. Si bien es posible implementar SIFT desde cero, el uso de bibliotecas optimizadas como OpenCV es generalmente más práctico para aplicaciones en tiempo real.

El artículo original también incluye una sección sobre la comparación de la implementación manual con la versión de OpenCV, destacando la superioridad en velocidad y número de características detectadas de la implementación de OpenCV.

Adicionalmente, se exploran escenarios como la rotación, traslación e inclinación de la cámara, analizando cómo afectan los gradientes de los puntos clave. Se concluye que la inclinación de la cámara tiende a generar mayores variaciones en los gradientes.

Nota: Las implementaciones de código proporcionadas aquí son esquemáticas y requieren la adición de las funciones matemáticas detalladas (como refine\_extremum\_location, compute\_orientation y compute\_descriptors) para una funcionalidad completa.