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 o 1.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)

~~~

Ahora 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 debemos deshacer la matriz de Bayer para procesar de forma individual cada uno de los 4 canales RAW (R, G₁, G₂, B). En nuestro caso escogemos el canal R del RAW asumiendo que la respuesta de los demás será equivalente. Superado el efecto del filtro Bayer no hay motivo para pensar que los fotocaptores de un color puedan tener un rendimiento diferente a los demás (más sobre esto al final del artículo).

Por simplicidad y también por rigor, la corrección la vamos a hacer sin interpolar valores. Realizaremos un muestreo nearest neighbour que implicará tanto desechar algunos píxeles origen como duplicar otros. Dado que la selección es aleatoria no afecta a la precisión del cálculo posterior.

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
• La Relación S/N: cociente entre los anteriores

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 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++.

En la toma ISO25600 por ejemplo descartamos los primeros 3 parches, usando el dato en los restantes 74 del total de 77 parches:

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. En el eje X representamos los niveles de Señal en escala logarítmica respecto a la saturación (EV) y en el eje Y los valores de relación S/N igualmente en escala logarítmica (dB) (hacer clic para ver en alta resolución):

~~~

Ahora solo queda calcular el rango dinámico para cada ISO como el margen desde la saturación (0EV) hasta el nivel de Señal en que la relación S/N caiga a cierto valor mínimo. Para ello es preciso establecer un criterio umbral de relación S/N mínima aceptable, y aquí 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, es decir sin necesariamente pasar por ellos, las series de puntos correspondientes a cada valor ISO. 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 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):

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 ISO16000)".

La siguiente tabla resume el RD calculado para cada ISO y criterio, indicando además las muestras o puntos de la curva (parches) que participaron del total de 77 parches:

~~~

Vamos a terminar el ejercicio haciendo un estudio de sensibilidad del método descrito a tres variables que hemos elegido en el proceso:

• El canal RAW analizado (R, G₁, G₂, B).
• El nivel de negro del sensor usado para linealizar y normalizar los niveles RAW.
• El tamaño de los parches sobre los que hemos hecho las mediciones de ruido.

Sobre el canal RAW analizado, todo el ejercicio se ha hecho con el canal R, casualmente el que mayores niveles RAW obtuvo en los parches magenta. El resultado sería diferente usando los otros canales?. Aquí puede verse la superposición del cálculo para los 4 canales RAW (R, G₁, G₂, B), arrojando todos ellos curvas, y por tanto valores de RD, muy similares (hacer clic para ampliar):

Tiene mucha más influencia, y por eso debemos resaltarlo como punto crítico del método, el nivel de negro escogido (BLACK=256 en nuestro caso). El cálculo de RD a los ISOs más bajos puede fluctuar del orden de 0,1EV solo por usar el nivel de negro adyacente (BLACK=255) (hacer clic para ampliar):

Esta sensibilidad disminuye conforme aumentamos el ISO o el número de bits del sensor, al reducirse el error relativo que supone desplazar BLACK en una unidad. En conclusión es más complicado medir con precisión el RD de un sensor cuanto menor sea el ISO o sus archivos RAW se codifiquen con menos bits.

En el siguiente histograma de la captura ISO200 vemos que solo se tienen 58 valores tonales por canal para codificar toda la información de la escena. Mover en una sola unidad el nivel de negro tiene un gran impacto a la hora de situar correctamente los niveles de exposición respecto a la saturación, lo que trastoca todos los puntos de las gráficas:

El último parámetro de sensibilidad a evaluar es el tamaño de los parches considerados. Cuanto más grandes tendremos más muestras y por tanto mayor robustez estadística, pero también más peligro de interpretar como ruido efectos indeseables como la invasión de parches cercanos por mal alineamiento o por el desenfoque de la óptica, o directamente incorporar gradientes de luminosidad en la captura (viñeteos, reflejos,...).

Comparando parches de tamaños bastante dispares (separación entre parches de 100 píxeles vs 160 píxeles, suponiendo los primeros 4 veces más muestras) vemos que la diferencia es mínima, lo que valida las capturas usadas en este aspecto (hacer clic para ampliar):

~~~

Repositorio con el código R: GitHub.