Rango dinámico de un sensor de imagen con R

En este artículo vamos a desgranar paso a paso la metodología clásica para calcular el rango dinámico de un sensor digital de imagen, basándonos en mediciones precisas de ruido hechas sobre capturas RAW de parches de color uniforme. Se define el rango dinámico (RD) de un sensor como el intervalo de niveles de luminosidad (típicamente medido en pasos EV) en los que dicho sensor es capaz de capturar información útil. El RD está limitado en su extremo superior por la propia saturación del sensor y en el extremo inferior por el nivel de ruido máximo aceptable.

Para ello vamos a echar mano de una serie de capturas RAW amablemente compartidas por Hugo Rodríguez procedentes de su Olympus OM-1, una para cada ISO de la cámara. La carta de parches de color usada por Hugo es la que Bill Claff (Photons to Photos) propone en aquí para calcular el RD, y consiste en proyectar sobre el monitor una cuadrícula de 11x7 parches de color magenta:

El color magenta no es caprichoso, el astuto de mi tocayo lo ha escogido sabedor de que el canal G siempre obtiene mayores niveles en el RAW. Con parches magenta (pobres en verdes), los niveles RAW del canal G se hacen más parejos a los R y B proporcionando un set de muestras más uniforme. La toma a ISO12800 tenía este aspecto:

Que la toma esté desenfocada no solo no es un error sino que es imprescindible para no capturar ningún tipo de textura (en este caso de los píxeles del monitor) que sería interpretada como ruido adicional, falseando todos los cálculos de relación S/N.

También es importante que no exista ningún gradiente de iluminación causado por viñeteo de la óptica o por reflejos en la pantalla, el cual igualmente incrementaría la medición de ruido. Siendo el área útil que aprovecharemos en cada parche pequeña en relación a la imagen total, estos gradientes sería raro que supusieran un problema.

~~~

El siguiente paso es disponer de los niveles de negro y de saturación RAW del sensor bajo estudio. En cuanto al primero, seguramente por adecuación a la electrónica del conversor A/D y a la aritmética digital con enteros, los fabricantes sitúan deliberadamente el nivel de negro (valor medio que tendría un RAW obtenido haciendo una foto con el objetivo tapado) en las potencias de 2:

• Sensores de 10 bits (1.024 niveles): nivel de negro 64 o 128
• Sensores de 12 bits (4.096 niveles): nivel de negro 256
• Sensores de 14 bits (16.384 niveles): nivel de negro 512 o1.024

Con el comando dcraw -v -D -t 0 *.DNG extraemos en formato TIFF los datos de todos los archivos RAW. Siendo la Olympus OM-1 una cámara de 12 bits lo esperable sería encontrar el nivel de negro en el valor 256, lo que comprobamos con el histograma RAW de la toma ISO800:

Comentamos que el rango dinámico del sensor está limitado en el extremo derecho por la saturación, por lo que necesitamos de forma imperativa conocer con precisión en qué valor satura para referir todos los demás niveles RAW a dicha saturación. Podría pensarse que los sensores saturan en el fondo de su escala de bits pero no siempre es así, aunque éste sí lo hace:

Que la toma ISO25600 sature en 4.095 no garantiza que las de ISOs inferiores lo hagan, pero ningún RAW por debajo de éste tiene píxeles saturados así que asumiremos este valor de saturación para todos. El error en cualquier caso será mínimo.

Los niveles RAW los leemos como valores enteros y los normalizamos al rango 0..1 en coma flotante, permitiendo valores normalizados negativos:

 BLACK=256  # Olympus OM-1 black level

 SAT=4095  # Olympus OM-1 sat level

 img=readTIFF("iso25600.tiff", as.is=TRUE)

 img=img-BLACK

 img=img/(SAT-BLACK)

~~~

En el siguiente paso vamos a realizar la corrección de perspectiva de la captura original. Esta corrección no es imprescindible pero tener los parches alineados facilita mucho la tarea de leer los datos en ellos.

Previamente a la corrección deberemos quedarnos de forma individual con cada uno de los 4 canales RAW (R, G1, G2, B) deshaciendo la matriz de Bayer. En nuestro caso escogemos el canal R del RAW asumiendo que la respuesta de los demás será idéntica. Superado el efecto del filtro Bayer no hay motivo para pensar que los fotocaptores de un color deban tener distintas prestaciones a los demás.

Por facilidad y también por exactitud, la corrección la vamos a hacer sin interpolar valores en ningún momento. Realizaremos un muestreo nearest neighbour que implicará tanto desechar algunos píxeles origen como duplicar otros. Como la selección de estos efectos va a ser totalmente aleatoria no afecta en nada a la precisión del cálculo que haremos a continuación.

Para la corrección de perspectiva se ha utilizado el método visto en 'Transformación trapezoidal de imágenes con R (I). Algoritmo', pero optimizado en C++ con Rcpp reduciendo el tiempo de cálculo por un factor 240!!!

~~~

El siguiente paso es el meollo de todo el ejercicio donde recorremos individualmente los parches calculando en cada uno:

• El nivel de Señal útil: promedio de los valores RAW
• El nivel de Ruido: desviación estándar de los valores RAW
• El cociente entre ambos será la Relación S/N

Desechamos los parches en los que el nivel de Señal resulte un valor negativo, la relación S/N sea inferior a -10dB (zonas tan ruidosas no nos aportan nada) o el número de píxeles cercanos a saturar (>90% en escala lineal) llegue al 1% del total. Este análisis de los parches también lo hemos acelerado en C++.

A partir de las parejas de Señal y relación S/N obtenidas, cada una proveniente de un parche, ya podemos intuir la forma de las Curvas de relación S/N en el scatterplot (hacer clic para ver en alta resolución):

~~~

Ahora solo queda calcular el rango dinámico para cada ISO y para ello es preciso establecer un criterio umbral de relación S/N mínima aceptable. Usaremos los dos baremos habituales:

• RD fotográfico (práctico): pasos desde saturación hasta S/N=12dB
• RD ingenieril (teórico): pasos desde saturación hasta S/N=0dB

Para poder automatizar este cálculo, parametrizamos curvas splines cúbicas que aproximen suavemente las series de puntos correspondientes a cada valor ISO, sin necesariamente pasar por ellos. Usamos la relación S/N como variable independiente y el nivel de Señal como variable dependiente, ya que los umbrales de 12dB y 0dB son los valores de relación S/N de entrada para los que queremos obtener sus correspondientes niveles de Señal (el RD).

Las splines serán más efectivas si las calculamos en el dominio logarítmico para ambos ejes, tanto el del nivel de Señal (EV) como el de relación S/N (dB). La transformación genera curvas más suaves que evitan los saltos bruscos de la segunda derivada que penalizan el cálculo de splines.

Sobre la misma gráfica se marcan los dos criterios de rango dinámico y los valores obtenidos de rango dinámico como intersección entre las curvas y cada uno de los criterios umbral (hacer clic para ver en alta resolución):

Si nos fijamos un poco, veremos que las curvas hasta ISO800 tienen un comportamiento parejo, el cual cambia en ISOs superiores a ISO800 donde mejora el rendimiento en las sombras profundas. Esto puede constatarse en la gráfica de RD para esta cámara de Bill Claff que precisamente presenta un salto en el RD tras ISO800, ver aquí. En sus propias palabras: "(...) the Olympus System OM-1 looks to have High Conversion Gain (HCG) starting at ISO1000 (and Noise Reduction starting at ISO 16000)".

Todo el ejercicio se ha hecho con el canal R, casualmente el que mayores niveles RAW obtuvo en los parches magenta, pero aquí puede verse la superposición de los 4 canales RAW arrojando todos curvas, y por tanto valores de RD, muy parecidos.

~~~

Repositorio con el código R: GitHub.