Algoritmos de integración numérica y convolución de imágenes en C++

Implementación de Algoritmos Numéricos en C++

1. Procesamiento de Imágenes: Método mn_mascara3x3

A continuación se presenta la implementación de una función para aplicar una máscara de convolución 2D de 3x3 sobre una matriz de datos:

mn_mascara3x3 () { 
  // P51_Mascara2D
  if(F.dim1() < 1 || F.dim2() < 1 || m.dim1() != 3 || m.dim2() != 3) return Array2D< real >();

  Array2D< real > M(F.dim1(), F.dim2(), 0.);

  for(int i = 1; i < M.dim1() - 1; i++) {
    for(int j = 1; j < M.dim2() - 1; j++) {
      for(int k = -1; k < 2; k++) {
        for(int l = -1; l < 2; l++) {
          M[i][j] += m[k+1][l+1] * F[k+i][l+j]; 
        } 
      } 
    } 
  }

  for(int j = 1; j < M.dim2() - 1; j++) {
    for(int k = -1; k < 2; k++) {
      for(int l = -1; l < 2; l++) {
        M[0][j] += m[k+1][l+1] * F[max(k, 0)][l+j];
        M[M.dim1()-1][j] += m[k+1][l+1] * F[M.dim1()-1 + min(k, 0)][l+j]; 
      } 
    } 
  }

  for(int i = 1; i < M.dim1() - 1; i++) {
    for(int k = -1; k < 2; k++) {
      for(int l = -1; l < 2; l++) {
        M[i][0] += m[k+1][l+1] * F[k+i][max(l, 0)];
        M[i][M.dim2()-1] += m[k+1][l+1] * F[k+i][M.dim2()-1 + min(l, 0)]; 
      } 
    } 
  }

  for(int k = -1; k < 2; k++) {
    for(int l = -1; l < 2; l++) {
      M[M.dim1()-1][M.dim2()-1] += m[k+1][l+1] * F[M.dim1()-1 + min(k, 0)][M.dim2()-1 + min(l, 0)];
      M[0][M.dim2()-1] += m[k+1][l+1] * F[max(k, 0)][M.dim2()-1 + min(l, 0)];
      M[M.dim1()-1][0] += m[k+1][l+1] * F[M.dim1()-1 + min(k, 0)][max(l, 0)];
      M[0][0] += m[k+1][l+1] * F[max(k, 0)][max(l, 0)]; 
    } 
  }

  return M;
}

El main del programa lee una imagen en formato .png como una matriz. Tras esto, se realiza el cálculo de los bordes de la imagen usando el método mn_calculo_bordes (ya implementado). Dicho método llama a la función mn_mascara3x3 (implementada por el estudiante); este proceso se repite aumentando el número de iteraciones y, finalmente, devuelve la imagen procesada al usuario.

2. Integración Numérica: Método del Rectángulo

Este método se utiliza para aproximar la integral definida de una función, aplicado aquí a un modelo epidemiológico:

real rectangulo() { 
  // P52_COVID_rectangulo
  /// HACER ALUMNO
  real integral = 0;
  real h = (b - a) / N;
  real x0 = a + 0.5 * h;

  for(int k = 0; k < N; k++) {
    integral += f(x0 + h * k);
  }
  return h * integral;
}

El main del programa crea la función del avance epidemiológico. Posteriormente, se calculan los compartimentos en un bucle de hasta 30 iteraciones, se determina el pico de incidencia diaria y se muestran los resultados reales comparados con los obtenidos mediante el método rectángulo.

3. Integración Numérica: Método del Trapecio

Implementación de la regla del trapecio para mejorar la precisión de la integración:

real trapecio() { 
  // P53_COVID_rectangulo 
  /// HACER ALUMNO
  real integral = 0;
  real h = (b - a) / N;
  real fk = f(a);

  for(int k = 0; k < N; k++) {
    real fk1 = f(a + (k + 1) * h);
    integral += fk + fk1;
    fk = fk1;
  }
  return 0.5 * h * integral;
}

El main del programa define la función de avance epidemiológico y calcula los compartimentos en un bucle de 30 pasos. Se evalúa el pico de incidencia diaria y se comparan los resultados reales con los calculados mediante el método trapecio.

4. Regla de Simpson para Integración

La regla de Simpson ofrece una mayor exactitud al utilizar polinomios de segundo grado:

real simpson() { 
  // P54_Simpson
  real integral = 0; 
  real h = (b - a) / N; 
  real xk, xk1; 
  xk = a; 
  real fk = f(xk); 
  real S = (b - a) / (6. * N); 

  for(int k = 0; k < N; k++) {
    xk1 = a + (k + 1) * h; 
    real xm = (xk + xk1) / 2; 
    real fm = f(xm); 
    real fk1 = f(xk1); 
    integral += (fk + fk1 + 4. * fm);
    xk = xk1; 
    fk = fk1; 
  }
  return S * integral;
}

El main del programa gestiona la función epidemiológica, calculando los compartimentos y el pico de incidencia diaria para comparar los valores reales con los obtenidos por el método Simpson.

5. Integración Doble: Método Simpson 2D

Extensión de la regla de Simpson para funciones de dos variables:

real mn_simpson_2D() { 
  // P55_Simpson2D
  real integral = 0; 
  real hx = (b - a) / M; 
  real hy = (d - c) / N;

  for(int k = 0; k < M; k++) {
    real xk = a + k * hx; 
    real xk1 = xk + hx; 
    real xm = xk + 0.5 * hx;
    for(int n = 0; n < N; n++) {
      real yn = c + n * hy; 
      real yn1 = yn + hy; 
      real ym = yn + 0.5 * hy;
      real fxk = (F(xk, yn) + 4 * F(xk, ym) + F(xk, yn1));
      real fxm = (F(xm, yn) + 4 * F(xm, ym) + F(xm, yn1));
      real fxk1 = (F(xk1, yn) + 4 * F(xk1, ym) + F(xk1, yn1));
      integral += (fxk + 4. * fxm + fxk1);
    }
  }
  return integral * hy * hx / 36.;
}

El main del programa define tres funciones distintas y establece el número de iteraciones para calcular la integral de cada una. Se muestra en pantalla el resultado real esperado y el obtenido mediante el método mn_simpson_2D.

6. Cuadratura de Gauss-Legendre 2D

Implementación de la cuadratura de Gauss para una integración numérica altamente eficiente:

real mn_cuadratura() { 
  // P56_Legendre2D
  /// HACER ALUMNO
  real integral = 0;
  for(int k = 0; k < N; k++) {
    real x2 = 0.5 * ((b - a) * x[N-1][k] + b + a);
    for(int m = 0; m < N; m++) {
      real y2 = 0.5 * ((d - c) * x[N-1][m] + c + d);
      integral += f(x2, y2) * w[N-1][k] * w[N-1][m];
    }
  }
  return (b - a) * (d - c) * 0.25 * integral;
}

El main del programa crea cinco funciones y declara los intervalos para los ejes x e y, junto con las matrices de pesos y puntos de la cuadratura de Gauss. Se cargan los datos mediante el método mn_leer_matriz y se comparan los resultados reales con diversas pruebas ejecutadas mediante el método mn_cuadratura.