---
title: ''
output:
  pdf_document: default
  html_notebook: default
  html_document: default
---

## A2 -- El modelo multiplicativo y distribuciones de intensidad

Según el modelo multiplicativo elemental (para áreas sin textura), la intensidad se puede modelar con una distribución gama de media $\sigma^2$ y parámetro de forma $L$: $Z = \sigma^2 Y\sim \sigma^2 \Gamma(1,L)$.
Nos interesa $\sigma^2$, mientras que $Y$ interfiere y, por ende, se considera ruido.

Diversos factores pueden hacer que $\sigma^2$ no se mantenga constante para todos los píxeles del área observada, aún cuando se trate de un mismo tipo de blanco de la escena (blanco homogéneo).
En ese caso, se puede modelar $\sigma^2$ como otra variable aleatoria, independiente de $Y$; basta que la variable aleatoria tenga soporte positivo (backscatter negativo no tiene sentido).

Se han propuesto diversos modelos para el backscatter $X$.
Uno de los primeros, que tiene buen fundamento físico, es una distribución gama.
Este modelo surge haciendo que el número de dispersores elementales $N$ siga una distribución binomial negativa.
Si suponemos que $X$ sigue una distribución gama, y manteniendo el modelo gama para el speckle $Y$, y suponiendo que $X$ e $Y$ son independientes, llegamos a una distribución $K$ [@yueh89].
Esta distribución es útil para modelar blancos con textura, esto es, con más variabilidad que la que el modelo gama puro puede explicar.
Tiene como desventaja que la densidad, y por ende todo su tratamiento analítico, involucra la función de Bessel modificada $K$ (de ahí su denominación).
Esta función no tiene forma explícita (es una integral), y aún buenas implementaciones padecen de problemas numéricos.

Otro modelo para el backscatter es la distribución recíproca de gama.
La distribución resultante para el retorno se denomina $\mathcal G_I^0$.
Fue propuesta en 1997 [@frery96], y se la considera un ``Modelo Universal'' para datos SAR.
Su densidad es
$$
f_Z(z;\alpha,\gamma,L) = \frac{L^L \Gamma(L-\alpha)}{\gamma^\alpha \Gamma(-\alpha) \Gamma(L)}
  \frac{z^{L-1}}{(\gamma+Lz)^{L-\alpha}},
$$
en que $\alpha<0$ es una medida de rugosidad, $\gamma>0$ es el parámetro de escala, y $L\geq1$ es el número de looks.
Esta densidad no involucra funciones especiales, y es capaz de modelar una buena variedad de blancos.

Una implementación de esta densidad es
```{r}
dGI0 <- function(z, p_alpha, p_gamma, p_Looks) {
 return( 
  (p_Looks^p_Looks * gamma(p_Looks-p_alpha) * z^(p_Looks-1) ) / 
  (p_gamma^p_alpha * gamma(-p_alpha) * gamma(p_Looks) * (p_gamma + z*p_Looks)^(p_Looks-p_alpha))
       ) 
}
```

El momento de orden $r$ de una variable aleatoria que obedece esta distribución está dado por
$$
E(Z^r) =\Big(\frac{\gamma}{L}\Big)^r\frac{\Gamma ( -\alpha-r )}{ \Gamma (-\alpha) }\cdot  
\frac{\Gamma (L+r )}{\Gamma (L)},
$$
si $-\alpha>r$, e infinito caso contrario.
Con esta condición podemos elegir el valor de $\gamma$ que impone esperanza unitaria, digamos $\gamma^*$, para las ilustraciones que siguen:
$$
\gamma^*=-\alpha-1.
$$
Veamos primero tres densidades de modelos para datos extremamente heterogéneos, esto es, con $\alpha=-1.5$, primero en escala lineal

```{r DensityPlotGI0ExtHet_Linear}
plot(c(0,5), c(0,3), type="n", axes = FALSE, xlab = "Intensidad", ylab = "Densidades GI0 de Media Unitaria",
     main = "Densidades del Retorno Extremamente Heterogéneo - Escala Lineal")
axis(1)
axis(2)
curve(dGI0(x, -1.5, 0.5, 1), col=colores[3], lwd=3, n = 500, from=0, to=5, add=TRUE)
curve(dGI0(x, -1.5, 0.5, 3), col=colores[2], lwd=3, n = 500, add=TRUE)
curve(dGI0(x, -1.5, 0.5, 8), col=colores[1], lwd=3, n = 500, add=TRUE)
legend("topright", lwd = c(3,3,3), col = colores[c(3,2,1)], legend = c("L = 1", "L = 3", "L = 5"))
```

y en escala semilogarítmica
```{r DensityPlotGI0ExtHet_Semilog}
plot(c(0,5), c(10^(-2),10), type="n", axes = FALSE, log = "y", xlab = "Intensidad", ylab = "Densidades GI0 de Media Unitaria", main = "Densidades del Retorno Extremamente Heterogéneo - Escala Semilogarítmica")
axis(1)
axis(2)
curve(dGI0(x, -1.5, 0.5, 1), col=colores[3], lwd=3, n = 500, from=0, to=5, add=TRUE)
curve(dGI0(x, -1.5, 0.5, 3), col=colores[2], lwd=3, n = 500, add=TRUE)
curve(dGI0(x, -1.5, 0.5, 8), col=colores[1], lwd=3, n = 500, add=TRUE)
legend("topright", lwd = c(3,3,3), col = colores[c(3,2,1)], legend = c("L = 1", "L = 3", "L = 5"))
```

Veamos ahora modelos análogos para un caso de textura intermedia, esto es, para $\alpha=-5$
```{r DensityPlotGI0Het_Linear}
plot(c(0,3), c(0,1.5), type="n", axes = FALSE, xlab = "Intensidad", ylab = "Densidades GI0 de Media Unitaria",
     main = "Densidades del Retorno Heterogéneo - Escala Lineal")
axis(1)
axis(2)
curve(dGI0(x, -3, 2, 1), col=colores[3], lwd=3, n = 500, from=0, to=5, add=TRUE)
curve(dGI0(x, -3, 2, 3), col=colores[2], lwd=3, n = 500, add=TRUE)
curve(dGI0(x, -3, 2, 8), col=colores[1], lwd=3, n = 500, add=TRUE)
legend("topright", lwd = c(3,3,3), col = colores[c(3,2,1)], legend = c("L = 1", "L = 3", "L = 8"))
```

y en escala semilogarítmica
```{r DensityPlotGI0Het_Semilog}
plot(c(0,3), c(10^(-2),1.5), type="n", axes = FALSE, log = "y", xlab = "Intensidad", ylab = "Densidades GI0 de Media Unitaria", main = "Densidades del Retorno Heterogéneo - Escala Semilogarítmica")
axis(1)
axis(2)
curve(dGI0(x, -3, 2, 1), col=colores[3], lwd=3, n = 500, from=0, to=5, add=TRUE)
curve(dGI0(x, -3, 2, 3), col=colores[2], lwd=3, n = 500, add=TRUE)
curve(dGI0(x, -3, 2, 8), col=colores[1], lwd=3, n = 500, add=TRUE)
legend("topright", lwd = c(3,3,3), col = colores[c(3,2,1)], legend = c("L = 1", "L = 3", "L = 8"))
```
