Búsqueda del Vecino más Próximo

Recorre el array desde el final hacia el principio:

for(int i=f.dim()-2; i>=0; i=i-1){
    /// Busca el intervalo donde x[i] es menor que x0 (asumiendo x ordenado)
    if( x[i] < x0 ){
        /// Compara distancias
        if( (x0-x[i]) <= mn_abs(x0-x[i+1]) ) return f[i]; /// Retorna el valor a la izquierda
        else return f[i+1]; /// Retorna el valor a la derecha
    }
}
return f[0];

Interpolación Lineal por Tramos

Bucle que recorre los intervalos desde el final hacia el principio:

for(int i=f.dim()-2; i>=0; i=i-1){
    /// Busca el primer x[i] que sea menor que x0 (encuentra el intervalo [x_i, x_i+1])
    if( x[i] < x0 ){
        /// Aplica la fórmula de la recta
        return f[i]+(x0-x[i])*(f[i+1]-f[i])/(x[i+1]-x[i]);
    }
}
/// Si x0 es menor que todos los puntos, extrapola usando los dos primeros puntos (x0 y x1)
return f[0]+(x0-x[0])*(f[1]-f[0])/(x[1]-x[0]);

Cálculo de Spline Cuadrático

Validación y Configuración Inicial

Los vectores deben tener igual tamaño y al menos 2 dimensiones:

if(f.dim()!=x.dim() || f.dim() < 2 ) return -1;

for(int k=1; k < x.dim(); k++){
    // Los puntos x deben estar ordenados ascendentemente
    if(x[k] <= x[k-1]) return -1;
}

a = Array1D<real>(x.dim()-1);
b = Array1D<real>(x.dim()-1);
c = Array1D<real>(x.dim()-1);

Cálculo del Primer Tramo

Cálculo del primer tramo:

a[0]=f[0]; // Debe pasar por el punto inicial

// Se despeja b[0]
b[0]=(f[1]-f[0]-c0*(x[1]-x[0])*(x[1]-x[0]))/(x[1]-x[0]);
c[0]=c0;

Bucle para el Resto de Tramos

for(int i=1; i < a.dim(); i++){
    a[i]=f[i];

    // La pendiente al final del tramo anterior debe ser igual a la pendiente al inicio del actual
    b[i]=b[i-1]+2*c[i-1]*(x[i]-x[i-1]);

    // Se despeja c[i] para que el spline llegue a f[i+1] en el punto x[i+1]
    c[i]=(f[i+1]-a[i]-b[i]*(x[i+1]-x[i]))/((x[i+1]-x[i])*(x[i+1]-x[i]));
}
return 0;

Evaluación del Spline Cuadrático

Recorre los tramos desde el último (derecha) hacia el primero (izquierda):

for(int i=a.dim()-1; i>=0; i--){
    // Evaluación del polinomio cuadrático
    if(x[i] < x0) return a[i]+(b[i]+c[i]*(x0-x[i]))*(x0-x[i]);
}
return a[0]+(b[0]+c[0]*(x0-x[0]))*(x0-x[0]); // Caso de extrapolación

Spline Cuadrático con Condición de Borde Final

Verifica que el tamaño sea correcto y que x esté ordenado:

if(f.dim()!=x.dim() || f.dim() < 2 ) return -1;

for(int k=1; k < x.dim(); k++){
    if(x[k] <= x[k-1]) return -1;
}

a = Array1D<real>(x.dim()-1);
b = Array1D<real>(x.dim()-1);
c = Array1D<real>(x.dim()-1);

/// Asigna el valor de la función al coeficiente a del último tramo
a[a.dim()-1]=f[a.dim()-1];

/// Calcula el coeficiente b del último tramo usando la condición de borde cLast
b[a.dim()-1]=(f[a.dim()]-f[a.dim()-1]-cLast*(x[a.dim()]-x[a.dim()-1])*(x[a.dim()]-x[a.dim()-1]))/(x[a.dim()]-x[a.dim()-1]);

/// Asigna el valor de cLast al coeficiente c del último tramo
c[a.dim()-1]=cLast;

/// Bucle para calcular los coeficientes de los tramos restantes hacia atrás
for(int i=a.dim()-2; i>=0; i--){
    a[i]=f[i];
    b[i]=2*(f[i+1]-f[i])/(x[i+1]-x[i])-b[i+1];
    c[i]=(b[i+1]*(x[i+1]-x[i])+f[i]-f[i+1])/((x[i+1]-x[i])*(x[i+1]-x[i]));
}
return 0;

Evaluación del Spline (Hacia Atrás)

Bucle que recorre los tramos del spline desde el último hacia el primero:

for(int i=a.dim()-1; i>=0; i--){
    if(x[i] < x0) return a[i]+b[i]*(x0-x[i])+c[i]*(x0-x[i])*(x0-x[i]);
}
return a[0]+b[0]*(x0-x[0])+c[0]*(x0-x[0])*(x0-x[0]);

Implementación de Spline Cúbico

Define el vector h y la diagonal principal M, la parte derecha B y las diagonales L y U del sistema tridiagonal:

Array1D<real> h(x.dim()-1);
Array1D<real> M(x.dim(),1.);
Array1D<real> B(x.dim());
Array1D<real> U(x.dim()-1,0.);
Array1D<real> L(x.dim()-1,0.);

/// Calcula el ancho de cada intervalo entre nodos consecutivos
for(int k=0; k < h.dim(); k++) h[k]=x[k+1]-x[k];

/// Establece las condiciones de contorno para el primer y último punto del sistema
B[0]=c0;
B[B.dim()-1]=cN;

/// Construye las ecuaciones internas del sistema tridiagonal para los coeficientes c
for(int k=1; k < M.dim()-1; k++){
    L[k-1]=h[k-1];
    M[k]=2*(h[k-1]+h[k]);
    U[k]=h[k];
    B[k]=3.*(f[k+1]-f[k])/h[k] - 3.*(f[k]-f[k-1])/h[k-1];
}
c = solucion_sistema(L,M,U,B);

a = Array1D<real>(x.dim()-1);
b = Array1D<real>(x.dim()-1);
d = Array1D<real>(x.dim()-1);

/// Calcula los coeficientes restantes
for(int k=0; k < a.dim(); k++){
    a[k]=f[k];
    d[k]=(c[k+1]-c[k])/(3*h[k]);
    b[k]=(f[k+1]-f[k])/h[k]-h[k]*(2*c[k]+c[k+1])/3.;
}
return 0;

Evaluación del Spline Cúbico

for(int i=a.dim()-1; i>=0; i--){
    if(x[i] < x0) return ((d[i]*(x0-x[i])+c[i])*(x0-x[i])+b[i])*(x0-x[i])+a[i];
}
return ((d[0]*(x0-x[0])+c[0])*(x0-x[0])+b[0])*(x0-x[0])+a[0];

Regresión Lineal por Mínimos Cuadrados

Obtiene el número de puntos de la nube a partir de la dimensión del vector x:

int N=x.dim();

/// Verifica que los vectores x e y tengan el mismo tamaño y al menos 2 puntos para operar
if(y.dim()!=N || N < 2) return(-1);

real suma_x=0, suma_y=0, suma_xy=0, suma_x2=0;

/// Recorre todos los puntos para calcular las sumas de x, y, el producto xy y el cuadrado x^2
for(int k=0; k < N; k++){
    suma_x+=x[k];
    suma_y+=y[k];
    suma_xy+=x[k]*y[k];
    suma_x2+=x[k]*x[k];
}

/// Calcula el denominador común para las fórmulas de los coeficientes a y b
real denominador=N*suma_x2-suma_x*suma_x;

/// Si el denominador es 0, los puntos están alineados verticalmente y no hay solución única
if(denominador==0.) return(-1);

/// Calcula la pendiente 'a' y el término independiente o punto de corte 'b' con el eje y
a=(N*suma_xy-suma_x*suma_y)/denominador;
b=(suma_x2*suma_y-suma_xy*suma_x)/denominador;
return(0);

Interpolación Polinómica de Newton

Cálculo de Coeficientes (Diferencias Divididas)

Verifica que los vectores X y F tengan el mismo número de elementos:

if(X.dim()!=F.dim()) return(Array1D<real>());
Array1D<real> A(X.dim());

/// Define N como el grado máximo del polinomio
int N=X.dim()-1;
Array1D<real> B(A.dim());

/// Inicializa el vector B con los valores de la función
for(int k=0; k <= N; k++){
    B[k]=F[k];
}

/// El primer coeficiente del polinomio es el primer valor de la función
A[0]=F[0];

/// Bucle principal para calcular las diferencias divididas de orden 1 hasta N
for(int k=1; k <= N; k++){
    for(int l=0; l <= (N-k); l++){
        if(X[k+l]==X[l]){
            return(Array1D<real>());
        }

        /// Aplica la fórmula de diferencias divididas:
        B[l]=(B[l+1]-B[l])/(X[k+l]-X[l]);
    }

    /// Almacena la primera diferencia dividida de cada orden en el vector de coeficientes A
    A[k]=B[0];
}
return(A);

Evaluación del Polinomio de Newton

int N=A.dim()-1;
real E=A[N];
for(int k=N-1; k>=0; k--){
    E=E*(x0-X[k])+A[k];
}
return E;

Escalado de Imágenes y Remuestreo

Algoritmo del Vecino más Próximo

Calcula las nuevas dimensiones de la imagen multiplicando las originales por el factor z:

int dim1New=A.dim1()*z;
int dim2New=A.dim2()*z;

/// Si las dimensiones resultantes son nulas o negativas, devuelve una imagen vacía
if(dim1New <= 0 || dim2New <= 0) return Array2D<real>();

Array2D<real> B(dim1New,dim2New);

// Inicia el recorrido por cada píxel de la nueva imagen
for(int ip=0; ip < B.dim1(); ip++){
    for(int jp=0; jp < B.dim2(); jp++){
        real x=ip/z;
        real y=jp/z;
        int i=x;
        int j=y;

        /// Implementa el redondeo al vecino más próximo para la coordenada x (fila)
        if((x-i)>0.5 && i+1 < A.dim1()) i=i+1;

        /// Implementa el redondeo al vecino más próximo para la coordenada y (columna)
        if((y-j)>0.5 && j+1 < A.dim2()) j=j+1;

        /// Copia el valor del píxel más cercano de la imagen original a la posición actual en B
        B[ip][jp]=A[i][j];
    }
}
/// Devuelve la imagen escalada tras completar todos los píxeles
return B;

Interpolación Bilineal

Calcula la nueva dimensión de filas y columnas multiplicando la original por el factor de zoom:

int dim1New=A.dim1()*z;
int dim2New=A.dim2()*z;

/// Si alguna dimensión resultante es no positiva, devuelve una matriz vacía
if(dim1New <= 0 || dim2New <= 0) return Array2D<real>();

Array2D<real> B(dim1New,dim2New);

/// Inicia el bucle para recorrer cada fila de la imagen de salida
for(int ip=0; ip < B.dim1(); ip++){
    for(int jp=0; jp < B.dim2(); jp++){

        /// Mapeo a coordenadas originales y cálculo de índices base
        real x=ip/z;
        real y=jp/z;
        int i=x;
        int j=y;

        /// Selección de los 4 vecinos (limitando bordes)
        int i1=min(i+1,A.dim1()-1);
        int j1=min(j+1,A.dim2()-1);

        /// Cálculo de pesos (distancia fraccionaria)
        real dx=x-i;
        real dy=y-j;

        /// Fórmula bilineal: promedio ponderado de los 4 vecinos
        B[ip][jp] = (1.-dy) * ( (1.-dx)*A[i][j] + dx*A[i1][j] ) +
                    dy      * ( (1.-dx)*A[i][j1] + dx*A[i1][j1] );
    }
}
return B;

Etiquetas: C++, Interpolación, Procesamiento de Imágenes, Regresión Lineal, Splines