¿Cómo calculo numéricamente una función a partir de su gradiente ruidoso?

8

Tengo el modelo . s(x,y)=x2+y2,0x1,0y1

En lugar de observar el modelo directamente, estoy observando las derivadas del modelo + algo de ruido (e):

 p(x,y)=sx+e,q(x,y)=sy+e

A partir de las mediciones de p (x, y y q (x, y) quiero estimar s (x). Digamos que sé que s (0,0) = 0.

Según el teorema del gradiente:  s(x,y)=(0,0)(x,y)[sx,sy]dr

independientemente de qué camino integramos.

Como un pequeño experimento (en Matlab) agregué ruido distribuido normal, N (0,1), a p = 2x y q = 2y. Luego integré primero a lo largo de x seguido de a lo largo de y: SXY. Luego integré primero a lo largo de y seguido de a lo largo de x: SYX.

Los resultados muestran que el teorema del gradiente no se cumple en este caso (debido al ruido):

S

SXY

SYX

Los errores de raíz cuadrática media relativos al modelo son:

ErmsXY =
    0.1125
ErmsYX =
    0.0920

¿Cómo puedo encontrar una mejor estimación (menos error RMS y más suave) de s de py q?

EDITAR:

Por lo que leí; El uso de la curva integral se conoce como integración local. También hay métodos de integración global en los que uno intenta elegir una S (x, y) que minimice:

 0101[|SxP|2+|SyQ|2]dxdy

Se supone que los métodos de integración global dan mejores resultados cuando el gradiente es ruidoso, pero ¿cómo hago esto en la práctica?

EDITAR 2:

Un enfoque que he usado es este:

primero presentamos operadores de derivación lineal: . sx=Dxs,sy=Dys

El resultado es el siguiente sistema de ecuaciones lineales:

 Dxs=p,Dys=q

Luego encuentre una solución de Mínimo error cuadrado para estas ecuaciones. Se supone que una solución LSE para estas ecuaciones es equivalente a minimizar la integral desde arriba. ¿Cómo se puede mostrar esto?

Los resultados son buenos: ingrese la descripción de la imagen aquí

El error RMS es aproximadamente 1/5 del de SXY y SYX y la solución también es más fluida.

Sin embargo, hay algunos inconvenientes en este enfoque:

  1. es difícil de implementar; debe usar las diferencias centrales y "aplanar" la matriz 2D s en un vector, etc.

  2. Las matrices de derivación son muy grandes y dispersas, por lo que pueden consumir mucha RAM.

Otro enfoque que parece potencialmente más simple de codificar, menos consumo de RAM y más rápido es usar FFT. En el espacio de Fourier, estos pdes se convierten en una ecuación algebraica. Esto se conoce como el algoritmo Frankot-Chellappa, pero desafortunadamente no lo tengo funcionando en mis datos de ejemplo.

Andy
fuente

Respuestas:

1

Puede filtrar el gradiente o el resultado, . Necesitaría conocer bastante bien las características de los gradientes verdaderos para saber cuál es el ancho de banda de frecuencia. En ese punto, podría diseñar un filtro de paso bajo que preservaría la señal pero eliminaría el ruido de frecuencia más alta.s

Jim Clay
fuente
Gracias Jim Entonces, por ejemplo, puedo tomar SXY y reemplazar cada valor SXY (xi, yj) por una suma ponderada sobre el valor y sus vecinos, donde los pesos pueden ser, por ejemplo, un gaussiano 2D?
Andy
Lo siento Jim Olvidé enfatizar que también quiero un pequeño error RMS en relación con el modelo. Edité mi pregunta para tener esto en cuenta. ¿El suavizado da un resultado más suave, pero no un error RMS más pequeño?
Andy
@Andy Sí, "una suma ponderada sobre el valor y sus vecinos" es una descripción bastante sucinta del filtrado, y un gaussiano 2D es una forma de filtro de paso bajo.
Jim Clay
Para @Andy error menor que estimaría el ancho de banda por FFT'ing múltiple "limpio" (sin ruido añadido) resultados, y ver donde existe la mayor atenuación de frecuencia es (estoy asumiendo que no todos son iguales). Diseñe un LPF con ese mismo rolloff (la "fdatool" de Matlab puede ayudar con esto) y luego use ese filtro. Debería mejorar su RMS. Todavía habrá un error, por supuesto, pero debería reducirse. s
Jim Clay
Gracias Jim Pero, ¿no hay forma de combinar los resultados de SXY y SYX para obtener un error RMS más pequeño?
Andy