Incertidumbre de medición del centro de gravedad calculado a partir de la imagen de píxeles

Cada valor de píxel$S_i$en el detector de$\vec x_i$tiene algun error$N_i$: Los CCD, por ejemplo, tienen un ruido de fondo$N_\text{bkg}$de la electrónica de lectura, el ruido térmico y el fondo del cielo, además de un ruido fotónico poissoniano$N_S=\sqrt{S}$. En muchos casos, este ruido sigue razonablemente bien una distribución gaussiana. Después de restar el fondo, una medición de posición

$$\vec x = \frac{\sum{S_i \vec x_i}}{\sum{S_i}}$$

tiene la incertidumbre

$$\operatorname{std} \vec x = \left[\sum_i \left(\frac{\partial x}{\partial S_i} N_i\right)^2 \right]^{1/2} = \frac{\left[\sum_i \left( \sum_j\frac{(x_i - x_j) S_j}{\sum_k S_k} N_i \right)^2 \right]^{1/2}}{\sum_k S_k}$$

$$N_i^2 = N_\text{bkg}^2 + N_{S_i}^2$$

utilizando la propagación de errores gaussianos... Espero haber ejecutado las matemáticas correctamente. La forma en que lo he escrito es un poco extraña, pero muestra una propiedad importante de este método: si miras la suma sobre$j$, puede ver cómo el ruido de un píxel se pondera básicamente por la distancia del píxel al resultado del centroide. El mismo error en el valor de un píxel tiene más impacto si el píxel está más alejado del centroide.

Los mejores métodos explican la incertidumbre$N_i$en los valores de píxeles en la primera ecuación ya. Puede hacer esto introduciendo pesos adicionales en su centroide o yendo a modelos de ajuste, que es el método "usual" en astrometría, hasta donde he visto.

Este enfoque más complicado para medir posiciones utiliza la función de dispersión de puntos$P(\vec x_i - \vec x_\text{obj})$del instrumento en determinadas condiciones de observación. El modelo puede ser una aproximación, por ejemplo, una función de Moffat, o un modelo empírico construido a partir de estrellas brillantes no confundidas en la imagen, por ejemplo, una interpolación spline. Para fuentes puntuales, el ajuste típico de mínimos cuadrados de posición y flujo del modelo a una imagen produce fácilmente resultados cercanos al óptimo estadístico de estimación de parámetros en términos de incertidumbre. Con nuestra potencia informática moderna, la forma más fácil de obtener la incertidumbre para un modelo dado y el ruido de datos suele ser un algoritmo de arranque .

Por supuesto, los objetos extendidos requieren un poco más de trabajo en el modelo, por ejemplo, suposiciones sobre su forma, como lo ha indicado en su pregunta.

Puede echar un vistazo a este artículo de 2006 de Thomas et al. sobre algoritmos de centroide para sistemas astronómicos de óptica adaptativa (AO), que incluye una discusión detallada de las estimaciones de error para la posición del centroide usando diferentes algoritmos. El enfoque que describe en su respuesta corresponde a lo que llaman "centroide simple" (Sección 3); se refieren a un capítulo de libro de Rousset (1999) para el análisis detallado (que creo que incluye contribuciones de Poisson y ruido de lectura, por lo que no es idéntico a su resultado).

En términos más generales, el enfoque del "centro de gravedad" parece usarse en situaciones en las que se necesitan estimaciones rápidas y computacionalmente baratas, como en los sistemas AO (donde es necesario determinar un centroide estelar muchas veces por segundo), o como un primera conjetura cruda para proporcionar un punto de partida para un análisis más sofisticado. El análisis posterior a la observación de imágenes astronómicas generalmente usa enfoques más complejos/sofisticados, dependiendo de si la fuente es una estrella u otras fuentes puntuales versus una fuente extendida, si tiene un modelo preciso de la función de dispersión de puntos (que puede ser no circular ), eliminación de mezcla de fuentes vecinas, cuáles son las características de ruido de sus datos, etc.

En la práctica, supongo que la mayoría de estos análisis utilizan algún tipo de análisis no lineal de mínimos cuadrados o de máxima verosimilitud que implica ajustar un modelo a los datos. Los errores en los parámetros del modelo ajustado (incluida la posición del centroide) se pueden derivar de suposiciones simplistas sobre el paisaje de ajuste (p. ej., Levenberg-Marquardt y otros algoritmos de minimización basados ​​en gradientes a veces proporcionan matrices de covarianza basadas en el tratamiento de los parámetros locales).$\chi^{2}$paisaje como una parábola), de remuestreo bootstrap, o de enfoques de Markov Chain Monte Carlo. Esto se puede complementar ejecutando simulaciones del proceso de ajuste en imágenes artificiales de modelos simples de estrellas o galaxias, para obtener algunas estimaciones o correcciones casi empíricas para las incertidumbres del centroide (y otros parámetros).

No sé sobre "... generalmente determinado ..." pero en un trabajo similar he realizado (y he sido testigo) de que se realiza un ajuste de spline bidimensional en los datos de intensidad de píxeles para determinar el centroide de potencia a la resolución de subpíxeles. Como notó, hacemos la suposición razonable de que el objeto es casi esférico y tiene una densidad azimutal constante (es decir, la densidad puede variar con el radio pero no con el ángulo).

La incertidumbre (varianza) en este cálculo generalmente se calcula aplicando métodos estadísticos estándar a las variaciones observadas en la señal recibida en cada píxel, después de tener en cuenta la fluctuación en la posición de la línea de visión. (y, por supuesto, teniendo en cuenta el ruido electrónico, etc.). Esencialmente, se puede observar la variación en el pico calculado por spline en N fotogramas, o observar la variación en todos los píxeles, determinar su contribución al ajuste de spline y ponderar sus contribuciones en consecuencia.

Este interesante problema es claramente relevante también en astronomía. Dentro de las imágenes médicas, la determinación del centro de un objeto también entra en juego en varias áreas de aplicación especiales.

Planteé la pregunta aquí para informarme sobre qué fórmulas se están utilizando dentro de la astronomía: ¿Cómo se estima la incertidumbre del 'COG', en presencia de ruido de medición (aditivo)?

En 2002, derivamos una fórmula general para la varianza del centro de gravedad estimado, en presencia de ruido aditivo normalmente distribuido, asociado con cada valor de píxel. La referencia es la siguiente:

HC van Assen, M. Egmont-Petersen, JHC Reiber. "Localización precisa de objetos en imágenes de nivel de gris utilizando la medida del centro de gravedad: precisión frente a precisión", IEEE Transactions on Image Processing , vol. 11, núm. 12, págs. 1379-1384, 2002.

Primero daré la fórmula general, que solo asume que el centro de la ventana que rodea el objeto es$(0,0)$. El objeto no tiene que estar colocado en el centro para que esta fórmula se mantenga.

Defina el ruido de medición aditivo asociado con cada píxel como$\epsilon \sim U(0,\sigma^2)$, con$\sigma$siendo su desviación estándar.

Definir una (sub)imagen cuadrada${\cal W}$de dimensiones$(d+1) \times (d+1)$píxeles, con coordenadas,$x=-\frac{d}{2},\frac{d}{2}+1,\ldots,\frac{d}{2}$, y$y=-\frac{d}{2},\frac{d}{2}+1,\ldots,\frac{d}{2}$. Dejar${\bf w}_{x,y}$ser la verdadera intensidad de píxel de píxel$(x,y)$en ausencia de cualquier ruido, en${\cal W}$. Definir la imagen de señal más ruido$W$como:$w_{x,y}={\bf w}_{x,y} + \epsilon$. Defina el número total de píxeles en la (sub)imagen como$N=(d+1)^2$, que incluye la 0-ésima fila central y la 0-ésima columna de$W$.

los valores de$w_{x,y}$, con ($w_{x,y} \geq 0$), son los que realmente se están observando. cerca del centro$(x=0,y=0)$se ha localizado un objeto brillante de interés.

El centro de gravedad estimado$cog$de este objeto se calcula de la siguiente manera:

$$\widehat{cog}(x,y)=\left(\dfrac{\sum_{x,y} \; x\,w_{x,y}}{\sum_{x,y} \; w_{x,y}},\dfrac{\sum_{x,y} \; y\,w_{x,y}}{\sum_{x,y} \; w_{x,y}}\right)$$ donde $x$y $y$están ejecutando índices tales que cada píxel en $W$entra en el cálculo de una suma (sigma) exactamente una vez. $cog$es la medida que está influenciada por el ruido de medida.

Usando la regla delta dos veces sucesivas, derivamos una fórmula general aproximada para la varianza del cog, dado un nivel de ruido conocido.

Definir$x2$como:

$$x2 = \sum_{x} \; \sum_{y} \; x^2$$ y de manera similar $y2$como: $$y2 = \sum_{x} \; \sum_{y} \; y^2$$ Finalmente, el 'peso' promedio (intensidad de píxel promedio) viene dado por: $$\hat{\mu}_w = N^{-1} \; \sum_{x} \; \sum_{y} \; w_{x,y}$$ Dejar $\widehat{cog}(x)$denota el centro de gravedad estimado $x$-coordinar y $\widehat{cog}(y)$el centro de gravedad estimado $y$-coordinar.

Las estimaciones de varianza derivadas de MLE de$\widehat{cog}(x)$y$\widehat{cog}(y)$están:

$$\text{var}(\widehat{cog}(x))=\left(\frac{\sigma^2 \; x2}{\left[N \; \hat{\mu}_w \; \widehat{cog}(x)\right]^2} + \frac{\sigma^2}{N \; (\hat{\mu}_w)^2} \right) \; (\widehat{cog}(x))^2$$ y $$\text{var}(\widehat{cog}(y))=\left(\frac{\sigma^2 \; y2}{\left[N \; \hat{\mu}_w \; \widehat{cog}(y)\right]^2} + \frac{\sigma^2}{N \; (\hat{\mu}_w)^2} \right) \; (\widehat{cog}(y))^2$$ Resulta que cuando el verdadero engranaje es precisamente $(0,0)$, entonces se mantienen las fórmulas de límite (simplificadas) para la varianza de cog: $$\lim_{cog(x) \to 0} \; \; \lim_{cog(y) \to 0} \; \text{var}(\widehat{cog}(x)) = \frac{\sigma^2 \; x2}{(N \; \hat{\mu}_w)^2}$$ y $$\lim_{cog(x) \to 0} \; \; \lim_{cog(y) \to 0} \; \text{var}(\widehat{cog}(y)) = \frac{\sigma^2 \; y2}{(N \; \hat{\mu}_w)^2}$$ ya que el segundo término de varianza (aditivo) dentro de los paréntesis exteriores desaparece.

Mis simulaciones recientes muestran que arriba$95\%$de la varianza real resulta de nuestra fórmula de varianza definida, como se presenta aquí. Este resultado de simulación también es válido cuando el verdadero cogse desvía más de una coordenada de la posición central$(0,0)$.

El gráfico de simulación que muestra la precisión de la estimación de la varianza se agregará a esta respuesta, uno de los próximos días.

Ruido multiplicativo

El ruido de Poisson ocurre en las cámaras CCD, su influencia se puede observar especialmente en las transiciones de gradiente de la intensidad de la imagen desde áreas muy oscuras a áreas realmente brillantes. Se sabe que la magnitud del ruido de Poisson es proporcional a la intensidad de la señal. En este caso, el ruido multiplicativo está presente en la imagen observada:

$$w_{x,y}={\bf w}_{x,y} \, \, \cdot \alpha \, \cdot \, \epsilon$$ con $\alpha$siendo una constante de proporcionalidad y $\epsilon$el término de ruido distribuido normalmente (que se aproxima bien a una distribución de Poisson para $w_{i,j} >0$).

Realizando un sencillo$\ln(\cdot)$transformación,$lw_{i,j}=\ln(w_{i,j})$produce una imagen transformada$lw$perturbado por el ruido aditivo . Posteriormente, el estimador de varianza cog se puede aplicar a la$lw$-imagen.