1. Introducción

La reconstrucción de imágenes es el proceso de generar una representación visual a partir de datos medidos indirectamente (por ejemplo, proyecciones de rayos X, datos de resonancia magnética o señales de tomografía por emisión de positrones). Este proceso es crítico en áreas como la medicina, la inspección industrial y la visión por computadora.

En los últimos años, la combinación de métodos clásicos con técnicas de aprendizaje profundo ha permitido alcanzar una calidad de imagen sin precedentes y reducir drásticamente los tiempos de cálculo.

2. Fundamentos Matemáticos

En la mayoría de los sistemas de adquisición se modela la relación entre la imagen original f(x,y) y los datos medidos g(s,θ) mediante la Transformada de Radón:

g(s,θ) = ∫∫ f(x,y) δ(s - x cosθ - y sinθ) dx dy

La tarea de reconstrucción consiste en invertir esta integral. Las técnicas más usadas son:

• Reconstrucción directa (Filtrado y Retroproyección)
• Iterativa (SART, MLEM, OSEM)
• Compresión y Reconstrucción (Compressed Sensing)
• Redes Neuronales Profundas (U‑Net, AUTOMAP)

3. Comparativa de Algoritmos (Dos Columnas)

4. Implementación en Python

Recomendamos usar entornos conda o venv para aislar dependencias.

# Crear entorno con Conda
conda create -n img-recon python=3.11 numpy scipy scikit-image matplotlib pytorch torchvision tqdm
conda activate img-recon

4.1. Reconstrucción por Retroproyección Filtrada (FBP)

import numpy as np
import matplotlib.pyplot as plt
from skimage.transform import radon, iradon
# Generar una imagen de prueba
image = np.zeros((256, 256))
rr, cc = np.ogrid[:256, :256]
mask = (rr-128)**2 + (cc-128)**2 < 40**2
image[mask] = 1
# Simular proyecciones
theta = np.linspace(0., 180., max(image.shape), endpoint=False)
sinogram = radon(image, theta=theta, circle=True)
# Reconstrucción FBP
reconstruction_fbp = iradon(sinogram, theta=theta, filter_name='ramp', circle=True)
# Visualizar
fig, axes = plt.subplots(1, 3, figsize=(12, 4))
axes[0].set_title('Original')
axes[0].imshow(image, cmap='gray')
axes[1].set_title('Sinograma')
axes[1].imshow(sinogram, cmap='gray', aspect='auto')
axes[2].set_title('Reconstrucción FBP')
axes[2].imshow(reconstruction_fbp, cmap='gray')
plt.tight_layout()
plt.show()

Este ejemplo utiliza scikit-image para la Transformada de Radón y su inversa. Es ideal para prototipos y datos bien muestreados.

4.2. Reconstrucción Iterativa (SART)

from skimage.transform import iradon_sart
# SART necesita un número de iteraciones y un ángulo inicial
reconstruction_sart = iradon_sart(sinogram, theta=theta)
for _ in range(4):
    reconstruction_sart = iradon_sart(sinogram, theta=theta, image=reconstruction_sart)
plt.figure(figsize=(5,5))
plt.title('Reconstrucción SART (5 iteraciones)')
plt.imshow(reconstruction_sart, cmap='gray')
plt.axis('off')
plt.show()

SART (Simultaneous Algebraic Reconstruction Technique) converge rápidamente con pocos ángulos, pero es más costoso que FBP.

4.3. Reconstrucción basada en Deep Learning (U‑Net)

Para casos donde la calidad de la adquisición es limitada, entrenar una red U‑Net que aprenda la transformación de sinograma → imagen puede superar a los métodos clásicos.

import torch
import torch.nn as nn
from torch.utils.data import DataLoader, Dataset
class UNet(nn.Module):
    def __init__(self, in_channels=1, out_channels=1, features=[64,128,256,512]):
        super(UNet, self).__init__()
        # Encoder
        self.encoders = nn.ModuleList()
        for f in features:
            self.encoders.append(nn.Sequential(
                nn.Conv2d(in_channels, f, kernel_size=3, padding=1),
                nn.ReLU(inplace=True),
                nn.Conv2d(f, f, kernel_size=3, padding=1),
                nn.ReLU(inplace=True)
            ))
            in_channels = f
        # Decoder
        self.decoders = nn.ModuleList()
        for f in reversed(features):
            self.decoders.append(nn.Sequential(
                nn.ConvTranspose2d(f*2, f, kernel_size=2, stride=2),
                nn.ReLU(inplace=True),
                nn.Conv2d(f, f, kernel_size=3, padding=1),
                nn.ReLU(inplace=True),
                nn.Conv2d(f, f, kernel_size=3, padding=1),
                nn.ReLU(inplace=True)
            ))
        self.bottleneck = nn.Sequential(
            nn.Conv2d(features[-1], features[-1]*2, kernel_size=3, padding=1),
            nn.ReLU(inplace=True),
            nn.Conv2d(features[-1]*2, features[-1]*2, kernel_size=3, padding=1),
            nn.ReLU(inplace=True)
        )
        self.final_conv = nn.Conv2d(features[0], out_channels, kernel_size=1)
    def forward(self, x):
        enc_features = []
        for enc in self.encoders:
            x = enc(x)
            enc_features.append(x)
            x = nn.MaxPool2d(2)(x)
        x = self.bottleneck(x)
        for dec, enc_feat in zip(self.decoders, reversed(enc_features)):
            x = dec(torch.cat([x, enc_feat], dim=1))
        return self.final_conv(x)
# Dummy dataset (replace with real sinograms & ground‑truth images)
class SinogramDataset(Dataset):
    def __init__(self, sinograms, images):
        self.sinograms = sinograms
        self.images = images
    def __len__(self):
        return len(self.sinograms)
    def __getitem__(self, idx):
        return torch.from_numpy(self.sinograms[idx]).unsqueeze(0).float(), \
               torch.from_numpy(self.images[idx]).unsqueeze(0).float()
# Entrenamiento rápido (ejemplo)
model = UNet().cuda()
criterion = nn.MSELoss()
optimizer = torch.optim.Adam(model.parameters(), lr=1e-4)
# DataLoader placeholder
train_loader = DataLoader(SinogramDataset(train_sinograms, train_images), batch_size=8, shuffle=True)
for epoch in range(10):
    model.train()
    epoch_loss = 0.0
    for sinogram_batch, image_batch in train_loader:
        sinogram_batch = sinogram_batch.cuda()
        image_batch = image_batch.cuda()
        optimizer.zero_grad()
        output = model(sinogram_batch)
        loss = criterion(output, image_batch)
        loss.backward()
        optimizer.step()
        epoch_loss += loss.item()
    print(f'Epoch {epoch+1}, Loss: {epoch_loss/len(train_loader):.4f}')

El entrenamiento requiere cientos de pares sinograma‑imagen. Se recomienda usar GPU con CUDA y aplicar data augmentation (rotaciones, ruido Poisson) para robustez.

5. Buenas Prácticas, Optimización y Seguridad

5.1. Rendimiento

• Vectorización: Utiliza numpy y torch en modo batch para evitar bucles en Python.
• GPU: Para algoritmos iterativos y redes neuronales, habilita torch.cuda.amp (mixed precision) para reducir la memoria y acelerar el cómputo.
• Compresión de datos: Emplea formatos .npz o HDF5 para almacenar sinogramas sin pérdidas.

5.2. Seguridad y Privacidad

• En entornos médicos, cifra los sinogramas en reposo usando AES‑256.
• Aplica de‑identificación antes de compartir datos con modelos de IA (remover metadatos PHI).
• Control de accesos mediante OAuth2 y tokens JWT en APIs que sirvan los datos de reconstrucción.

5.3. Depuración (Troubleshooting)

• Artefactos de cuarteo: Verifica la uniformidad del ángulo de adquisición; usa np.unwrap para corregir desalineaciones.
• Sobre‑suavizado en TV: Ajusta el parámetro λ; valores demasiado altos eliminan detalles finos.
• Instabilidad de entrenamiento DL: Normaliza los sinogramas (media 0, var 1) y emplea torch.nn.BatchNorm2d.

6. Escalabilidad y Compatibilidad

Los algoritmos pueden escalar desde desktops hasta clusters HPC usando MPI4Py o Dask para distribuir batches de sinogramas.

# Ejemplo de paralelismo con Dask
from dask import array as da
sinograms = da.from_array(all_sinograms, chunks=(8, 512, 512))
recon = sinograms.map_blocks(lambda s: iradon(s, theta=theta, filter_name='ramp'))
result = recon.compute()

Esta estrategia permite procesar terabytes de datos en minutos cuando se combina con nodos GPU.