Análise numérica de algoritmos de tratamento de imagens científicas - perigos e potencialidades

Análise numérica de Algoritmos de Tratamento de Imagens Científicas – Perigos e Potencialidades Por Inês Ferreira Figueiredo Tese apresentada à Faculdade de Engenharia da Universidade do Porto para obtenção do grau de Mestre Porto, Julho 2012 Análise numérica de Algoritmos de Tratamento de Imagens Científicas – Perigos e Potencialidades Por Inês Ferreira Figueiredo Tese apresentada à Faculdade de Engenharia da Universidade do Porto para obtenção do grau de Mestre Orientadoras: Professora Filomena Dias d’Almeida Professora Maria João Rodrigues Porto, Julho 2012 i "The future will depend on what we do in the present." Mahatma Gandhi ii Agradecimentos A realização deste trabalho só foi possível graças à participação direta e indireta de um grande conjunto de pessoas. Como tal, quero deixar uma palavra de profundo agradecimento a todos eles. Em particular, agradeço às minhas orientadoras, pela força e confiança que dedicaram à minha pessoa e ao meu trabalho e ainda, ao Doutor Duarte Pignatelli, que atenciosamente nos forneceu as imagens médicas sobre as quais trabalhei. Um agradecimento muito especial, aos meus pais, que me permitiram chegar até aqui, contribuindo não só economicamente mas, acima tudo, afetivamente, incentivando o meu crescimento pessoal. Aos meus amigos e de mais família, um muito obrigada pelo carinho e amizade. E a todos aqueles que foram passando pela minha vida e de uma forma ou de outra foram tendo o seu papel nela, ajudando a tornar-me naquilo que sou hoje, um agradecimento sincero. iii Resumo Análise numérica de algoritmos de tratamento de imagens cientificas – perigos e potencialidades A reconstituição de imagens desfocadas é um problema de interesse para diversas áreas como, por exemplo, a Medicina, pelo que surge a necessidade de se estudarem modelos de desfocagem artificial de imagens. Em processos de Reconstituição de Imagem, a Decomposição em Valores Singulares, SVD, é uma ferramenta que pode ser útil como ferramenta alternativa para resolução do modelo linear usado no tratamento de imagens desfocadas. O uso da SVD truncada é importante no sentido de permitir reduzir a dimensão dos problemas e permitir também reconstituir a imagem com boa qualidade e em menos tempo quando as dimensões são grandes. Palavras-chave: Desfocagem de Imagem, SVD, Reconstituição de Imagem, SVD truncada. iv Abstract Numerical analysis of algorithms for Treatment of Scientific Images – advantages and misunderstandings Reconstitution of blurred images is an issue of concern for several areas such as the Medicine, so that there is a need to study models of artificial image blur. In proceedings of Image Reconstruction, the Singular Value Decomposition, SVD can be a useful alternative tool for solving the linear model used in the treatment of blurred images. The use of SVD truncated is important in order to reduce the dimension of problems and also allow reconstruct an image with good quality and in less time when the dimensions are large. Keywords: Deblurring Process, SVD, Image Reconstitution, SVD truncated v Conteúdo Agradecimentos .............................................................................................................. ii Resumo ........................................................................................................................... iii Abstract .......................................................................................................................... iv Lista de Figuras ............................................................................................................ vii Lista de Tabelas ............................................................................................................. ix Lista de Acrónimos ......................................................................................................... x Capítulo 1. Introdução ............................................................................................. 1 1.1. Objetivo e Motivação ......................................................................................... 1 1.2. Estrutura do Relatório ........................................................................................ 3 Capítulo 2. Processo de Desfocagem e Reconstituição de uma Imagem .............. 4 2.1. Descrição do Problema ...................................................................................... 4 2.2. Modelo Linear de Desfocagem de uma Imagem ............................................... 5 2.3. O recurso à SVD para refocar ............................................................................ 7 2.4. Exemplo ............................................................................................................. 8 Capítulo 3. Aplicação a Imagens Científicas ........................................................ 10 3.1. Imagem “Joaninha” .......................................................................................... 10 3.2. Imagens médicas de mitocôndrias ................................................................... 20 3.2.1. Imagem de mitocôndrias de boa qualidade .............................................. 21 Principais Resultados .............................................................................................. 22 Resultados recorrendo à SVD truncada.................................................................. 32 3.2.2. Imagem desfocada de mitocôndrias ......................................................... 33 Capítulo 4. Decomposição em Valores Singulares ............................................... 35 4.1. Interpretação geométrica da SVD .................................................................... 36 4.2. Definição de SVD ............................................................................................ 36 vi 4.3. Decomposição em Valores Singulares Truncada ............................................ 38 4.4. Cálculo da Decomposição em Valores Singulares .......................................... 39 4.4.1. Cálculo da SVD a partir dos valores próprios .......................................... 40 4.4.2. Cálculo da SVD sem calcular A*A e usando a Bidiagonalização de Golub Kahan 41 4.4.3. Cálculo da SVD sem calcular A*A e usando a Bidiagonalização de Lawson-Hanson-Chan ............................................................................................. 45 4.4.4. Cálculo da SVD sem calcular A*A e usando a Bidiagonalização em 3 passos 46 4.4.5. Comparação entre os diferentes métodos de Bidiagonalização................ 46 4.5. Aplicações da Decomposição em Valores Singulares ..................................... 48 4.5.1. Aplicações da SVD no tratamento de Imagem ......................................... 48 Capítulo 5. Conclusão e Trabalho Futuro ............................................................ 52 Referências Bibliográficas ........................................................................................... 55 Apêndice A – Programas e Scripts .............................................................................. 56 Apêndice B – Tabelas de Resultados .......................................................................... 61 B-1. Alguns dos resultados obtidos, aplicando Processo de Desfocagem e Recuperação da Imagem, para a imagem “Joaninha” ................................................. 62 B-2. Alguns dos resultados obtidos, aplicando Processo de Desfocagem e Recuperação da Imagem, para a imagem de mitocôndrias de melhor qualidade ....... 66 B-3. Resultados obtidos para a Imagem de mitocôndrias desfocada ....................... 72 vii Lista de Figuras Figura 1 Imagem original antes de submetida a um processo de Desfocagem ................ 8 Figura 2 Imagem desfocada e posteriormente recuperada utilizando um modelo linear . 8 Figura 3 Imagem convertida para níveis de cinzento, X (à esq.); Imagem desfocada, B, (PSF propriamente dita, raio 3 e dimensão da imagem) (ao centro); Imagem recuperada, Xreconstruído (à dir.) .............................................................................................................. 9 Figura 4 Imagem “Joaninha”, [500,500] ........................................................................ 11 Figura 5 Gráfico da norma da diferença entre a imagem original e a recuperada, considerando ou não o valor do erro e usando ou não a SVD, para uma PSF de raio 2 13 Figura 6 Imagem original (à esq.), Imagem desfocada, B, usando PSF de dimensão [500x500], raio 2 e condições fronteira “periodic” (ao centro), Imagem recuperada (à dir.) ................................................................................................................................. 13 Figura 7 Imagem original (à esq); Imagem desfocada ................................................... 13 Figura 8 Imagem recuperada a partir do uso da SVD truncada, k =400, raio da PSF igual a 4 ................................................................................................................................... 20 Figura 9 Imagem original das mitocôndrias, [350, 350] ................................................ 21 Figura 10 Imagem original (à esq.) e Imagem recuperada após desfocagem usando uma PSF de raio 2 (à dir.) ....................................................................................................... 24 Figura 11 Imagem original (à esq.), Imagem recuperada após desfocagem com PSF de dimensão [350,350], Raio = 3 e condições fronteira “periodic” ................................... 24 Figura 12 Gráfico da norma da diferença entre a imagem original e a imagem recuperada, introduzindo ou não o erro e usando ou não o erro, a partir de uma PSF de raio = 3 e condições fronteira “periodic” ....................................................................... 26 Figura 13 Imagem original (à esq. ); Imagem recuperada sem a subtracção do parâmetro ruído invertido (valor do erro 1x10-6), para uma PSF de raio 3 e dimensão [350, 350]. 26 4 2.1. Descrição do Problema A necessidade de recordar um momento leva-nos a recorrer, continuamente, ao uso de máquinas fotográficas, telemóveis e outras tecnologias que nos permitam adquirir imagens. No entanto, nem sempre a imagem que obtemos é aquela que desejaríamos. Na prática, uma imagem pode sofrer uma desfocagem por diversos motivos. Um desajuste no alinhamento das lentes, uma alteração atmosférica, erros por parte do utilizador, através de movimentos involuntários, são alguns deles. Em microscopia confocal, por exemplo, pode ocorrer uma distorção no alinhamento da luz incidente, que leva à obtenção de uma imagem esbatida. Capítulo 2. Processo de Desfocagem e Reconstituição de uma Imagem 5 O passo inicial deste trabalho consistiu em tentar simular alguma destas situações. Como tal, submetemos uma fotografia a um processo artificial de desfocagem e posteriormente tentamos recuperá-la. É importante realçar, que num processo artificial de desfocagem, torna-se mais fácil recuperar a imagem, uma vez que nós conhecemos a causa da desfocagem e temos a imagem original para medir a qualidade da reconstituição, o que infelizmente na prática não acontece. Neste capítulo, pretendemos fazer uma abordagem aos processos de reconstituição de imagens. 2.2. Modelo Linear de Desfocagem de uma Imagem De uma forma geral, num processo de desfocagem representa-se a imagem através de um modelo linear. A imagem é uma matriz X. Como sofreu desfocagens ou ficou afetada de erro com a digitalização, em vez de X temos efetivamente B, imagem desfocada. Sejam as matrizes Ac e Ar a matriz de alteração das colunas da imagem e das linhas da imagem, respetivamente. Então a matriz desfocada é dada por [1]: (1) Alternativamente, este modelo permite escrever (2) Dado que a imagem armazenada também é afetada por ruído, tanto mais não seja devido à precisão máquina finita com que os pixels são representados na versão digital, é necessário acrescentar uma matriz ruído, E, ao modelo linear (1): (3) 6 O que implica que (4) isto é, acrescenta-se um novo parâmetro ao nosso modelo a que chamamos ruído invertido. Uma vez que a desfocagem de uma imagem pode ter diferentes causas é importante percebermos que cada imagem desfocada tem a sua própria história. Assim, é fácil de entender que a desfocagem é um processo complexo que requer um estudo cuidadoso e minucioso. Atualmente existe um conjunto de modelos matemáticos que nos permitem representar vários tipos de desfocagem de uma imagem. Consoante a situação que queremos simular, existe um conjunto de alternativas às quais podemos apelar. Segundo [1] podemos considerar uma desfocagem propriamente dita, ou seja, atribuímos à imagem uma distorção radial centrada num único ponto da imagem ou em vários pontos, uma desfocagem que simula a turbulência atmosférica ou a uma desfocagem horizontal (ou vertical) ao longo de uma linha. Todos estes métodos têm como base a criação de uma matriz PSF, “Point Spread Function”, que não é mais do que uma função na qual está descrita a informação relativa ao processo de desfocagem. A PSF depois de criada é aplicada à imagem de forma a obtermos uma nova imagem, agora desfocada. A criação da PSF é o que nos permite extrair as matrizes Ac e Ar, que caraterizam a forma da desfocagem, através de uma Decomposição de Kronecker [1, pág. 48]. Um outro aspeto a considerar é o tipo de condições fronteira que atribuímos à alteração da imagem. Habitualmente consideram-se três tipos de condições fronteira, isto é, três tipos de parâmetros que nos permitem descrever o formato da periferia da desfocagem [1]. Em imagens provenientes de áreas como a astronomia, a condição fronteira mais comummente utilizada consiste em considerar que a periferia da imagem é toda ela sem informação (a 7 preto). Outro tipo de condição fronteira consiste em atribuir aos bordos da imagem um carácter periódico em que a informação da imagem se reproduz em diferentes direções. E por último pode definir-se o exterior da imagem como sendo um espelho do seu interior, ou seja, a informação da imagem está replicada no seu exterior. Depois de definirmos as condições fronteira também procedemos à determinação da dimensão da matriz PSF que traduz o raio de alcance da desfocagem. A influência que a dimensão e o raio da PSF têm na desfocagem e posterior reconstituição da imagem será abordada mais à frente neste trabalho. Nos casos em que as componentes vertical e horizontal de desfocagem podem ser separadas, a matriz PSF pode ser representada como um produto de Kronecker de duas matrizes Ac e Ar (já referidas) que dependem apenas do tipo de PSF escolhido e respectivos parâmetros (centro e raio) e do tipo de condições fronteira [1]. Outro dos aspetos que consideramos neste trabalho consistiu em sujeitar a imagem não só a uma desfocagem pela PSF, mas também, a diferentes valores de ruído e testar a capacidade do modelo em recuperar a imagem com a intervenção deste parâmetro. 2.3. O recurso à SVD para refocar Depois de desfocarmos a imagem, uma alternativa para a recuperação da imagem passa pelo uso da SVD. No quarto capítulo, faremos uma análise mais aprofundada sobre esta técnica, no entanto, é de referir que esta técnica tem importantes aplicações não só na recuperação da imagem, mas também na otimização da duração do processo. Estudaremos ainda, o recurso à SVD truncada, que permite recuperar a imagem a partir de uma só parte da informação sobre a imagem. O uso da SVD é mais um aspeto interessante de analisar, uma vez que se levantam questões sobre a sua eficácia. Maior rapidez, igual eficiência? O uso de menos informação, resulta em igualmente bons resultados? Estes serão algumas das dúvidas que discutiremos mais à frente. 8 2.4. Exemplo Consideremos a seguinte imagem “Joaninha”: Em Tratamento de Imagem o processo torna-se mais simples se trabalharmos com imagens a níveis de cinzento. Depois disso, escolhemos um processo de desfocagem a aplicar à imagem. Neste exemplo, consideramos uma desfocagem propriamente dita. Um aspeto importante que consideramos neste trabalho foi perceber que parâmetros podem aumentar ou diminuir a semelhança entre a imagem original e a recuperada. Vejamos o que foi acontecendo à imagem “Joaninha”. Depois da conversão para níveis de cinzento, como já referimos, criamos a PSF com uma dimensão igual à da imagem e um raio de desfocagem de, por exemplo, 3 unidades. Figura 1 Imagem original antes de submetida a um processo de Desfocagem Figura 2 Imagem desfocada e posteriormente recuperada utilizando um modelo linear 9 Atribuindo à imagem Condições Fronteira do tipo ‘zero’, fomos obtendo os seguintes resultados (figura 3): Visualmente, a recuperação da imagem parece excelente, e também o podemos comprovar numericamente. Para isso, recorrermos à norma euclidiana da diferença entre as matrizes das duas imagens, ‖ ‖. No exemplo apresentado, o que fizemos foi gerar uma matriz aleatória, a qual multiplicamos por diferentes valores e somámo-la à imagem desfocada, B. Pelo que, passamos a ter uma nova matriz que, para além da desfocagem, sofreu uma alteração mínima provocada por causa desconhecida. Mais à frente, veremos a forma como a subtração, ou não, deste parâmetro afeta a recuperação da imagem. Quando tratamos imagens de grandes dimensões, os processos acabam sempre por se tornar ligeiramente morosos, pelo que existe uma necessidade de encontrar soluções que minimizem o custo do processo. Figura 3 Imagem convertida para níveis de cinzento, X (à esq.); Imagem desfocada, B, (PSF propriamente dita, raio 3 e dimensão da imagem) (ao centro); Imagem recuperada, Xreconstruído (à dir.) 10 A melhor forma de compreender uma metodologia é trabalhar sobre ela. Assim, neste capítulo apresentaremos alguns dos testes que fizemos sobre diferentes imagens científicas, cuja qualidade, no que diz respeito à nitidez e ao contraste, varia. 3.1. Imagem “Joaninha” Numa primeira abordagem ao problema da desfocagem e reconstituição de uma imagem recorremos à imagem “Joaninha”. De seguida, abordaremos o processo de forma mais aprofundada, apresentando alguns dos resultados dos testes que realizamos sobre esta mesma imagem. Pelo facto de a imagem apresentar um elevado peso computacional, fizemos um reajuste à sua dimensão, sendo que trabalhamos com uma imagem de 500 por 500 pixels (figura 4). Capítulo 3. Aplicação a Imagens Científicas 11 Figura 4 Imagem “Joaninha”, [500,500] Como já referimos, o modelo de desfocagem e reconstituição de imagens é complexo, e como tal, surgiu a necessidade de estudar a forma como os diferentes parâmetros influenciam o processo. Assim, começamos por testar o modelo, variando a dimensão da PSF, o valor do raio, as condições fronteira e o valor escalar do erro contido no ruído invertido. Começamos o processo por fixar o tipo de condição fronteira, variando alternadamente os outros parâmetros. Dimensão da PSF constante Numa primeira fase, desfocamos a imagem a partir de uma PSF com dimensão igual à da nossa imagem, 500 por 500 pixels, e condições fronteira do tipo “periodic”. Os parâmetros alternáveis foram o raio e o ruído invertido. A recuperação da imagem consistiu em dois focos fundamentais: o uso do modelo linear e o mesmo modelo reajustado ao uso da SVD. Para cada um destes modelos, o que fizemos foi considerar a ausência ou presença do parâmetro ruído invertido e determinar a sua menor ou maior influência nos resultados, de acordo com o seu valor. Mais à frente, apresentaremos alguns dos resultados obtidos com posterior análise dos mesmos. 12 Depois de realizarmos vários testes alternando apenas o raio e o valor do ruído invertido, o que fizemos foi considerar diferentes dimensões para a PSF, tendo sempre o cuidado de lhe atribuir uma dimensão duas vezes superior ao valor do raio no momento considerado. E por fim, criamos um conjunto de testes para estes mesmos parâmetros, porém atribuindo à nossa PSF condições fronteira do tipo “zero”. Ao longo de todos estes testes fomos considerando os tempos de duração de cada um dos processos, de forma a perceber se existem diferenças significativas entre os mesmos. Um aspeto importante a considerar em todo o processo é o condicionamento de cada uma das matrizes, Ac e Ar, uma vez que é relevante perceber a influência que cada um destes parâmetros tem sobre esta característica da matriz. A segunda fase de testes sobre esta imagem consistiu no uso da SVD truncada para reconstruir a imagem. Também para este processo recorremos a diferentes tipos de condições fronteira, valores de raio, ruído invertido e diferentes dimensões da PSF. A grande diferença relativamente a esta metodologia está no facto de nesta situação fazermos uma aproximação da matriz, recorrendo apenas a informação até k valores singulares (os maiores em módulo). O conjunto de testes que executámos sobre as imagens foi vasto, pelo que surgiu a necessidade de criarmos algumas folhas de cálculo Excel para armazenarmos os resultados. Em seguida, apresentaremos algumas tabelas de resultados com as respetivas análises. O restante dos resultados poderá ser consultado nos anexos deste trabalho. Na realização destes testes, o primeiro aspeto a observar foi a ineficácia do processo quando usamos os parâmetros condição fronteira do tipo “periodic” e raio igual a 2, uma vez que obtivemos diferenças entre as imagens, original e reconstruida (figura 5), da ordem dos 1015 e, portanto, não conseguimos reconstituir a imagem (figura 6). 13 À medida que fomos progredindo no valor do raio obtivemos melhores resultados. Após vários testes, verificamos que os melhores resultados foram obtidos utilizando raio igual a 3 (tabela 1) e 4, uma vez que as diferenças, em norma, entre as imagens foram da ordem dos 10-10. A partir deste valor as diferenças começam a aumentar ligeiramente. Figura 7 Imagem original (à esq); Imagem desfocada 1,00E+04 1,00E+06 1,00E+08 1,00E+10 1,00E+12 1,00E+14 1,00E+16 1,00E+18 1,00E+20 1,00E+22 1,00E+24 1,00E+26 1,00E+28 1,00E+30 Gráfico da norma da Diferença entre Imagem original e recuperada Xreconstruido_com_ruído Xreconstruido_sem_ruído Xreconstruido_svd_com_ruído Xreconstruido_svd_sem_ruído Figura 5 Gráfico da norma da diferença entre a imagem original e a recuperada, considerando ou não o valor do erro e usando ou não a SVD, para uma PSF de raio 2 Figura 6 Imagem original (à esq.), Imagem desfocada, B, usando PSF de dimensão [500x500], raio 2 e condições fronteira “periodic” (ao centro), Imagem recuperada (à dir.) 20 Embora, à partida, a imagem recuperada possa parecer de boa qualidade, é necessário analisá-la cuidadosamente. Se observarmos a imagem num plano maior, podemos reparar em certos pormenores e perceber, portanto, que a recuperação não foi assim tão boa, nomeadamente nas zonas da periferia em que se nota um ligeiro esbatimento (figura 8). Para outros valores de raio e dimensões da PSF, os resultados foram similares, embora utilizando condições fronteira do tipo “periodic” se tenham obtido resultados ligeiramente inferiores. 3.2. Imagens médicas de mitocôndrias As mitocôndrias são a principal fonte de energia das células. Estes organelos, produzem, a partir de um processo aeróbio, uma molécula chamada ATP (Trifosfato de Adenosina). O predomínio deste organelo no nosso organismo dá-se nas células musculares, cerebrais e células responsáveis pela produção de hormonas, uma vez que estas são as que necessitam de energia permanentemente. Figura 8 Imagem recuperada a partir do uso da SVD truncada, k =400, raio da PSF igual a 4 21 A análise da estrutura destes organelos é muito importante no estudo de doenças e, portanto, é necessário ter ao nosso dispor fotografias, da mais elevada qualidade, para interpretar da melhor forma possível os problemas que afetam a mitocôndria. Nem sempre, contudo, a imagem possui a qualidade que desejaríamos e essa situação impede a correta interpretação do problema. A solução nestes casos é repetir o procedimento porém, uma vez que este processo é invasivo, essa não será a melhor solução. Assim, surge a utilidade de estudar alternativas que permitam corrigir a nitidez e o contraste da imagem. No que segue, aplicamos os algoritmos do capitulo 2 a imagens de mitocôndrias obtidas com o microscópio eletrónico, armazenadas como fotografias analógicas e depois digitalizadas para o nosso estudo. No primeiro caso, usamos uma imagem boa para servir de referência, desfocamos e recuperamos como em 3.1. No segundo caso, usamos uma imagem desfocada e tentamos recuperar com a SVD truncada. 3.2.1. Imagem de mitocôndrias de boa qualidade A primeira imagem da mitocôndria (figura 9), sobre a qual trabalhamos foi, juntamente com a outra que estudaremos, obtida a partir de um microscópio eletrónico. Figura 9 Imagem original das mitocôndrias, [350, 350] 22 Sobre esta primeira imagem aplicamos diversos testes variando os parâmetros que cada metodologia engloba. Consideramos uma desfocagem propriamente dita sobre a qual aplicamos diferentes valores de raio e condições fronteira. Também adicionamos à imagem algum ruído aleatório de diferentes grandezas. Principais Resultados Nesta primeira tabela, apresentamos alguns dos resultados que obtivemos para esta imagem quando lhe aplicamos um processo de desfocagem. Principais características:  Imagem redimensionada para uma matriz quadrada de 350 por 350 pixels;  Desfocagem propriamente dita;  Dimensão da PSF igual à da imagem original;  Raio variado entre 2 e 180;  Ruído Invertido, introduzido a partir de uma matriz aleatória, de dimensão igual à da imagem, com valores de erro escalar que variam entre 1x10-1 e 1x10-10;  Condições fronteira do tipo ‘periodic’. Uma vez que o valor do ruído varia dez vezes para cada valor de raio, é fácil percebermos a grande quantidade de resultados e imagens que obteremos. Comecemos pelo uso do Modelo Linear sem recurso à SVD. Para um Raio de dimensão 2, temos um Condicionamento de Ac igual a 2,28x1016 e um Condicionamento de Ar igual a 8,16x1015, valores bastantes elevados. Sem introduzirmos o parâmetro erro invertido, a nossa imagem não recupera o aspecto original, apresentando um tempo de execução de aproximadamente 0,0624 segundos e uma diferença em relação à original da ordem dos 9,60x1015. Da mesma forma, quando introduzimos o ruído invertido à imagem, os resultados também não são satisfatórios, quer posteriormente subtraiamos esse valor à imagem ou não. Porque o condicionamento das matrizes Ac e Ar é demasiado elevado (Tabela 5). 23 Tabela 6 Resultados obtidos a partir de uma desfocagem propriamente dita, com condições fronteira “periodic” e raio 2 Ruído invertido (valor erro) Recupera a imagem (Sem subtrair o ruído)? Diferença com a original Tempo (seg.) Recupera a imagem (Subtraindo o ruído)? Diferença com a original Tempo (seg.) 1,00E-01 Não 5,67E+30 0.0624 Não 5,35E+16 0.0624 1,00E-02 Não 4,88E+29 0.0312 Não 1,28E+16 0.0312 1,00E-03 Não 6,98E+28 0.0624 Não 7,00E+15 0.0312 1,00E-04 Não 4,17E+27 0.0312 Não 8,06E+15 0.0624 1,00E-05 Não 4,66E+26 0.0468 Não 1,30E+16 0.0312 1,00E-06 Não 3,98E+25 0.0624 Não 1,43E+16 0.0468 1,00E-07 Não 5,02E+24 0.0624 Não 1,86E+16 0.0624 1,00E-08 Não 5,42E+23 0.0624 Não 8,74E+15 0.0624 1,00E-09 Não 7,53E+22 0.0312 Não 1,60E+16 0.0312 1,00E-10 Não 7,61E+21 0.0624 Não 1,58E+16 0.0312 Os resultados, quando aplicamos o processo de refocagem recorrendo à SVD, como seria de esperar, não diferem muito dos anteriores, sendo que a recuperação da imagem sem a introdução do parâmetro ruído invertido foi de 1,47E+16, com um tempo de 0.0624 segundos (ver apêndice B). Se analisarmos atentamente os resultados, podemos ver que à medida que diminuímos o valor do ruído invertido, a qualidade da imagem melhora ligeiramente, no entanto, os valores continuam muito discrepantes, o que não nos permite recuperar a imagem. Como veremos mais à frente, sempre que aplicamos à imagem uma desfocagem com um raio de 2, os resultados são maus. A imagem abaixo mostra um dos resultados obtidos quando recuperamos a imagem recorrendo a uma PSF da dimensão da imagem original com raio 2. 24 Figura 10 Imagem original (à esq.) e Imagem recuperada após desfocagem usando uma PSF de raio 2 (à dir.) Quando realizamos o mesmo teste sobre a imagem recorrendo agora a um raio de dimensão 3, os resultados foram diferentes, melhorando significativamente (figura 11), uma vez que o condicionamento das matrizes diminuiu bastante, para aproximadamente 872 (Tabela 7). Figura 11 Imagem original (à esq.), Imagem recuperada após desfocagem com PSF de dimensão [350,350], Raio = 3 e condições fronteira “periodic” 25 Tabela 7 Resultados obtidos na recuperação da imagem, após desfocagem com uma PSF de raio 3, dimensão [350, 350] e condições fronteira “periodic” Ruído invertido (valor erro) Recupera a imagem (Sem subtrair o ruído)? Diferença com a original Tempo (seg.) Recupera a imagem (Subtraindo o ruído)? Diferença com a original Tempo (seg.) 1,00E-01 Não 4.2990E+004 0.0624 Sim 1.4285E-010 0.0624 1,00E-02 Não 4.1185E+003 0.0312 Sim 1.2966E-010 0.0468 1,00E-03 Não 442.0884 0.0624 Sim 9.9397E-011 0.0312 1,00E-04 Não 47.6699 0.0156 Sim 1.0643E-010 0.0624 1,00E-05 Não 4.6360 0.0312 Sim 7.8978E-011 0.0624 1,00E-06 Sim 0.6008 0.0624 Sim 1.1651E-010 0.0312 1,00E-07 Sim 0.0419 0.0624 Sim 1.0354E-010 0.0624 1,00E-08 Sim 0.0036 0.0468 Sim 1.1807E-010 0.0624 1,00E-09 Sim 4.3222E-004 0.0624 Sim 7.1269E-011 0.0624 1,00E-10 Sim 3.6431E-005 0.0312 Sim 1.4653E-010 0.0468 Com este valor do raio sem considerar o ruído obtivemos uma diferença entre as imagens de 9,19E-11 para o modelo linear e de 6,09E-11 para o modelo linear recorrendo à SVD, ambos com um tempo de execução de 0,0624 segundos. Sempre que somamos o ruído invertido, os resultados foram ligeiramente piores. No gráfico seguinte, podemos analisar a influência que o ruído tem na recuperação da imagem, uma vez que, quando o introduzimos e não o retiramos posteriormente, os resultados sofrem graves perdas. 26 Figura 12 Gráfico da norma da diferença entre a imagem original e a imagem recuperada, introduzindo ou não o erro e usando ou não o erro, a partir de uma PSF de raio = 3 e condições fronteira “periodic” Podemos verificar que com um valor do erro escalar inferior a 1x10-6, a diferença entre as imagens não é visível, pelo que podemos assumir uma boa recuperação da imagem independentemente desse valor (figura 13) Figura 13 Imagem original (à esq. ); Imagem recuperada sem a subtracção do parâmetro ruído invertido (valor do erro 1x10-6), para uma PSF de raio 3 e dimensão [350, 350] O processo repetiu-se para valores de raio de 4, 5, 6, 7, 8, 10, 20, 50, 100, 175 e 180. O que se verificou ao longo destes testes foi uma ligeira perda de precisão à medida que aumentamos o valor do raio, que se pode medir com a norma da matriz diferença mas que 1,000E-11 1,000E-10 1,000E-09 1,000E-08 1,000E-07 1,000E-06 1,000E-05 1,000E-04 1,000E-03 1,000E-02 1,000E-01 1,000E+00 1,000E+01 1,000E+02 1,000E+03 1,000E+04 1,000E+05 Gráfico da norma da Diferença entre Imagem original e recuperada Xreconstituido Xreconstituido_sem_ruido Xreconstituido_svd Xereconstituido_svd_sem_ruido 27 não vemos à vista desarmada. Outro aspeto a realçar tem que ver com o facto de se verificarem melhores resultados quando definimos um raio ímpar (figura 14). A diferença de resultados quando utilizamos um raio igual a 175 e um raio igual a 180 é bastante significativa, uma vez que quando utilizamos o raio 180 já atingimos um valor superior a metade da dimensão da imagem e portanto, obtemos resultados cada vez menos satisfatórios até um ponto em que a imagem já não pode ser recuperada (figura 15). 1,00E-11 1,00E-10 1,00E-09 1,00E-08 1,00E-07 1,00E-06 1,00E-05 1,00E-04 1,00E-03 1,00E-02 1,00E-01 1,00E+00 Gráfico da norma da Diferença entre Imagem original e recuperada Xreconstituido_sem_ruido_PSF_R=3 Xereconstituido_svd_sem_ruido_PSF_R=3 Xreconstituido_sem_ruido_PSF_R=4 Xreconstituido_svd_sem_ruido_PSF_R=4 Figura 14 Resultados obtidos para a reconstituição da imagem usando uma PSF com diferentes valores de raio 28 Figura 15 Gráfico da norma da diferença entre a imagem original e a recuperada, após desfocagem com PSF de raio 175 e 180 PSF de dimensão variável De seguida, aplicamos à imagem uma PSF de dimensão variável, ou seja, atribuímos diferentes dimensões para a PSF e sobre cada uma delas aplicamos diferentes valores de raio. Nesta etapa, verificamos novamente um inviabilidade do processo para um raio igual a 2, independentemente da dimensão da PSF. Por exemplo, para uma dimensão da PSF de 5 por 5 pixels os resultados obtidos não permitiram recuperar a imagem, o que já seria de esperar pelo facto de o único valor possível de raio para esta dimensão de PSF ser igual a 2 e como já vimos este valor não funciona no processo. Por outro lado, quando definimos uma PSF de 10 por 10 pixels os resultados melhoraram substancialmente. Para um raio igual a 3, obtivemos resultados da ordem dos 10-11 e novamente, verificamos que a partir de um determinado valor de erro, nomeadamente 1x10-6, conseguimos recuperar a imagem ainda que sem a subtrair esse valor ao resultado da imagem. 1,000E-08 1,000E-06 1,000E-04 1,000E-02 1,000E+00 1,000E+02 1,000E+04 1,000E+06 1,000E+08 1,000E+10 1,000E+12 Gráfico da norma da Diferença entre Imagem original e recuperada Xreconstituido_com_ruido_PSF_R=175 Xereconstituido_sem_ruido_PSF_R=175 Xreconstituido_com_ruido_PSF_R=180 Xreconstituido_sem_ruido_PSF_R=180 29 Continuamos os testes para dimensões de PSF de 20 por 20, 50 por 50, 100 por 100 e 300 por 300 pixels, variando o raio para cada uma (tabela ). Tabela 8 Resultados obtidos na reconstituição da imagem após uma desfocagem com uma PSF com condições fronteira “periodic” Dimensão da PSF Valor do Raio Norma da diferença entre as imagens (com o ruido 1x10-10) Recupera a imagem? Norma da diferença entre as imagens (subtraindo o ruido) Recupera a imagem? [5,5] 2 3 2.1906e+021 1.7755e+022 Não Não 4.5337e+015 2.0260e+016 Não Não [10,10] 3 4 5 5.7550e-005 0.0012 3.6794e-006 Sim Sim Sim 8.6987e-011 3.2255e-009 7.9255e-012 Sim Sim Sim [20,20] 3 5 10 5.1022e-005 3.8777e-005 2.2346e-005 Sim Sim Sim 6.2544e-011 5.3112e-011 4.5657e-011 Sim Sim Sim [50,50] 5 10 25 4.7528e-005 3.7291e-005 0.0098 6.9866e-004 Sim Sim Sim Sim 7.4723e-011 5.5490e-011 1.5982e-008 1.9603e-009 Sim Sim Sim Sim [300,300] 3 5 50 100 150 7.6576e-005 4.4346e-005 0.1083 0.3974 0.7522 Sim Sim Sim Sim Sim 9.3832e-011 5.3672e-011 3.4155e-007 1.2600e-006 2.0534e-006 Sim Sim Sim Sim Sim De uma forma geral, a prestação dos modelos foi semelhante aos resultados obtidos utilizando a PSF de dimensão [350,350], mantendo-se a tendência de resultados ligeiramente mais próximos do desejado quando recorremos a raios de dimensão ímpar e 36 4.1. Interpretação geométrica da SVD Embora a Decomposição em Valores Singulares se aplique a matrizes reais e complexas, a sua interpretação geométrica é melhor compreendida se apenas considerarmos o domínio real. Consideremos como [2] uma esfera de raio unitário. A sua imagem, por uma transformação de uma qualquer matriz mxn é uma hiperelipse. No domínio , uma hiperelipse pode ser descrita como uma superfície obtida por um alongamento de uma esfera de raio unitário, em determinadas direções ortogonais. Os vetores {σi,ui} são os semieixos principais da hiperelipse, com comprimentos σ1,…, σm. A esfera de raio unitário, usualmente definida como a esfera euclidiana no espaço n, ou seja, a esfera de norma 2, é representada por S. Assim, a hiperelipse AS resulta da imagem de S sob a transformação A (figura 23). 4.2. Definição de SVD Seja A uma matriz real mxn (vamos supor que m ≥ n por uma questão de descrição mais simples) então existem 2 matrizes, U (matriz real mxm) e V (matriz real nxn), ortogonais (UTU = Imxm e VTV = Inxn) tais que: Figura 23 Decomposição em Valores Singulares de uma matriz 2x2 [1] 37 (5) Onde os σi são iguais ou maiores que zero e a matriz diagonal é de dimensão mxn. Os valores σ1, σ2, … , σn chamam-se valores singulares de A e costumam estar ordenados de forma decrescente σ1 ≥ σ2 ≥ … ≥ σn > 0 [2]. As colunas de U são os vetores singulares à esquerda e as colunas de V são os vetores singulares à direita de A. Chama-se decomposição SVD de A à igualdade , onde é uma matriz diagonal com valores singulares iguais a σ1, σ2, … , σn. A partir de (5) vem (6) A relação entre os vetores singulares à direita e à esquerda pode ser descrita da seguinte forma: (7) Ou seja, (8) Sendo, A Demonstração da existência e unicidade subjacente à Decomposição em Valores Singulares pode ser encontrado em [2, capítulo 4]. Pode ainda descrever-se esta relação como uma equação matricial, a qual se denomina por Decomposição em Valores Singulares Reduzida: 38 Se em (8) considerarmos apenas as n primeiras colunas de U, formando a matriz , mxn, vem   (9) (mxm)(nxn) (mxn)(nxn) Esta expressão pode ainda ser estruturada da seguinte forma:  (10) Embora esta forma simplificada da SVD tenha a sua utilidade, a verdade é que o recurso à sua forma completa ainda é muito importante. As colunas de Û têm n vetores ortonormais de dimensão m num subespaço . A menos que m = n, os vetores não formam uma base de e, por conseguinte Û também não é uma matriz quadrada. No entanto, se acrescentarmos m – n colunas ortonormais, então é possível converter Û numa matriz quadrada unitária. Se for feito um ajuste em Û, então também é necessário realizar um ajuste em  [2]. Uma das utilidades da decomposição em valores singulares é que revela a característica da matriz A: o número de elementos não nulos, r, de Σ é a característica da matriz. Esta característica é o número de linhas ou colunas da matriz linearmente independentes [3]. Outra utilidade da SVD refere-se ao cálculo da norma de uma matriz, cujo valor é dado pelo maior destes valores singulares [3]. 4.3. Decomposição em Valores Singulares Truncada O uso da SVD na sua forma completa será porventura a forma mais correta de proceder à reconstituição de uma imagem contudo, é possível recorrer a uma truncatura dos valores singulares e obter uma imagem relativamente próxima da original. Numa primeira abordagem, a ideia de utilizar apenas parte da informação de uma matriz pode suscitar algumas dúvidas quanto à qualidade do resultado, ou seja, se nós, por exemplo, comprimirmos uma imagem não guardando a totalidade dos seus valores 39 singulares, esperamos que quando a formos analisar novamente, ela não terá a qualidade inicial. A resposta a esta dúvida prende-se com um aspeto muito importante de uma matriz, a sua característica, cujo cálculo nos permite saber o número máximo de linhas ou colunas linearmente independentes da matriz. Uma imagem cuja informação foi truncada pode permanecer inalterada se nós considerarmos toda a informação relativa à sua característica. Como seria previsível, se nós realizarmos uma truncatura mais pequena, a nossa imagem perde detalhes. Senão vejamos. Se tivermos uma matriz X, de dimensão N (mxn), nós podemos definir um k, menor que a característica de matriz X, cuja reconstituição da imagem a partir de Xreconst =A-1b, seja dada por [1]: ∑ Sendo, [ ][ ] [ ] ∑ Mais à frente neste trabalho, veremos algumas aplicações desta técnica. 4.4. Cálculo da Decomposição em Valores Singulares Inicialmente, a SVD não era muito utilizada por problemas de precisão. Os valores singulares estão relacionados com os valores próprios da matriz A*A (sendo A de dimensão mxn, com m> n). Quando se aplica a SVD a A*A (ou AA*, se fosse n> m) existe sempre uma grande perda da precisão na determinação dos valores singulares mais pequenos [3]. O primeiro algoritmo, considerado estável, para o cálculo da SVD, foi desenvolvido por Golub e Kahan. Ambos sugeriram a redução da matriz à forma bidiagonal, com uma transformação Householder baseada no método de Lanczos, como o passo inicial do cálculo da SVD. 40 Neste método, os valores e os vetores singulares podem ser calculados como valores próprios e vetores próprios de uma determinada matriz tridiagonal, segundo um método baseado numa sequência de Sturm. Golub e Reinsch sugeriram ainda uma adaptação do método de Factorização QR para o cálculo da matriz bidiagonal, ao invés da tridiagonal dada pelo método de Lanczos [3]. Embora existam outros algoritmos esta adaptação continua a ser, ainda hoje, a mais utilizada. 4.4.1. Cálculo da SVD a partir dos valores próprios A SVD de uma matriz pode ser calculada a partir da Decomposição em Valores Próprios de uma matriz simétrica correspondente A*A. Para entender este processo comecemos por considerar o seguinte teorema [3, capítulo 6]: Teorema. Os valores singulares não-nulos de uma matriz A são as raízes quadradas dos valores próprios não-nulos de A*A (ou AA*). Demonstração. Temos que Então é semelhante a , e portanto possuem os mesmos n valores próprios. Sendo assim, é compreensível que os valores próprios da matriz diagonal sejam com n-r valores próprios nulos adicionais quando n>r [2] . Considere-se agora uma matriz A , com m ≥ n, sendo a sua Decomposição em Valores Singulares dada por . Sabe-se que . A SVD de A pode ser calculada partindo dos pressupostos do teorema anterior, ou seja (cf [2]): 1. Determina-se C = A*A. 2. Recorrendo ao algoritmo QR, calculam-se os vetores próprios de C, que arrumamos nas colunas da matriz V e os valores próprios de C que arrumamos na diagonal de lambda, ou seja, C = VΛV*. 3. Considera-se como sendo uma matriz, m x n, que contém as raízes quadradas dos elementos de Λ. 4. Resolve-se o sistema Σ*U* = V*A*, equivalente a UΣ = AV, para obter U. 41 Apesar de se recorrer frequentemente a este tipo de algoritmo, a verdade é que este método acarreta um problema no que se refere à estabilidade, uma vez que a transformação de um problema de SVD num problema de Decomposição em Valores próprios faz com que a sensibilidade a perturbações aumente significativamente [4] [5]. 4.4.2. Cálculo da SVD sem calcular A*A e usando a Bidiagonalização de Golub Kahan Uma forma alternativa de reduzir a SVD a um problema de Valores Próprios de uma matriz quadrada consiste em substituirmos A por uma matriz quadrada com ela relacionada de maior dimensão. Considere-se (seguindo a referência 2) H uma matriz Hermítica (no caso de matrizes reais H será simétrica) de dimensão 2m x 2m, dada por: [ ] (11) Sendo, e Estas igualdades podem ser rearranjadas para substituir em (11): [ ][ ] [ ][ ] (12) Esta igualdade corresponde justamente à Decomposição em Valores Próprios de H. Os valores singulares de A representam os valores absolutos dos valores próprios de H. Os vetores singulares de A podem ser obtidos a partir dos vetores próprios de H. Vamos agora resolver o problema de valores próprios de (12) começando por o transformar num problema equivalente em que a matriz H é substituída por uma bidiagonal. A segunda etapa deste método consiste na aplicação do algoritmo QR. Em seguida, será explicada com maior detalhe cada uma das etapas. Nesta opção 3.2 vamos fazer a bidiagonalização com o algoritmo descrito por Golub e Kahan [2]. 42 4.4.2.1. 1ª etapa - Bidiagonalização A decomposição SVD envolve o cálculo das matrizes U e V, tais que U*AV é uma matriz bidiagonal. O método conhecido como a Bidiagonalização de Golub-Kahan utiliza reflectores de Householder, à direita e à esquerda alternadamente. Cada reflector à esquerda, coloca a zeros uma coluna abaixo da diagonal, e cada reflector à direita coloca a zeros os elementos da linha à direita da supradiagonal correspondente [2]. O processo é exemplificado no esquema seguinte. A matriz U que procuramos é o produto dos reflectores Ui acima definidos, U = U1, …, Un. De modo análogo, temos que V = V1, … , Vn-2. É de realçar que o custo do algoritmo é dado por flops. Figura 24 Bidiagonalização da matriz A [2] 43 4.4.2.2. 2ª Fase - método QR Considere-se B (ver [6]) a matriz bidiagonal obtida do passo anterior, retirando-se os elementos nulos, de forma que [ ] . O problema da decomposição de A fica agora reduzido a uma Decomposição em Valores Singulares de B. Define-se d1, d2,… , dn e f1, f2, … , fn-1 como sendo os elementos da diagonal e da supradiagonal de B, respectivamente. A aplicação ímplicita do algoritmo QR à matriz tridiagonal, T = B*B, é apresentada em seguida. O método QR para calcular os valores próprios de uma matriz B consiste em fazer uma factorização QR dessa matriz T= Q0R0 e em seguida criar uma matriz semelhante T1=R0Q0, que por sua vez será factorizada em QR, T1=Q1R1 e assim sucessivamente até Ti convergir para uma matriz triangular cujos valores diagonais são os que pretendemos. A convergência é acelerada se forem feitas translações de origem e, para manter a aritmética real, convém fazer as translações duplas e implícitas. Para isso, o primeiro passo consiste em calcular os valores próprios de: [ ] A formação explícita da matriz não garante a estabilidade do processo. Assim, surge uma alternativa que consiste em aplicar as rotações de Givens directamente em B. O resultado desta alternativa é apresentado em seguida, conforme [6]. Note-se que, se não fosse pelo sinal +, a matriz seria bidiagonal. No entanto, para a reduzirmos a essa forma, basta determinarmos as rotações de Givens U1, V2, U2, … , Vn-1, Un-1. Todo o processo é seguidamente esquematizado [6]. 44 Uma condição essencial para que se possa aplicar este método pressupõe que a matriz tridiagonal seja irredutível. Quando aparece um elemento nulo na subdiagonal, o problema original fica dividido em dois outros sub-problemas, que envolvem as matrizes B1 e B2, de dimensão menor [6]. [ ] Se para algum , coloca-se o elemento a zero, através das rotações de Givens. Na figura seguinte, apresenta-se uma esquematização do processo, com rotações no plano das linhas (2,3), (2,4) e (2,5) [6]. Figura 25 Redução da matriz à forma bidiagonal recorrendo a rotações de Givens [2] 45 Deste modo, os casos em que ou , para um qualquer k, são facilmente tratados [6]. 4.4.3. Cálculo da SVD sem calcular A*A e usando a Bidiagonalização de LawsonHanson-Chan Uma outra alternativa para o processo consiste em utilizar a Bidiagonalização Lawson-Hanson-Chan (ou LHC) que pressupõe o uso da factorização QR de A, A=QR, onde Q, mxm, é uma matriz ortogonal (se a fosse complexa seria uma matriz unitária) e R, mxn, é uma matriz triangular no sentido de que todos os elementos abaixo da diagonal são nulos. Seguidamente recorre-se à bidiagonalização de R. O custo da factorização QR é de flops e a bidiagonalização de GolubKahan requer flops. Assim, o custo geral do processo é de flops [2]. A LHC só é um processo vantajoso quando requer ou seja para [2]. 4.4.3.1. 2ª Fase - método QR Uma vez que a segunda fase do processo de Decomposição em Valores Singulares é comum para todos os métodos realizados a partir a Bidiagonalização da matriz, também nesta variante recorremos ao método QR como segunda fase do prcesso. Figura 26 Bidiagonalização LHC [2] 52 Ao longo deste trabalho foi possível compreender alguns dos conceitos inerentes ao estudo dos processos de desfocagem e reconstituição de uma imagem e, como tal, compreender a necessidade de possuir informação sobre a causa da desfocagem. Artificialmente é possível desfocarmos e reconstituirmos uma imagem com um elevado grau precisão, uma vez que podemos determinar as matrizes Ac e Ar, responsáveis pela desfocagem. Existem diferentes metodologias que nos permitem desfocar uma imagem, pelo que é possível tentar simular algumas das desfocagens mais comuns na aquisição de imagens, por exemplo, deslocamentos seguindo uma linha, turbulência atmosférica, desfocagem devida ao sistema de entes da máquina fotográfica, ou outro aparelho de aquisição de Capítulo 5. Conclusão e Trabalho Futuro 53 imagem, que aqui chamamos desfocagem propriamente dita e foi a que usamos mais e ocorre quando utilizamos o microscópio eletrónico. O processo de desfocagem engloba um conjunto de variáveis que são definidas de acordo com a situação que queremos simular e que requerem algum conhecimento especializado do modo como as imagens foram adquiridas. Relativamente a algumas imagens médicas de mitocôndrias conseguimos melhorar imagens desfocadas e determinar alguns valores aconselháveis para o raio da matriz da PSF (Point Source Function) que traduz o processo de desfocagem bem como dimensões desta mesma matriz que permitem uma reconstituição ótima da imagem. Determinamos também valores do raio que não vale a pena tentar usar porque dão sistematicamente maus resultados. O uso da SVD truncada como ferramenta alternativa para resolução do modelo linear usado no tratamento de imagens desfocadas é importante no sentido de permitir reduzir a dimensão dos problemas e permite também reconstituir a imagem com boa qualidade e em menos tempo quando as dimensões são grandes, no entanto em alguns casos não foi possível recuperar uma boa imagem e o uso da SVD truncada embora tenha produzido uma imagem reconstruida onde nos parecia ver melhor que na imagem original os organelos, quando analisadas por um especialista foram reprovadas e consideradas artificiais. Previamente tinha sido feito um estudo aprofundado da SVD (Singular Value Decomposition) que permitiu concluir que a SVD feita a partir de duas fases, Bidiagonalização e método QR, é uma forma mais eficiente de proceder à determinação dos valores singulares de uma matriz, do que, como se considera por vezes, o uso de formas simplicistas para o cálculo dos valores singulares. É o caso da decomposição em valores próprios, de ATA, que provoca uma grande variação na estabilidade dos resultados, uma vez que quando estamos perante matrizes cujos valores singulares são muito pequenos ou próximos, esta técnica sofre uma grande perda de precisão. Este trabalho permitiu-nos perceber que o processo de Desfocagem e Reconstituição de uma Imagem é um tema com muito interesse e utilidade em diferentes áreas e muito vasto, mas que por isso mesmo ainda necessita de muita investigação e trabalho. Como trabalho futuro vamos criar novas “functions”, ligadas à construção das matrizes Ar e Ac, de alteração das linhas e colunas da imagem verdadeira para representar 54 a desfocada, no modelo linear, tendo em conta o seu número de condicionamento, de forma a controlarmos melhor o processo de desfocagem e reconstituição. Prevemos ainda o tratamento de outras imagens científicas que possam ter sofrido desfocagens por deslocamento do aparelho de aquisição das mesmas segundo uma linha, horizontal, vertical ou outra. 55 Referências Bibliográficas [1] Hansen, C., Nagy, J. G., O’Leary, D. P.: Deblurring Images: Matrices, Spectra and Filtering, Siam, 2006 [2] Trefethen, L. N., Bau III, D.: Numerical Linear Algebra, Siam, 1997 [3] Björck, Å.: Numerical Methods For Least Squares Problems, Siam, 1996 [4] Stewart, G. W., Matrix Algorithms, Vol. II: Eigensytems, Siam, 2001 [5] Figueiredo, Inês F., Monografia: Técnicas de Cálculo da Decomposição em Valores singulares SVD, 2011 [6] Golub, G H., Van Loan, C. F., Matrix Computations, 3ª Edição, Siam, 1996 [7] N. Muller, L. M., B. M. Herbst.: Singular value decomposition, eigenfaces, and 3D reconstructions, SIAM Review, 46, p. 518-545, 2004. [8] Ashin, R., Morimoto, A., Nagase, M., Vaillancourt, R., Image Compression with Multiresolution Singular Value Decomposition and other methods, Mathematical and Computer Modelling 41 (2005) 773-779, disponível online em http://www.sciencedirect.com/. [9] Cao, Lijie, Singular Value Decomposition Applied To digital Image Processing, Division of Computing Studies, 85212, Arizona State University Polytechnic Campus, Arizona [10] Lu, Jiwen, Zhao, Yongwei, Dominant value decomposition representation for face recognition, (2009), disponível online em www.elsevier.com/locate/sigpro 56 Apêndice A – Programas e Scripts 57 A componente experimental deste trabalho foi toda ela realizada recorrendo ao software MatLab® versão 2009. Os códigos apresentados foram elaborados com recurso à referência [1]. A-1. Implementação de um processo de Desfocagem e posterior Recuperação da Imagem %% Ler imagem I=imread('mitocBem1quarto.JPG'); X=rgb2gray(I); X=im2double(X); figure (1), imshow(X), title ('imagem original preto-branco') [m,n]=size(X); %% construir a PSF for R=3:150; % raio da psf for dim1=350; % dimensão da psf dim=[dim1,dim1]; [PSF, center] = psfDefocus(dim, R); % tipo de desfocagem P=zeros(m,n); P(1:dim1, 1:dim1)=PSF; BC='periodic'; % tipo de condição fronteira % Decomposição de Kronecker para obtenção das matrizes Ac e Ar [Ar, Ac] = kronDecomp(P, center, BC); % Determinação do condicionamento de Ac e Ar cond(Ar) cond(Ac) % criação da matriz desfocada Be= Ac*X*Ar'; figure (2), imshow(Be), title ('imagem desfocada sem adiçao ruido') % recuperar a imagem sem adição ruido t=cputime(); Xreconstruido= Ac\Be/Ar'; tempo=cputime()-t figure (3), imshow(Xreconstruido), title ('imagem desfocada sem adiçao ruido') 58 % introdução do parâmetro ruído invertido for nk=1:10; % valor erro noise=1*10^(-nk) E=noise*rand(m,n); B=Be+ E; % imagem desfocada com adição de ruido figure (4), imshow(B), title ('imagem desfocada + ruído') t=cputime(); Xreconstruido_com_ruido= Ac\B/Ar'; % imagem recuperada sem a subtração do ruído tempo=cputime()-t figure (5), imshow(Xreconstruido_com_ruido), title ('imagem recuperada sem adiçao ruído') t=cputime(); einv=Ac\E/Ar'; % recuperar a imagem com adição de ruido Xreconstruido_sem_ruido=Xreconstruido_com_ruido-einv; tempo=cputime()-t figure (6), imshow(Xreconstruido_sem_ruido), title ('imagem recuperada com adição ruido') % recorrendo à svd [Uc Sc Vc]=svd(Ac); val_sing=diag(Sc); [Ur Sr Vr]=svd(Ar); val_singr=diag(Sr); S1=diag(Sc)*diag(Sr)'; t=cputime(); % imagem recuperada, recorrendo à SVD, sem adição e sem subtraçao erro Xreconstruido_svd =Vc*((Uc'*Be*Ur)./S1)*Vr'; tempo_svd=cputime()-t figure (7), imshow(Xreconstruido_svd), title ('imagem recuperada sem adição ruido') t=cputime(); % imagem recuperada, recorrendo à SVD, sem subtraçao erro 59 Xreconstruido_com_ruido_svd =Vc*((Uc'*B*Ur)./S1)*Vr'; tempo_svd=cputime()-t figure (8), imshow(Xreconstruido_com_ruido_svd), title ('imagem recuperada sem adição ruido e sem subtraçao') t=cputime(); Xreconstruido_sem_ruido_svd =Vc*((Uc'*B*Ur)./S1)*Vr'- Vc*((Uc'*E*Ur)./S1)*Vr'; tempo_svd=cputime()-t figure (9), imshow(Xreconstruido_sem_ruido_svd), title ('imagem recuperada com adição ruido') display('norma XXreconstruido') norm(XXreconstruido) display('norma XXreconstruido_com_ruido') norm(XXreconstruido_com_ruido) display('norma XXreconstruido_sem_ruido') norm(XXreconstruido_sem_ruido) display('norma XXreconstruido_svd ') norm(XXreconstruido_svd) display('norma XXreconstruido_sem_ruido_svd ') norm(XXreconstruido_sem_ruido_svd) display('norma XXreconstruido_com_ruido_svd ') norm(XXreconstruido_com_ruido_svd) pause end end end A-2. Para a reconstrução das imagens recorrendo à SVD truncada foi usado o seguinte script (com a devida atualização dos parâmetros) : %% ler imagem I=imread('mitocBem1quarto.JPG'); X=rgb2gray(I); X=im2double(X); X=X(1:350, 1:350); %% Construir PSF [m,n]=size(X); 60 for R=3; for dim1=6; dim=[dim1,dim1]; [PSF, center] = psfDefocus(dim, R); P=zeros(m,n); P(1:dim1, 1:dim1)=PSF; BC='periodic'; [Ar, Ac] = kronDecomp(P, center, BC); Be= Ac*X*Ar'; [Uc Sc Vc]=svd(Ac); val_sing=diag(Sc); [Ur Sr Vr]=svd(Ar); val_singr=diag(Sr); for k=100; for nk=1:10; noise=1*10^(-nk) E=noise*rand(m,n); S=diag(Sc(1:k, 1:k))*diag(Sr(1:k, 1:k))'; B=Be+ E; Xk=Vc(1:m, 1:k)*((Uc(1:m, 1:k)'*B*Ur(1:n, 1:k))./S)*Vr(1:n, 1:k)'-Vc(1:m, 1:k)*((Uc(1:m, 1:k)'*E*Ur(1:n, 1:k))./S)*Vr(1:n, 1:k)'; imagesc(Xk); axis image; colormap(gray) pause norm(X-Xk) end end end end 61 Apêndice B – Tabelas de Resultados 68 Xrecons_svd Xrecons_svd_com_ruido Xrecons_svd_sem_ruido (recupera?) (recupera?) (recupera?) Xrecons_svd Xrecons_svd_com_ruido Xrecons_svd_sem_ruido DIM da Imagem original 1,00E-01 Sim 0.0624 Não 0.0624 Não 0.0936 4,21E-11 3,19E+04 5,93E-11 1,00E-02 \ \ Não 0.0624 Não 0.1248 \ 3,19E+03 4,52E-11 Condições Fronteira: 1,00E-03 \ \ Não 0.0468 Não 0.0936 \ 324,513 4,275E-11 Periodic' 1,00E-04 \ \ Não 0.0312 Não 0.0936 \ 30,412 4,978E-11 1,00E-05 \ \ Sim 0.0624 Sim 0.0936 \ 3,366 4,844E-11 1,00E-06 \ \ Sim 0.0780 Sim 0.0936 \ 0.3673 5,07E-11 1,00E-07 \ \ Sim 0.0624 Sim 0.1248 \ 0.0341 5,40E-11 1,00E-08 \ \ Sim 0.0312 Sim 0.1248 \ 0.0034 4,68E-11 1,00E-09 \ \ Sim 0.0624 Sim 0.1248 \ 4,04E-04 5,47E-11 1,00E-10 \ \ Sim 0.0312 Sim 0.0936 \ 3,50E-05 4,59E-11 Dimensão PSF Raio Cond Ar Cond Ac Erro 566,106 566,106 5 diferença norma (XTempo(s) Tempo(s) Tempo(s) Xreconst Xrecons_com_ruido Xrecons_sem_ruido (recupera?) (recupera?) (recupera?) Xreconst Xrecons_com_ruido Xrecons_sem_ruido DIM da Imagem original 10 1,00E-01 Sim 0.0624 Não 0.0468 Sim 0.0624 2,77E-08 7,35E+06 4,18E-08 1,00E-02 \ \ Não 0.0312 Sim 0.0312 \ 1,37E+06 1,81E-08 Condições Fronteira: 1,00E-03 \ \ Não 0.0624 Sim 0.0312 \ 5,06E+04 3,55E-08 Periodic' 1,00E-04 \ \ Não 0.0312 Sim 0.0624 \ 1,83E+04 2,70E-08 1,00E-05 \ \ Não 0.0312 Sim 0.0624 \ 1,865E+03 4,14E-08 1,00E-06 \ \ Não 0.0312 Sim 0.0780 \ 94,517 3,55E-08 1,00E-07 \ \ Não 0.0468 Sim 0.0624 \ 13,296 1,71E-08 1,00E-08 \ \ Sim 0.0624 Sim 0.0624 \ 1,0715 4,16E-08 1,00E-09 \ \ Sim 0.0312 Sim 0.0312 \ 0.1548 2,97E-08 1,00E-10 \ \ Sim 0.0624 Sim 0.0624 \ 0.0049 2,57E-08 1,43E+04 1,43E+04 diferença norma (XTempo(s) Tempo(s) Tempo(s) Dimensão PSF Raio Cond Ar Cond Ac Erro Xrecons_svd Xrecons_svd_com_ruido Xrecons_svd_sem_ruido (recupera?) (recupera?) (recupera?) Xrecons_svd Xrecons_svd_com_ruido Xrecons_svd_sem_ruido DIM da Imagem original 10 1,00E-01 Sim 0.0624 Não 0.0624 Não 0.1248 1,13E-08 7,35E+06 2,15E-08 1,00E-02 \ \ Não 0.0624 Não 0.0936 \ 1,37E+06 1,27E-08 Condições Fronteira: 1,00E-03 \ \ Não 0.0312 Não 0.0936 \ 5,06E+04 1,29E-08 Periodic' 1,00E-04 \ \ Não 0.0624 Não 0.0936 \ 1,83E+04 2,17E-08 1,00E-05 \ \ Não 0.0624 Não 0.1560 \ 1,865E+03 1,462E-08 1,00E-06 \ \ Não 0.0624 Não 0.1404 \ 94,517 1,096E-08 1,00E-07 \ \ Não 0.0624 Não 0.0936 \ 13,296 1,653E-08 1,00E-08 \ \ Sim 0.0624 Sim 0.1248 \ 1,0715 1,95E-08 1,00E-09 \ \ Sim 0.0936 Sim 0.1560 \ 0.1548 7,12E-09 1,00E-10 \ \ Sim 0.0624 Sim 0.1248 \ 0.0049 1,23E-08 Dimensão PSF Raio Cond Ar Cond Ac Erro diferença norma (XTempo(s) Tempo(s) Tempo(s) 1,43E+04 1,43E+04 69 Xreconst Xrecons_com_ruido Xrecons_sem_ruido (recupera?) (recupera?) (recupera?) Xreconst Xrecons_com_ruido Xrecons_sem_ruido DIM da Imagem original 175 1,00E-01 Sim 0.0624 Não 0.0624 Sim 0.0624 9,00E-08 1,033E+07 4,36E-07 1,00E-02 \ \ Não 0.0624 Sim 0.0624 \ 1,07E+06 1,08E-07 Condições Fronteira: 1,00E-03 \ \ Não 0.0468 Sim 0.0624 \ 1,03E+05 9,36E-08 Periodic' 1,00E-04 \ \ Não 0.0624 Sim 0.0624 \ 1,12E+04 9,57E-08 1,00E-05 \ \ Não 0.0468 Sim 0.0624 \ 1,07E+03 9,95E-08 1,00E-06 \ \ Não 0.0312 Sim 0.0624 \ 1,02E+02 9,27E-08 1,00E-07 \ \ Não 0.0780 Sim 0.0624 \ 10,348 9,45E-08 1,00E-08 \ \ Sim 0.0624 Sim 0.0624 \ 1,0672 8,91E-08 1,00E-09 \ \ Sim 0.0312 Sim 0.0624 \ 0.1070 9,38E-08 1,00E-10 \ \ Sim 0.0624 Sim 0.0624 \ 0.0109 9,29E-08 diferença norma (XTempo(s) Tempo(s) Tempo(s) Dimensão PSF Raio Cond Ar Cond Ac Erro 6,37E+03 6,37E+03 Xrecons_svd Xrecons_svd_com_ruido Xrecons_svd_sem_ruido (recupera?) (recupera?) (recupera?) Xrecons_svd Xrecons_svd_com_ruido Xrecons_svd_sem_ruido DIM da Imagem original 175 1,00E-01 Sim 0.0624 Não 0.0312 Não 0.0936 3,25E-08 1,033E+07 3,40E-08 1,00E-02 \ \ Não 0.0624 Não 0.1560 \ 1,07E+06 3,26E-08 Condições Fronteira: 1,00E-03 \ \ Não 0.0624 Não 0.0936 \ 1,03E+05 3,20E-08 Periodic' 1,00E-04 \ \ Não 0.0312 Não 0.1248 \ 1,12E+04 3,35E-08 1,00E-05 \ \ Não 0.0624 Não 0.0936 \ 1,07E+03 3,30E-08 1,00E-06 \ \ Não 0.0624 Não 0.0936 \ 1,02E+02 3,17E-08 1,00E-07 \ \ Não 0.0624 Não 0.1248 \ 10,348 3,18E-08 1,00E-08 \ \ Sim 0.0312 Sim 0.1248 \ 1,0672 3,22E-08 1,00E-09 \ \ Sim 0.0312 Sim 0.1248 \ 0.1070 3,19E-08 1,00E-10 \ \ Sim 0.0936 Sim 0.0936 \ 0.0109 3,15E-08 6,37E+03 6,37E+03 Dimensão PSF Raio Cond Ar Cond Ac Erro diferença norma (XTempo(s) Tempo(s) Tempo(s) Xreconst Xrecons_com_ruido Xrecons_sem_ruido (recupera?) (recupera?) (recupera?) Xreconst Xrecons_com_ruido Xrecons_sem_ruido DIM da Imagem original 180 1,00E-01 Sim 0.0468 Não 0.0312 Sim 0.0624 0.0086 1,782E+12 0.4429 1,00E-02 \ \ Não 0.0312 Sim 0.0468 \ 1,83E+11 0.0494 Condições Fronteira: 1,00E-03 \ \ Não 0.0624 Sim 0.0624 \ 1,42E+10 0.0102 Periodic' 1,00E-04 \ \ Não 0.0624 Sim 0.0312 \ 2,91E+09 0.0132 1,00E-05 \ \ Não 0.0312 Sim 0.0624 \ 1,37E+08 0.0101 1,00E-06 \ \ Não 0.0468 Sim 0.0312 \ 1,80E+07 0.0128 1,00E-07 \ \ Não 0.0312 Sim 0.0312 \ 2,038E+06 0.0106 1,00E-08 \ \ Não 0.0312 Sim 0.0624 \ 3,1577E+05 0.0117 1,00E-09 \ \ Não 0.0468 Sim 0.0624 \ 3,05E+04 0.0174 1,00E-10 \ \ Não 0.0468 Sim 0.0624 \ 3,10E+03 0.0206 diferença norma (XTempo(s) Tempo(s) Tempo(s) Dimensão PSF Raio Cond Ar Cond Ac Erro 6,35E+06 6,35E+06 70 Xrecons_svd Xrecons_svd_com_ruido Xrecons_svd_sem_ruido (recupera?) (recupera?) (recupera?) Xrecons_svd Xrecons_svd_com_ruido Xrecons_svd_sem_ruido DIM da Imagem original 180 1,00E-01 Sim 0.0312 Não 0.0312 Não 0.0936 0.0080 1,782E+12 0.0097 1,00E-02 \ \ Não 0.0312 Não 0.1092 \ 1,83E+11 0.0077 Condições Fronteira: 1,00E-03 \ \ Não 0.0624 Não 0.1248 \ 1,42E+10 0.0075 Periodic' 1,00E-04 \ \ Não 0.0312 Não 0.0936 \ 2,91E+09 0.0068 1,00E-05 \ \ Não 0.0312 Não 0.0936 \ 1,37E+08 0.0072 1,00E-06 \ \ Não 0.0312 Não 0.0936 \ 1,80E+07 0.0070 1,00E-07 \ \ Não 0.0624 Não 0.1248 \ 2,038E+06 0.0059 1,00E-08 \ \ Não 0.0780 Não 0.1248 \ 3,1577E+05 0.0071 1,00E-09 \ \ Não 0.0468 Não 0.0936 \ 3,05E+04 0.0069 1,00E-10 \ \ Não 0.0624 Não 0.1248 \ 3,10E+03 0.0063 6,35E+06 6,35E+06 Dimensão PSF Raio Cond Ar Cond Ac Erro diferença norma (XTempo(s) Tempo(s) Tempo(s) 71 Imagem Original: Raio Dimensao PSF intervalo k Melhor aproximação k dif original intervalo k Melhor aproximação k dif intervalo k Melhor aprox dif [350,350] 3 [6,6] [1,100] ≈100 ≈1.7 [101,200] ≈200 ≈1.5 [201,300] ≈300 ≈1.2 [8,8] [1,100] ≈100 ≈2 [101,200] ≈200 ≈2 [201,300] ≈300 ≈1.5 [20,20] [1,100] ≈100 ≈1.99 [101,200] ≈200 ≈2 [201,300] ≈300 ≈1.5 [100,100] [1,100] ≈100 ≈1.99 [101,200] ≈200 ≈2 [201,300] ≈300 ≈1.5 [200,200] [1,100] ≈100 ≈2 [101,200] ≈200 ≈2 [201,300] ≈300 ≈1.5 [350,350] [1,100] ≈100 ≈3 [101,200] ≈200 ≈2 [201,300] ≈300 ≈1.5 4 [8,8] [1,100] ≈100 ≈3 [101,200] ≈200 ≈2 [201,300] ≈300 ≈1.5 [10,10] [1,100] ≈100 ≈3 [101,200] ≈200 ≈2 [201,300] ≈300 ≈1.5 [20,20] [1,100] ≈100 ≈3 [101,200] ≈200 ≈2 [201,300] ≈300 ≈1.5 [100,100] [1,100] ≈100 ≈3 [101,200] ≈200 ≈2 [201,300] ≈300 ≈1.5 [350,350] [1,100] ≈100 ≈3 [101,200] ≈200 ≈2 [201,300] ≈300 ≈1.5 5 [10,10] [1,100] ≈100 ≈3 [101,200] ≈200 ≈2 [201,300] ≈300 ≈2.5 [12,12] [1,100] ≈100 ≈4 [101,200] ≈200 ≈2 [201,300] ≈300 ≈2.5 [50,50] [1,100] ≈100 ≈4 [101,200] ≈200 ≈2 [201,300] ≈300 ≈2.5 [100,100] [1,100] ≈100 ≈4 [101,200] ≈200 ≈2 [201,300] ≈300 ≈2.5 [350,350] [1,100] ≈100 ≈4 [101,200] ≈200 ≈2 [201,300] ≈300 ≈2.5 6 [12,12] [1,100] ≈100 ≈4 [101,200] ≈200 ≈3 [201,300] ≈300 ≈2.4 [50,50] [1,100] ≈100 ≈4 [101,200] ≈200 ≈3 [201,300] ≈300 ≈2.4 [100,100] [1,100] ≈100 ≈4 [101,200] ≈200 ≈3 [201,300] ≈300 ≈2.4 [350,350] [1,100] ≈100 ≈4 [101,200] ≈200 ≈3 [201,300] ≈300 ≈2.4 7 [14,14] [1,100] ≈100 ≈3.8 [101,200] ≈200 ≈3 [201,300] ≈300 ≈2.5 [50,50] [1,100] ≈100 ≈3.7 [101,200] ≈200 ≈3 [201,300] ≈300 ≈2.5 [100,100] [1,100] ≈100 ≈3.7 [101,200] ≈200 ≈3 [201,300] ≈300 ≈2.5 [350,350] [1,100] ≈100 ≈3.7 [101,200] ≈200 ≈3 [201,300] ≈300 ≈2.5 Recuperação da Imagem recorrendo à SVD truncada 72 Dimensão Imagem Raio Dimensao intervalo k Melhor aprox k dif Raio Dimensao intervalo k Melhor aprox k dif Original [500,500] 5 [10,10] [1,500] ≈70 ≈13 6 [12,12] [1,500] ≈65 ≈14 [12,12] [1,500] ≈70 ≈13 [15,15] [1,100] ≈60 ≈14 [15,15] [1,500] ≈65 ≈13 [20,20] [1,100] ≈60 ≈14 [20,20] [1,500] ≈65 ≈13 [50,50] [1,100] ≈60 ≈14 [50,50] [1,500] ≈65 ≈13 [100,100] [1,100] ≈60 ≈14 [100,100] [1,500] ≈65 ≈13 [200,200] [1,100] ≈58 ≈14 [200,200] [1,500] ≈65 ≈13 [500,500] [1,100] ≈58 ≈14 [500,500] [1,500] ≈65 ≈13 B-3. Resultados obtidos para a Imagem de mitocôndrias desfocada Dimensão Imagem Raio Dimensao intervalo k Melhor aprox k dif Raio Dimensao intervalo k Melhor aprox k dif Original [500,500] 3 [6,6] [1,500] ≈100 ≈10 4 [8,8] [1,500] ≈100 ≈10 [7,7] [1,500] ≈100 ≈10 [9,9] [1,500] ≈100 ≈10 [8,8] [1,500] ≈100 ≈10 [10,10] [1,500] ≈100 ≈10 [9,9] [1,500] ≈100 ≈10 [12,12] [1,500] ≈80 ≈11 [10,10] [1,500] ≈100 ≈10 [15,15] [1,500] ≈85 ≈11 [12,12] [1,500] ≈100 ≈10 [20,20] [1,500] ≈86 ≈11 [15,15] [1,500] ≈100 ≈10 [50,50] [1,500] ≈86 ≈11 [20,20] [1,500] ≈100 ≈10 [100,100] [1,500] ≈86 ≈11 [50,50] [1,500] ≈100 ≈10 [200,200] [1,500] ≈90 ≈11 [100,100] [1,500] ≈100 ≈10 [500,500] [1,500] ≈90 ≈11 [200,200] [1,500] ≈100 ≈10 [500,500] [1,500] ≈100 ≈10