Usando o filtro inverso para corrigir uma imagem espacialmente envolvida (deconvolução)

8

Como parte de uma tarefa de casa, estamos implementando o filtro inverso. Degrade uma imagem e depois recupere com um filtro inverso.

Convolo a imagem no domínio espacial com um filtro de caixa 5x5. FFT o filtro, FFT a imagem degradada e, em seguida, divido a imagem degradada pelo filtro. FFT inverso o resultado em uma imagem e recebo lixo.

Se eu FFT a imagem, FFT o filtro, multiplique os dois, divida o resultado pelo filtro da FFT, obviamente chego muito perto da imagem original. ((X * Y) / Y ~ == X)

Eu tenho uma idéia de que a matemática não é tão simples quanto "espacialmente convoluida == multiplicação FFT".

Qual é a maneira correta de usar o filtro inverso? Eu tenho o kernel exato usado degradar a imagem. Não estou adicionando nenhum ruído.

O livro de Bovik, The Essential Guide to Image Processing, é quase completamente desprezível do filtro inverso. Gonzalez & Woods é um pouco mais esperançoso, mas pula quase imediatamente para o Filtro Wiener.

Tenho uma pergunta semelhante no stackoverflow.com /programming/7930803/inverse-filter-of-spatially-convolved-versus-frequency-convolved-image

(Essas perguntas também devem ser marcadas como [lição de casa], mas a tag ainda não existe e eu não tenho o representante para criá-la.)

EDITAR. Para algumas das ótimas sugestões abaixo. @ dipan-mehta Antes de FFT, estou preenchendo o kernel de convolução com o mesmo tamanho da imagem. Estou colocando o kernel no canto superior esquerdo. Eu ifft (ifftshift ()), em seguida, salvo em uma imagem e recebo um bom resultado. Eu fiz o ifft (ifftshift ()) no kernel e na imagem. Bons (ish) resultados. (As imagens estão na minha /programming/7930803/inverse-filter-of-spatially-convolved-versus-frequency-convolved-image question.)

@ jason-r provavelmente está correto. Não entendo a matemática da convolução subjacente + transformar. "Deconvolução" era uma nova palavra para mim. Ainda tem muito a aprender. Obrigado pela ajuda!

Minha solução para a tarefa de casa é fazer tudo no domínio da frequência. Eu falei com o professor. Eu estava fazendo a tarefa mais difícil do que o necessário. Ela queria que adicionássemos ruído e experimentássemos o filtro inverso, o filtro Wiener e o filtro de mínimos quadrados restritos. O objetivo do exercício era ver como os filtros lidavam com o ruído.

David Poole
fonte
11
Você preenche o filtro com zeros para torná-lo do mesmo tamanho da imagem antes de tirar a FFT? Você está fazendo a divisão complexa corretamente?
Dima
Sim, preenchido o filtro com zeros, kernel no canto superior esquerdo. Todo o meu código Python / numpy está no link stackoverflow.com acima mencionado. A divisão complexa é provavelmente o meu problema.
David Poole

Respostas:

8

Existem algumas subquestões que abordarei separadamente:

  • Convolução no domínio espacial (ou correspondentemente no domínio do tempo para sinais amostrados no tempo) é equivalente à multiplicação no domínio da frequência. Nos sistemas amostrados, existem algumas sutilezas nos casos de fronteira (ou seja, ao usar o DFT, a multiplicação no domínio da frequência realmente oferece convolução circular, não convolução linear), mas, em geral, é realmente simples assim.

  • A filtragem inversa pura quase nunca é a solução certa na prática. Na maioria dos casos, você não tem acesso ao filtro exato que foi aplicado aos seus dados; portanto, não pode simplesmente invertê-lo. Mesmo se você conhece o filtro, ainda é problemático. Considere o fato de que o filtro pode ter zeros em determinadas frequências espaciais; se isso acontecer, depois de aplicar o filtro à sua imagem, todas as informações nessas frequências serão perdidas. Se você inverter ingenuamente esse filtro, ele terá um ganho infinito (ou pelo menos muito alto) nesses valores nulos. Se você aplicar o inverso ingênuo a uma imagem que tenha qualquer conteúdo aditivo nessas frequências (por exemplo, ruído, o que provavelmente será o caso), esse componente não assistido será bastante amplificado. Isso geralmente não é desejável.

    Esse problema de filtragem inversa é muito semelhante à equalização em sistemas de comunicação, onde esse fenômeno é chamado de aprimoramento de ruído . Nesse contexto, a abordagem de filtro inverso é referida como um equalizador de força zero , que raramente é realmente usado.

  • A área que você está explorando é conhecida em geral como deconvolução . Como regra geral, a deconvolução é uma operação complicada. Mesmo se você souber o filtro exato que foi aplicado e quiser desfazê-lo, nem sempre é fácil. Como você observou, a abordagem de filtro inverso geralmente é descartada em favor de um filtro Wiener ou de alguma outra estrutura que não visa exatamente inverter o sistema, mas sim estimar qual foi a entrada do sistema e minimizar algum critério de erro (minimizar a média erro quadrático é um objetivo comum). Como seria de esperar, aplicar um filtro Wiener a esse problema é conhecido como deconvolução Wiener .

Jason R
fonte
"Deconvolução" era uma nova palavra para mim. Eu ainda tenho muito que aprender. Obrigado!
David Poole
@JasonR não precisaríamos conhecer uma sequência de 'instrutores' na imagem para a deconvolução Wiener, de modo que o critério MMSE seja minimizado com relação a algo que se sabe ser verdade?
Spacey
11
Em geral, você precisa conhecer a densidade espectral de potência do sinal e a função de transferência que foi aplicada ao sinal para projetar o filtro Wiener. No entanto, no caso provável de você não conhecer todas essas informações, é possível fazer suposições educadas que produzem uma estrutura funcional mais robusta que o filtro inverso. Veja esta seção na página da Wikipedia para uma discussão.
Jason R
3

Espero que você não tenha cometido um erro na maneira como a computação é feita -

Convolo a imagem no domínio espacial com um filtro de caixa 5x5. FFT o filtro, FFT a imagem degradada e, em seguida, divido a imagem degradada pelo filtro. FFT inverso o resultado em uma imagem e recebo lixo.

Suponha que sua imagem tenha tamanho de 256x256 e filtro de 5x5 - para aplicar a filtragem multiplicando as FFTs, você deve primeiro converter o filtro em tamanho equivalente primeiro. Para isso, você deve manter o filtro da caixa 5x5 no "canto SUPERIOR" (não no centro da imagem) e preenchê- lo com zeros para preencher 256x256 - você deve obter uma FFT de 256x256 para o filtro.

Para ajudar a diagnosticar, na etapa 1 da programação - primeiro pegue apenas 256 x 256 FFT de filtro sozinho e verifique se a rotina IFFT - é capaz de devolver o filtro. Da mesma forma, teste se FFT -> IFFT da própria imagem funciona corretamente para trás.

Etapa 2 - se você aplicar apenas o filtro (e nenhum filtro inverso) no domínio FFT multiplicando - verifique a imagem resultante depois que o IFFT estiver bom. Deve ser basicamente imagem borrada.

Se toda a sua programação estiver correta - assegure-se de que quando você faz 1 / x para coeficiente de FFT, não haja erro de divisão por zero e, inversamente, quando houver muito pico nas multiplicações, haverá distorções pesadas.

Em geral - para qualquer filtro estável , o filtro inverso por definição instável - esse pode ser o principal motivo. No entanto, eu sempre gostaria de verificar a implementação antes de explorar os limites teóricos.

Se bem feito, tenho visto a multiplicação em FFT é convolução no espaço da amostra, tanto para imagens quanto para sinais de áudio.

Dipan.

Dipan Mehta
fonte