Curso: Procesamiento Digital de Imágenes Otoño 2019
Dr. G. Urcid - INAOE 1
Restauración de imágenes
mejoramiento de imágenes(i h t)
• métodos heurísticos• procesos subjetivos
(image enhancement) • orientados al sistemade percepción visualhumana
restauración de imágenes• criterios de optimalidad• procesos objetivosg
(image restoration)p j
• no orientados al sistemade percepción visualhumana
Modelado de degradacionesConsiste en considerar la presencia de un sistema lineal y
espacialmente invariante con o sin ruido (aditivo) quecausa una degradación en la imagen.
( , ) ( , ) ( , ) ( , )g x y f x y h x y x yη= ∗ +
( , ) H[ ( , )] ( , )g x y f x y x yη= +
?
Suposiciones: H es lineal y espacialmente invariante
1 2 1 2H[ ] H[ ] H[ ]H[ ( , )] ( , )
af bf a f b ff x y g x yα β α β+ = +− − = − −
( , ) 0x yη =
Curso: Procesamiento Digital de Imágenes Otoño 2019
Dr. G. Urcid - INAOE 2
Modelado de degradacionesCaracterización del sistema que causa la degradación:
respuesta al impulso
( , , , ) ( , ; , ) H[ ( , )]h x y h x y x yα β α β δ α β= = − −función de punto extendido (óptica)
Volvemos al modelo general,
imagendegradada
imagenoriginal
función dedegradación
ruido(aditivo)
dominio espacial ( , ) ( , ) ( , ) ( , )g x y f x y h x y x yη= ∗ +g g g ( )
dominio espectral ( , ) ( , ) ( , ) ( , )G u v F u v H u v N u v= +
dominio espacial ( , ) ( , ) ( , ) ( , )g x y f x y h x y x yη= ∗ +
imagendegradada
imagenoriginal
función dedegradación ruido
( , ) ( , ) ( , ) ( , )g y f y y yηˆ ( , ) [ ( , )]f x y R g x y=
imagenrestaurada
transformacioneso filtros de
restauración
aditivo
Curso: Procesamiento Digital de Imágenes Otoño 2019
Dr. G. Urcid - INAOE 3
( , ) H[ ( , )] ( , ) ; 0 ( , ) H[ ( , )]g x y f x y x y g x y f x yη η= + = ⇒ =
1 2 1 2H[ ( , ) ( , )] H[ ( , )] H[ ( , )]af x y bf x y a f x y b f x y+ = +
1a b= = ⇒
linearidad
1 2 1 2
1H[ ( , ) ( , )] H[ ( , )] H[ ( , )]a b
f x y f x y f x y f x y= = ⇒
+ = + aditividad
2 1 10, 0 y =0 H[ ( , )] H[ ( , )]a b f af x y a f x y≠ = ⇒ = homogeneidad
H[ ( , )] ( , ) ; ,f x y g x yα β α β α β− − = − − ∈R invariante a la posición,invariante espacialmente
( , ) ( , ) ( , )f x y f x y d dα β δ α β α β∞ ∞
−∞ −∞= − −∫ ∫
propiedad de “cribado” o “colado” (sifting) del impulso; valor puntual de f
seguimos suponiendo que
( , ) [ ( , ) ( , ) ]
[ ( ) ( )]
g x y H f x y d d
H f x y d d
α β δ α β α β
α β δ α β α β
∞ ∞
−∞ −∞
∞ ∞
= − −
= − −
∫ ∫∫ ∫
( , ) 0 ( , ) H[ ( , )]x y g x y f x yη = ⇒ =
[ ( , ) ( , )]
( , ) [ ( , )]
( , ) ( , ; , )]
H f x y d d
f H x y d d
f h x y d d
α β δ α β α β
α β δ α β α β
α β α β α β
−∞ −∞
∞ ∞
−∞ −∞
∞ ∞
−∞ −∞
= − −
=
∫ ∫∫ ∫∫ ∫
respuesta de H al( , , , ) ( , ; , ) H[ ( , )]h x y h x y x yα β α β δ α β= = − −
respuesta de H alimpulso unitario
función de punto extendido, PSF(point spread function)
respuesta a una fuente de luz puntual(point light source response)
Curso: Procesamiento Digital de Imágenes Otoño 2019
Dr. G. Urcid - INAOE 4
fuentes de luzpuntuales
funciones depunto extendido
PSF’s
( , ) ( , ) ( , ; , )]g x y f h x y d dα β α β α β∞ ∞
−∞ −∞= ∫ ∫
integral de Fredholmde primera especie;
integral desuperposición
Concepto fundamental en la teoría de sistemas lineales;cualquier sistema lineal H queda completamentecualquier sistema lineal H queda completamente
caracterizado por su respuesta al impulso
H[ ( , )] ( , )
( , ) ( , ) ( , )]
( , )* ( , )
x y h x y
g x y f h x y d d
f x y h x y
δ α β α β
α β α β α β∞ ∞
−∞ −∞
− − = − − ⇒
= − −
=∫ ∫ integral de convolución
producto de convolución
La respuesta de un sistema lineal espacialmente invariante H queda determinada por la convolución de la entrada con la respuesta al impulso
Curso: Procesamiento Digital de Imágenes Otoño 2019
Dr. G. Urcid - INAOE 5
En presencia de ruido, ( , ) 0x yη ≠ ⇒
( , ) ( , ) ( , ; , )] ( , )g x y f h x y d d x yα β α β α β η∞ ∞
−∞ −∞= +∫ ∫
Lineal
( , ) ( , ) ( , )] ( , )g x y f h x y d d x yα β α β α β η∞ ∞
−∞ −∞= − − +∫ ∫
Lineal Invariante
imagendegradada
imagenoriginal
función dedegradación
ruido(aditivo)
dominio espacial ( , ) ( , ) ( , ) ( , )g x y f x y h x y x yη= ∗ +degradada original degradación (aditivo)
dominio espectral ( , ) ( , ) ( , ) ( , )G u v F u v H u v N u v= +
Modelado de degradacionesDiversas clases de degradación admiten un modelo o representación
como sistemas lineales espacialmente invariantes (2D)
• degradación caracterizada por convolución• filtros de restauración para revertir el proceso de degradación• deconvolución de imágenes o filtros de “deconvolución”
Estimación de la función de degradación h(x,y) o H(u,v)
• por observación de la imagen d l ió i• por observación de la imagen• por experimentación• por modelado matemático
deconvolución ciega(blind deconvolution)
rara vez se conoce completamente la función degradadora
Curso: Procesamiento Digital de Imágenes Otoño 2019
Dr. G. Urcid - INAOE 6
( , )H u v =2 2 5/ 6exp[ ( ) ]k u v− +
“Borrosidad” por movimiento lineal uniforme
0 0( ), ( )x t y t componentes de velocidad en ambas direcciones
Hipótesis: abertura y cierre del diafragma es instantáneo, el proceso de formación de la imagen es perfecto y T es el tiempo de exposición.
0 00
2 ( )
2 ( )
( , ) [ ( ), ( )]
( , ) ( , )
[ ( ) ( )]
T
i ux vy
T i ux vy
g x y f x x t y y t dt
G u v g x y e dxdy
f x x t y y t dt e dxdy
π
π
∞ ∞ − +
−∞ −∞
∞ ∞ − +
= − − ⇒
=
⎡ ⎤
∫∫ ∫
∫ ∫ ∫ ( )0 00
2 ( )0 00
[ ( ), ( )]
[ ( ), ( )]
y
T i ux vy
f x x t y y t dt e dxdy
f x x t y y t e dxdy dtπ
−∞ −∞
∞ ∞ − +
−∞ −∞
⎡ ⎤= − −⎢ ⎥⎣ ⎦⎡ ⎤= − −⎢ ⎥⎣ ⎦
∫ ∫ ∫
∫ ∫ ∫
Curso: Procesamiento Digital de Imágenes Otoño 2019
Dr. G. Urcid - INAOE 7
0 0
0 0
2 [ ( ) ( )]
0
2 [ ( ) ( )]
0
( , ) ( , )
( , ) ( , ) ( , )
T i ux t vy t
T i ux t vy t
G u v F u v e dt
F u v e dt F u v H u v
π
π
− +
− +
=
= =
∫∫
0 02 [ ( ) ( )]
0( , )
T i ux t vy tH u v e dtπ− +⇒ = ∫
0 0( ) / , ( ) /x t at T y t bt T= = ⇒movimiento lineal uniforme
a, b son los desplazamientos máximos
2 [ ] /
0( , )
T i ua vb t TH u v e dtπ− +⇒ = ∫
cálculo de la integral2 / 2 / 2 0/
2 /
0
2
02 / 2 / 2 /
i t T i T T i TT i t T
i ii i
Te e ee dti T i T i T
T T e eT
π ξα π ξα π ξαπ ξα
π ξα π ξαπ ξα π ξα
π ξα π ξα π ξα
− − −−
−− −
= = −− − −
⎛ ⎞−⎜ ⎟
∫
2 2 2
sin( )2
i i
i ii i
e Tei i i
T e e e T ei
π ξα π ξα
π ξα π ξαπξα πξα
π ξα π ξα π ξα
πξαπξα πξα
−− −
= − = ⎜ ⎟⎝ ⎠
⎛ ⎞−= =⎜ ⎟
⎝ ⎠
2 / sinc( )T i t T ie dt T eπ ξα πξαξα− −⇒ =∫0 sinc( )e dt T eξα⇒ =∫
( )( , ) sinc( ) i ua vbH u v T ua vb e π− +⇒ = +
Curso: Procesamiento Digital de Imágenes Otoño 2019
Dr. G. Urcid - INAOE 8
( )( , ) sinc( ) ; , , 0i ua vbH u v T ua vb e a b Tπ− += + ∈ >R
Filtrado inversoEl filtrado inverso directo consiste en determinar laestimación espectral o espectro estimado de f(x,y).
( , )ˆ ( , )( , )
G u vF u vH u v
=estimaciónespectral función de degradación
estimada o modelada
Equivalentemente,
( , )ˆ ( , ) ( , ) N u vF u v F u v= + • ¿? desconocemos N(u,v)! i |H( )| 0 l 2d té i( , ) ( , )
( , )u v u v
H u v • ¡! si |H(u,v)| ≈ 0 el 2do términodominará la estimación
En general: DIF (direct inverse filtering) no se recomienda pues esinestable si H(u,v) tiene ceros o valores próximos a cero.
Curso: Procesamiento Digital de Imágenes Otoño 2019
Dr. G. Urcid - INAOE 9
Butterworth LPde orden n = 10
Filtrado de Wiener¿Cómo mejorar el desempeño del filtrado inverso directo?
• filtrado de Wiener (error cuadrático medio mínimo)• filtrado restringido (en el sentido de cuadrados mínimos)• filtrado de la media geométrica (generalización del filtro de Wiener)
Esta clase de filtros incorpora la función de degradación estimadad l ét d i d idcon uno de los métodos mencionados y considera
las características estadísticas del ruido.
Curso: Procesamiento Digital de Imágenes Otoño 2019
Dr. G. Urcid - INAOE 10
Consideraciones:
• la imagen y el ruido se consideran procesos aleatorios
• la determinación de una estimación de la imagen no corrompida( i i l) f li tili d did l d áti di(original) f se realiza utilizando como medida el error cuadrático medio(MSE – Mean Square Error) y minimizarlo
• la imagen y el ruido no están correlacionados
2 2ˆ{( ) } ; { } valor esperadoe E f f E= − •
• la imagen o el ruido tienen media cero
• los niveles de gris en la imagen estimada son una función linealde los niveles de gris en la imagen degradada g
Dominio espectral
*
2
( , ) ( , )ˆ ( , ) ( , )( , ) ( , ) ( , )
f
f
H u v S u vF u v G u v
S u v H u v S u vη
⎡ ⎤= ⎢ ⎥
+⎢ ⎥⎣ ⎦*
2
2
2
( , ) ( , )( , ) ( , ) / ( , )
( , )1 ( , )( , ) ( , ) ( , ) / ( , )
f
f
H u v G u vH u v S u v S u v
H u vG u v
H u v H u v S u v S u v
η
⎡ ⎤= ⎢ ⎥
+⎢ ⎥⎣ ⎦⎡ ⎤
= ⎢ ⎥+⎢ ⎥⎣ ⎦( , ) ( , ) ( , ) / ( , )fH u v S u v S u vη+⎢ ⎥⎣ ⎦
espectro depotencia del ruido
espectro depotencia de la
imagen original
Curso: Procesamiento Digital de Imágenes Otoño 2019
Dr. G. Urcid - INAOE 11
Otros nombres para el filtro de Wiener,
• filtro de error cuadrático medio mínimo (min MSE filter)
• filtro de error de cuadrados mínimos (LSE filter)
2
2
( , ) ( , )
( , ) ( , )f
S u v N u v
S u v F u v
η =
=
espectro de potencia del ruido
espectro de potencia de la imagen original
En el dominio espacial, 1ˆ ˆ( , ) { ( , )}Ff x y T F u v−=
id ( ) 0 ( ) 0S bti l DIF• ruido cero, ( , ) 0 ( , ) 0x y S u vηη = ⇒ = y se obtiene el DIF
• el denominador se anulapara las mismas u, v ( , ) 0 ( , )H u v S u vη⇔ = =
• ruido blanco (espectral)2( , )S N u vη⇒ = , sin embargo
( , )fS u v prácticamente es desconocido
• cuando ambos espectros de potencia, del ruido y la imagen originalNO SE CONOCEN la estimación del espectro se aproxima medianteNO SE CONOCEN, la estimación del espectro se aproxima mediante
2
2
( , )1ˆ ( , ) ( , ) ; 0 (const.)( , ) ( , )
H u vF u v G u v K
H u v H u v K
⎡ ⎤= >⎢ ⎥
+⎢ ⎥⎣ ⎦
primera aproximación delcociente de espectros de potencia
Curso: Procesamiento Digital de Imágenes Otoño 2019
Dr. G. Urcid - INAOE 12
Filtrado de Wiener
2 2 5/ 6( , ) exp[ ( ) ]H u v k u v= − +
RuidoGaussiano
2
0650
μ
σ
=
=
2
065
μ
σ
=
=
2
00.0065
μ
σ
=
<
Curso: Procesamiento Digital de Imágenes Otoño 2019
Dr. G. Urcid - INAOE 13
Filtrado restringido
• este método requiere solamente del conocimiento de la media yla varianza del ruido; pueden calcularse de la imagen degradada.
• el filtro de Wiener se basa en minimizar un criterio estadísticoy por tanto es óptimo en un sentido promedio.
• el algoritmo CLS (Constrained Least Squares) proporciona unresultado óptimo para cada imagen al que se aplique.
• recordar que estos criterios de optimalidad son satisfactoriosdesde un punto de vista teórico y no están relacionados con ladinámica de la percepción visual (humana).
( , ) ( , ) ( , ) ( , )g x y h x y f x y x yη= ∗ +
; dim( ) MN MN= ×g = Hf + n H; dim( )dim( ) 1 dim( ) dim( )
MN MNMN
= ×= × = =
g = Hf + n Hg f n
• manejar directamente los problemas de restauración usandomatrices y vectores no es trivial debido a sus dimensiones,y en particular la matriz H es sensitiva al ruido,
• sin embargo, la expresión matricial es útil para proponeresquemas de restauración
Curso: Procesamiento Digital de Imágenes Otoño 2019
Dr. G. Urcid - INAOE 14
Determinar el valor mínimode la función C sujeta a larestricción
1 1 22
0 0( , ) ( , )
M N
x yC x y f x y
− −
= =
⎡ ⎤= ∇⎣ ⎦∑∑
2 2 2ˆ donden
Tkx= =∑
2g - Hf = n x x x
1k=
Solución al problema de optimización en el dominio espectral
*
2( , )ˆ ( , ) ( , )
( , ) ( , )H u vF u v G u v
H u v P u vγ
⎡ ⎤= ⎢ ⎥
+⎢ ⎥⎣ ⎦
• γ es un parámetro ajustable para lograr la restricción dada
• P(u,v) es la transformada de Fourier del operador Laplaciano
0 1 0( , ) { ( , )} y ( , ) 1 4 1
0 1 0FP u v T p x y p x y
−⎡ ⎤⎢ ⎥= = − −⎢ ⎥⎢ ⎥−⎣ ⎦
Curso: Procesamiento Digital de Imágenes Otoño 2019
Dr. G. Urcid - INAOE 15
Filtrado deWiener
Filtrado CLS
Cómputo de γ por iteración
1. Dar un valor inicial a γ y un factor de precisión ‘a’2. Calcular la norma euclídea del vector residual r(γ),
2
se desea ajustar γ de modo que(Nota: ϕ(γ) = || r ||2 es una función monótona creciente)
3. Parar si se cumple la condición de ajuste en 2 ; de otra manera
2 ˆ=2
r g - Hf
2 a= ±2r n
p j ;recalcular el mismo paso (2), después de incrementar γ si
, o decrementar γ si . Usar elnuevo valor de γ y recalcular el estimado óptimo
2 a< −2r n 2 a> +2r nˆ ( , )F u v
Curso: Procesamiento Digital de Imágenes Otoño 2019
Dr. G. Urcid - INAOE 16
• determinación de || r ||2
1 11 2
ˆ( , ) ( , ) ( , ) ( , )
( , ) { ( , )} ( , )M N
F
R u v G u v H u v F u v
r x y T R u v r x y− −
−
= −
⇒ = ⇒ = ∑∑2r0 0x y= =
• determinación de || n ||2
es posible implementar un algoritmo de restauración óptimo conociendo sola-mente la media y la varianza del ruido.
1 1 1 12 2
0 0 0 0
1 1( , ) ; [ ( , ) ]M N M N
x y x y
m x y x y mMN MNη η ηη σ η
− − − −
= = = =
= = −∑∑ ∑∑
1 12 2 2
0 0
1 1 1 1 1 12 2
0 0 0 0 0 0
[ ( , ) 2 ( , ) ]
( , ) 2 ( , ) 1
M N
x y
M N M N M N
x y x y x y
MN x y m x y m
x y m x y m
η η η
η η
σ η η
η η
− −
= =
− − − − − −
= = = = = =
⇒ = − +
= − +
∑∑
∑∑ ∑∑ ∑∑
1 1M N− −
1 12 2 2
0 0
1 12 2
0 0
( , ) 2
( , )
y y y
M N
x y
M N
x y
x y MNm MNm
x y MNm
η η
η
η
η
− −
= =
− −
= =
= − +
= −
∑∑
∑∑1 1
2 2 2
0 0( , ) [ ]
M N
x yx y MN mη ηη σ
= =
⇒ = = +∑∑2n
Curso: Procesamiento Digital de Imágenes Otoño 2019
Dr. G. Urcid - INAOE 17
Filtrado restringido
5 6
2 5
10 , 10 , 0.25
0 ; 10i f aγ γ
μ σ
− −
−
= = =
= =
5 6
2 2
10 , 10 , 0.25
0 ; 10i f aγ γ
μ σ
− −
−
= = =
= =
imagen original imagen degradada imagen restaurada
2 2 5/ 6( , ) exp[ ( ) ]H u v k u v= − +
5 60k = 0.0025k = 5 6
2 5
10 , 10 , 0.25
0 ; 10i f aγ γ
μ σ
− −
−
= = =
= =filtrado CLS iterativo
Curso: Procesamiento Digital de Imágenes Otoño 2019
Dr. G. Urcid - INAOE 18
Filtro GM-WienerFiltro combinado de la media geométrica (GM) y de Wiener,
dá una familia de filtros de restauración paramétricos
1 αα −⎡ ⎤⎡ ⎤ *
21 ( , )ˆ ( , ) ( , )
( , ) ( , ) ( , ) / ( , )f
H u vF u v G u vH u v H u v S u v S u v
α
ηβ
⎡ ⎤⎡ ⎤= ⎢ ⎥⎢ ⎥
+⎢ ⎥⎣ ⎦ ⎣ ⎦
• α = 1 , filtro inverso (IF)
• α = 0 , filtro paramétrico de Wiener (PWF)
• α = 0, β = 1, filtro estándar de Wiener (SWF)
• α = ½, media geométrica de IF y PWF
• α = ½, β = 1, filtro de igualación espectral (SEF)