Elipsoides xeneralizados activos: aplicación á segmentación de imaxes médicas 3D
Full text
UNIVERSIDADE DE SANTIAGO DE COMPOSTELA ELIPSOIDES XENERALIZADOS ACTIVOS: APLICACIÓN Á SEGMENTACIÓN DE IMAXES MÉDICAS 3D RAQUEL DOSIL LAGO NOVEMBRO, 2000 DEPARTAMENTO DE ELECTRÓNICA E COMPUTACIÓN
Ós meus pais
Agradecementos Quero expresar o meu sincero agradecemento a tódalas persoas que contribuíron á realización desta tesiña. En primeiro lugar ó director desta memoria, Xosé Manuel Pardo López, pola súa colaboración, pero sobre todo polo interese e a confianza que depositou no meu traballo. A Diego Cabello e Antonio Mosquera por contar comigo, pois o apoio recibido por ambos facilitou enormemente o meu labor neste grupo de investigación. A Víctor Leborán e a Fernando Vilariño por inverter tanto tempo en axudarme a poñer esto en marcha. A Manolo Penedo pola información facilitada e as suxestións sobre algúns aspectos do traballo. A Paula López pola axuda que me mandou dende Italia a cerca dos algoritmos xenéticos. A Víctor Brea, Juanjo e Nando polos consellos e por darme ánimos e, en xeral, a tódolos compañeiros deste departamento, que realmente conseguen que viva un ambiente estupendo.
i ÍNDICE 1 INTRODUCCIÓN..................................................................................................... 1 2 RECONSTRUCCIÓN DE IMAXES 3D................................................................. 5 2.1 Procesamento de imaxes médicas........................................................................ 6 2.2 Segmentación baseada en modelos...................................................................... 7 2.3 Inicialización de superficies deformables.......................................................... 11 2.4 Modelo de forma global: o superelipsoide......................................................... 12 2.5 Recoñecemento de baixo nivel .......................................................................... 13 3 MODELOS DE SUPERFICIE............................................................................... 17 3.1 Esquemas de representación .............................................................................. 17 3.1.1 Distribución de puntos de fronteira .......................................................... 17 3.1.2 Ecuacións implícitas................................................................................. 19 3.1.3 Ecuacións explícitas ................................................................................. 20 3.1.4 Representacións “de varrido”................................................................... 21 3.1.5 Técnicas de subdivisión do espacio.......................................................... 22 3.1.6 Transformación do eixe medio................................................................. 25 3.1.7 Representacións cualitativas..................................................................... 26 3.1.8 Xeometría constructiva: CSG................................................................... 27 3.1.9 Esquemas de representación comparados................................................. 28 3.2 Métodos de reconstrucción ................................................................................ 30 3.2.1 Posta en correspondencia mediante transformacións ríxidas ................... 30 3.2.2 Métodos deformables................................................................................ 32 3.2.3 Técnicas de optimización ......................................................................... 39 4 SUPERCUÁDRICAS.............................................................................................. 45 4.1 Medidas de erro.................................................................................................. 46 4.2 Xeneralización para posición e orientación arbitrarias...................................... 48 4.3 Ambigüidade na representación mediante superelipsoides ............................... 49 4.4 Inicialización dos parámetros de transformación ríxida .................................... 50 4.5 Deformacións globais......................................................................................... 53 4.6 Optimización con algoritmos xenéticos............................................................. 55 4.6.1 Codificación.............................................................................................. 56 4.6.2 Rendemento do AX. ................................................................................. 57 4.6.3 Mecanismos de selección ......................................................................... 58 4.6.4 Función de adaptación.............................................................................. 59 4.7 Posta en correspondencia con modelos a priori................................................. 59 5 EXTRACCIÓN DE PARCHES DE SUPERFICIE ............................................. 61
ii 5.1 Filtros multiescala recursivos............................................................................. 62 5.1.1 Filtros multiescala..................................................................................... 62 5.1.2 Implementación recursiva dos filtros ....................................................... 64 5.2 Detección de puntos de fronteira ....................................................................... 66 5.3 Cálculo das curvaturas locais............................................................................. 67 5.4 Extracción de parches de superficie................................................................... 70 5.4.1 Agrupamento de puntos en parches de superficie .................................... 71 5.4.2 Selección de parches de superficie........................................................... 72 5.5 Observacións...................................................................................................... 75 6 RESULTADOS........................................................................................................ 79 6.1 Optimización da función implícita..................................................................... 79 6.2 Preprocesado...................................................................................................... 82 6.2.1 Detección de puntos de fronteira.............................................................. 82 6.2.2 Detección de parches de superficie .......................................................... 83 6.2.3 Etiquetado e selección de parches de superficie....................................... 83 6.3 Posta en correspondencia................................................................................... 84 7 CONCLUSIÓNS E PRINCIPAIS APORTACIÓNS........................................... 99 BIBLIOGRAFÍA........................................................................................................ 101
1 1 INTRODUCCIÓN A importante mellora que están a experimentar as imaxes médicas tridimensionais de alta resolución, está derivando nun maior desenvolvemento das técnicas de reconstrucción de estructuras anatómicas, que cada vez poden ofrecer resultados de mellor calidade. A extracción da xeometría a partir das imaxes médicas serve de apoio a actividades clínicas dende a diagnose ata a planificación cirúrxica e simulación do comportamento fisiolóxico. Pero a pesar dos recentes avances nas tecnoloxías de adquisición de imaxes, que proporcionan vistas excepcionais da anatomía interna, o uso do computador para cuantificar e analizar as estructuras involucradas con precisión e eficiencia é aínda limitada. A dificultade na análise das imaxes médicas está en dotar ós sistemas de reconstrucción de obxectos 3D da suficiente flexibilidade para tratar coa complexidade e variabilidade das formas anatómicas entre distintos individuos. Ademais, os inconvenientes típicos dos datos mostreados e o ruído provocan a miúdo que os contornos das estructuras resulten indistinguibles e desconectados. Durante os últimos anos, os modelos deformables [28] fixéronse moi populares en recoñecemento e delineación de obxectos, particularmente no campo das imaxes médicas [37]. O esquema básico dun modelo deformable consiste nunha curva, se traballamos en 2D, ou superficie, en 3D, que se deforma localmente ata axustarse á fronteira do obxecto máis próximo. Os modelos deformables 3D foron inicialmente propostos por Terzopoulos en [56]. Aquel primeiro modelo consistía nun esquema de elementos finitos guiado por técnicas dinámicas. Posteriores traballos foron contribuíndo a darlles flexibilidade. Exemplos da aplicación dos modelos deformables poden atoparse en [34, 38, 47, 50, 52, 53]. Cando os modelos deformables son aplicados a datos 3D existen dúas modalidades de segmentación: 2D (corte a corte) e 3D. En ambos casos o modelo inicial debe estar o suficientemente próximo á fronteira do obxecto que se quere segmentar para que o axuste sexa correcto. Unha das vantaxes da segmentación corte a corte é que se pode partir dunha inicialización semiautomática ou manual dun corte. Despois de segmentar ese corte mediante un modelo deformable 2D utilízase o contorno resultante como modelo inicial no corte veciño, repetindo o proceso para tódolos cortes da imaxe. Esta inicialización baséase na proximidade espacial entre cortes adxacentes, que asegura contornos próximos en forma e localización. O principal inconveniente é que non hai outra interacción entre contornos de cortes veciños no proceso de deformación, a parte da inicialización, o que pode levar a inconsistencias na superficie resultante. Se o contorno se perdese nun corte en particular, o método fallaría. Ademais, require unha etapa de postprocesado para conecta-la secuencia da contornos e obter a superficie continua. O uso de verdadeiros modelos de superficie podería conducir a métodos de segmentación máis robustos, xa que aseguraría a coherencia e a suavidade global nos
8 CAPÍTULO 2. RECONSTRUCCIÓN DE IMAXES 3D técnica top-down, é dicir, baseada en modelos. Partindo dunha descrición de alto nivel representada por un modelo definido a priori, e despois dunha fase de localización do mesmo na imaxe, aplicaríase un proceso de deformación local que ofrece unha descrición de máis baixo nivel. A outra tarefa é a do preprocesado da imaxe, necesaria para a posta en correspondencia co modelo. O enfoque que se lle deu a este problema foi o dun proceso bottom-up, como se verá no capítulo 1. Partindo dos datos da imaxe, vaise organizando a información ata obter unha representación que permita decidir que puntos da imaxe poden pertencer ó obxecto e cales non. A posta en correspondencia realizase entre o modelo a priori estes e puntos. Cabe pensar que esta fase de preprocesado podería aproveitarse directamente para obter unha descrición do modelo de forma global, en lugar de servirse dela simplemente para a posta en correspondencia cun modelo definido de antemán. Sen embargo, imos ver nesta sección que o procesado bottom-up non é unha técnica robusta da que se poida obter unha descrición de alto nivel fiable. A continuación imos facer unha comparación das características que presentan as técnicas baseadas en imaxes e as baseadas en modelos, para concluír que as primeiras son, polo xeral, as máis adecuadas para calquera tipo de aplicación. Como se dixo antes, uns dos problemas que afronta a análise de imaxes son o recoñecemento de formas e a reconstrucción de imaxes para visualización en 3D, onde os datos cos que se vai traballar representan un nivel de intensidade de gris ou ben compoñentes de cor. A segmentación de baixo nivel non baseada en modelos traballa unicamente sobre esta información, analizando as propiedades locais da imaxe para dividir a escena en obxectos. Esto non se pode acadar sen facer ningún tipo de suposición sobre as características da imaxe. Por exemplo, se se fai unha segmentación baseada en rexións, estarase establecendo a hipótese de que áreas ou volumes de distinta cor pertencen a obxectos distintos. Tamén podemos realizar unha detección de bordes, na que se asume que a variación do nivel de gris se corresponde cunha fronteira entre diferentes obxectos se supera un determinado umbral. O mesmo se pode facer operando no dominio da frecuencia, atribuíndolle ás baixas frecuencias a información de bordes ou superficies fronteirizas. Todas estas aproximacións son perfectamente válidas cando as hipóteses nas que están baseadas se cumpren. O feito é que estas case nunca se verifican de xeito estricto. Son moitos os efectos que se están obviando ó establecer as hipóteses anteriores, como as variacións de contraste ou iluminación en distintas partes dunha imaxe, a presencia doutras estructuras ademais das de interese, as texturas na superficie dos obxectos, as oclusións e, en xeral, calquera tipo de distorsión ou de degradación que poida sufrir unha imaxe no proceso de captura. As técnicas non baseadas en modelos non teñen a capacidade de manexar información de tipo estructural. Esto non só vai repercutir na calidade dos resultados, senón que ademais, dificulta a extracción de propiedades globais dos obxectos. Por exemplo, se aplicamos unha detección de bordes a unha imaxe ruidosa ou con texturas, o resultado vai ser un conxunto de segmentos de curva ou de superficie, segundo traballemos en dúas ou tres dimensións, inconexos e ramificados, e non a fronteira completa dun obxecto, e só dese obxecto, diferenciándoo así do resto da imaxe. Polo tanto, tampouco se van poder extraer propiedades xeométricas do obxecto, como posición, área ou
2.2 Segmentación baseada en modelos 9 volume, etc., e menos aínda, ir máis alá e pretender compara-lo obxecto con outros coñecidos e identificalo como membro dunha clase. Existen técnicas que, como se mencionou ó principio desta introducción, pretenden extraer información da estructura global das imaxes, tentando identifica-los obxectos que contén partindo dos resultados obtidos no procesado de baixo nivel, ás que chamabamos técnicas guiadas por datos. Hai desenrolados métodos, coma o que se presenta no libro de Ballard e Brown [1], que, partindo dunha detección de puntos de fronteira, tentan conectar os fragmentos de superficie resultantes para formar superficies pechadas que definan obxectos. O inconveniente que teñen estes métodos é que parten xa de resultados que poden ser erróneos, o que desvirtúa todo o proceso posterior. Por eso estas técnicas nunca van te-la robustez e a fiabilidade das técnicas baseadas en modelos. Tódalas consideracións anteriores son trabas que xorden a causa de traballar a nivel local sen empregar modelos. Hai outro aspecto que se pode cuestionar nesta clase de técnicas, e é que non se fai uso de coñecemento a priori sobre os obxectos, nin predeterminando certas propiedades, nin impoñendo limitacións, baseándose en hipóteses sobre o contido das imaxes. Ben é certo que, como se comentou anteriormente, sempre que se procesa unha imaxe estase a realizar algún tipo de suposición sobre a información que contén. Sen embargo, as premisas sobre as que se leva a cabo esta modalidade de segmentación non teñen en conta ningún tipo de información específica da clase particular de imaxe que se está a analizar, non sendo o axuste dalgún posible umbral de nivel de gris ou de gradiente de nivel de gris. Será raro o dominio de aplicación do tratamento de imaxes no que non se teña ningunha noción sobre a súa natureza ou do seu contido. Podemos descoñecer que é ou como é o que a imaxe representa exactamente, pero en moitas ocasións vai ser posible intuír certas características. Por exemplo, se os obxectos teñen contornos pechados ou abertos, se as súas superficies son planas ou curvas. Nalgunhas aplicacións incluso se sabe de antemán que é o que estamos vendo na imaxe. Tanto a falta de visión de conxunto, como o desaproveitamento do coñecemento a priori, como a imposibilidade de exercer ningún tipo de control sobre as propiedades xeométricas dos obxectos, poñen os métodos anteriores en desvantaxe con respecto ás técnicas baseadas en modelos, que por definición son o método axeitado para tratar con este tipo de información. O feito de poder controlar ou prefixar determinados parámetros da representación dunha superficie, implica poder forzar a que os obxectos se adapten ás características desexadas e tamén introducir coñecemento. A definición e o modo de utilización dun modelo vai depender do tipo de aplicación. Se o que queremos é obter a reconstrucción dun obxecto, o modelo vai consistir nun esquema de representación que caracterice ó obxecto de forma apropiada e nun método de reconstrucción para obte-los parámetros do esquema a partir da imaxe. Nestes casos, aínda o modelo non impón ningunha restricción sobre a forma concreta do obxecto, si pode limitar os graos de liberdade na representación da forma. Esto sucede cando, por exemplo, queremos modelar un obxecto mediante unha función implícita cuadrática. Os únicos parámetros que podemos variar son radios e os parámetros de transformacións afíns. Un esquema de representación sempre vai limitar en certa medida a xeometría dun obxecto, pero tamén se poden controlar determinadas propiedades no proceso de
10 CAPÍTULO 2. RECONSTRUCCIÓN DE IMAXES 3D reconstrucción. Un exemplo claro desta situación dáse no modelos de snakes ou contornos activos. O modelo non é tan só o esquema de representación, que non deixa de ser unha simple distribución de puntos no espacio, senón que incorpora unha técnica de reconstrucción que permite guiar a evolución dos puntos ata que se adaptan á fronteira do obxecto, definindo forzas entre os puntos, como se fosen partículas físicas. En aplicacións de recoñecemento o enfoque é distinto. Pártese dun modelo a priori, e entón o que hai que implementar é un método de busca que explore imaxe para atopar obxectos de características similares á primitiva dada. Deste xeito soluciónanse as limitacións do procesado de baixo nivel, pois podemos identifica-lo obxecto como un caso particular do modelo, ou sexa realizar unha clasificación dos obxectos que aparecen na escena. Ademais, podemos obter medidas das súas propiedades, que veñen dadas polo modelo de forma explícita, na súa parametrización, ou de forma implícita. Neste tipo de tarefas é preciso ter un certo coñecemento do tipo de obxectos que van aparecer na imaxe e dispor dun modelo a priori de cada un deles. Existen tamén aplicacións nas que, a pesar de teren como obxectivo a descrición detallada da forma dun obxecto, vaise precisar tamén dunha fase de recoñecemento dese obxecto nunha escena, que é o caso que comentabamos anteriormente, sobre o que versa este traballo. A maneira máis común de enfoca-lo problema é partir dun modelo a priori, que representa un prototipo do obxecto. Con el realizarase unha primeira etapa de recoñecemento do obxecto. Logo permítese a deformación do modelo para obter finalmente o modelo particular desa imaxe. Este é o caso do sistema de segmentación 2D desenvolto por Altman e Brinkley [10], que emprega un modelo de forma promedio, obtida a partir de imaxes de adestramento. O proceso de axuste consiste en realizar unha estimación dos puntos do contorno da imaxe cada vez máis refinada a medida que se van introducindo novos puntos do modelo a priori. A estimación realízase mediante un filtro de Kalman en función da información de correlación entre modelo e imaxe. En ocasións, os modelos a priori e a posteriori van vir representados mediante esquemas distintos. Na metodoloxía que máis adiante presentamos, esto vai ser inevitable, posto que o modelo de supercuádricas ofrece unha descrición de forma moi global. Esto faino moi axeitado como modelo a priori, pero moi pobre como descrición local, polo que se, cando sexa precisa unha fase de refinamento vai haber que realizar un cambio de representación. No caso de utilizar o modelo de supercuádricas como inicialización dun modelo deformable, unha simple mostraxe da superficie será suficiente. Á hora de expresa-las restriccións impostas, o tipo de información que hai que manexar vai ser, usualmente, aínda que non necesariamente, de baixo ou alto nivel en aplicacións de reconstrucción ou recoñecemento respectivamente. Podemos, por exemplo, pretender limitar a rugosidade dunha superficie nun proceso de reconstrucción, e entón estaremos operando a baixo nivel, ou podemos querer encontrar obxectos de forma aproximadamente esférica en tarefas de recoñecemento, e aquí estaremos utilizando información de máis alto nivel. O esquema de representación que escollamos debe ser capaz de parametrizar a información coa que vaiamos traballar. En consecuencia, a representación escolleita é determinante na construcción dun bo método de modelado. Ademais, condicionará ou limitará os posibles métodos de
2.3 Inicialización de superficies deformables 11 reconstrucción ou de busca. Para o obxectivo que nós propoñemos, a elección de supercuádricas parece ser a máis adecuada. 2.3 Inicialización de superficies deformables Aínda que a nosa pretensión é que o sistema sexa aplicable a calquera tipo de imaxes 3D, o dominio de aplicación que o motiva é o modelado de estructuras anatómicas, concretamente de estructuras óseas. Así, observarase que as imaxes coas que traballamos son todas deste tipo. A idea de facer este tipo de aproximación no ámbito da imaxe médica, xorde como consecuencia dun problema que presentan os modelos de deformación local. Os modelos deformables, coñecidos tamén como modelos activos, ofrecen un gran nivel de detalle nos resultados, que a penas se ven afectados polo ruído presente na imaxe, posto que permiten controlar a suavidade da forma do contorno ou superficie, de aí que sexa unha das técnicas máis empregadas neste campo. A forma de reconstruír o modelo é mediante a evolución da súa forma guiada pola información de gradiente de nivel de gris. As ecuacións que rexen a evolución dos puntos de fronteira resultan de resolve-lo problema de minimización dunha funcional que representa a enerxía asociada á fronteira. Entre os termos que contribúen a esa enerxía, está o orixinado polas forzas que atraen ós puntos cara ós máximos de gradiente. Este algoritmo de reconstrucción iterativo tende a estabilizarse en mínimos locais da enerxía, é dicir, os puntos de fronteira máis próximos ó punto de partida. Deste xeito, se na imaxe aparece máis dun obxecto, o resultado pode non ser o desexado. A superficie ou contorno debe estar inicialmente moi próxima á do obxecto de interese para que outras estructuras non interfiran. Ata o de agora unha das formas máis sinxelas de solucionar esto é mediante a interacción co usuario, ó que se lle permite inicializar manualmente o modelo. Como xa se comentou na introducción, o método de inicialización de cada corte por separado, deseñando manualmente o contorno inicial no primeiro corte e propagando os resultados ós cortes veciños, non proporciona resultados totalmente satisfactorios. Outra opción é definir unha primitiva simple, como por exemplo unha esfera, e permitir ó usuario elixir a posición do seu centro e tal vez o radio. Tampouco esta é, en xeral, unha boa solución en imaxe médica, pois a presencia de gran número de tecidos no volume esixe unha inicialización mellor [37]. Estas consideracións sobre os modelos deformables obrigan a pensar nunha maneira de inicializa-lo contorno activo, que neste caso será unha superficie pois traballaremos en 3D, mediante un modelo a priori. Nalgúns estudios [29] realizados óptase por obter estes modelos mediante a reconstrucción manual dun primeiro prototipo de obxecto mediante unha representación de distribución espacial de puntos, como unha triangulación, que representaría un modelo promedio dese tipo de obxecto. Nós consideramos que esto é innecesario, posto que a obtención e almacenamento deste modelo a priori non compensa se despois se vai someter a unha deformación, que é a que decide a forma final. Pensamos que unha descrición a groso modo é suficiente para a inicialización, especialmente no caso de estructuras óseas, e que resulta máis sinxela
12 CAPÍTULO 2. RECONSTRUCCIÓN DE IMAXES 3D de obter e manexar, ademais de que ten unha parametrización moito máis reducida. Esto levounos á elección dos modelos de elipsoides xeneralizados. 2.4 Modelo de forma global: o superelipsoide O esquema de representación que se escolleu para crear o modelo de forma global foi o de supercuádricas, tamén chamadas elipsoides xeneralizados ou superelipsoides. Vén descrito por unha función implícita non lineal. Os motivos desta escolla son dous: o bo compromiso que ofrece entre potencia de representación e parametrización reducida e a invariancia ante transformacións ríxidas. Esta representación emprega un número reducido de graos de liberdade, o que aporta facilidade de manexo e eficiencia no almacenamento en memoria. Esta simplicidade non se obtén a costa de reducir en exceso a variedade de formas que pode adoptar, pois ofrece unha potencia de representación aceptable para o dominio de aplicación, o que o fai moi atractivo para o noso propósito. Para aumenta-la potencia de representación destas primitivas, se así fose necesario, é posible aplicar sobre este modelo unha serie de transformacións de forma global de modo moi sinxelo (como se fixo neste traballo) obtendo así unha flexibilidade de representación moito maior, sen incrementar en exceso o número de parámetros do modelo. Por outro lado, o vector de parámetros do modelo, ademais dos parámetros de forma, debe incorporar variables que den información sobre a posición, a orientación e o tamaño con que aparece na imaxe. No proceso de obtención do modelo, ben sexa para recoñecemento ou para reconstrucción, vai haber que determinar estes parámetros. Esto vaise facer buscando a transformación ríxida que aplicada ó modelo proporcione a mellor aproximación ó obxecto. Por este motivo, o esquema de representación debe ser invariante a este tipo de transformacións, que engloban translación, rotación e escalado, que tamén se chaman transformacións ríxidas porque non alteran a forma do modelo. Non tódolos esquemas presentan invariancia a estas transformacións, un motivo máis a favor dos superelipsoides, que si teñen esta propiedade. Dispoñendo dun modelo a priori para un determinado obxecto podemos localizar dito obxecto nunha imaxe mediante a posta en correspondencia entrambas (ver figura 2.1). Tanto para a obtención do modelo de partida (figura 2.1, 1º paso), como para a súa posterior posta en correspondencia co obxecto na imaxe (figura 2.1, 4º paso), o xeito de calcular os parámetros do modelo que se axusta adecuadamente ós puntos de superficie é realizando un proceso de optimización, no que se tenta minimizar unha función que depende da distancia de cada punto do obxecto ó punto máis próximo da superficie da primitiva. No primeiro caso habería que obter tódolos parámetros do modelo de obxecto, tanto os que describen a súa forma, coma os de transformación ríxida. Para a posta en correspondencia os parámetros de forma están prefixados e só hai que determina-la transformación ríxida. Para extraer os puntos pertencentes á superficie do obxecto é necesario realizar un preprocesado que elimine da imaxe outros puntos de fronteira que poidan aparecer (pasos 2º e 3º), aínda que esto pódese evitar empregando técnicas de crecemento de modelos. Na aplicación que aquí se describe recórrese á primeira das dúas opcións. O crecemento de modelos o que fai é inicializar un gran número de superelipsoides por
2.5 Recoñecemento de baixo nivel 13 todo o volume e, a continuación, entrar nun proceso iterativo que fai evolucionar a algúns deles ata adaptarse ás formas existentes na imaxe, mestres que os demais van sendo eliminados. Como resultado aparecerán segmentados tódolos obxectos presentes na imaxe, e non só o de interese. Para poder adaptar esto ás necesidades da aplicación sería necesario realizar unha selección do modelo correcto, empregando para elo coñecemento a cerca da súa xeometría. Coa outra estratexia, o coñecemento aplícase antes do modelado, na selección dos puntos de fronteira. Temos entón que o método de construcción dos modelos vai ser o cálculo do mínimo dunha función non lineal nun espacio multidimensional. De entre o amplo conxunto de técnicas que se poden empregar para a optimización, a elixida foi un algoritmo xenético. Esta elección ten varias xustificacións, relacionadas principalmente coa converxencia a mínimos globais e coa súa adecuación a funcións complexas. Neste caso a complexidade radica na non linearidade e na gran dimensionalidade do problema, que non afectan en absoluto á codificación do vector de parámetros no algoritmo xenético, e tamén no manexo de restriccións no rango dos parámetros, que se tratan dun xeito moi sinxelo. De todos modos ó algoritmo de optimización non vai ser un factor crítico no resultado final. Lembremos que o obxectivo non é unha descrición detallada do obxecto. 2.5 Recoñecemento de baixo nivel O procesado da imaxe previo ó modelado ou á posta en correspondencia ten como obxectivo a extracción dun conxunto de puntos que, presumiblemente, pertenzan á superficie do obxecto. Para determinar esta distribución de puntos pódense utilizar técnicas de procesado de baixo nivel non baseado en modelos. O resultado que vai ofrecer, como cabe esperar dese tipo de tratamento da imaxe, será moi pobre, polos motivos que mencionamos no apartado 2.2, pero a etapa de modelado será capaz de extraer a información global implícita nesa nube de puntos. O preprocesado consta de dúas fases (ver figura 2.1). En primeiro lugar hai que facer unha detección de bordes (figura 2.1, 2º paso). Para esta fin pódense aplicar filtros no dominio da frecuencia, que actúan eliminando a información de baixa frecuencia, ou pódense empregar filtros espaciais, que detectan variacións locais do nivel de gris. Dentro desta última modalidade, existen varias estratexias. Por un lado están aquelas baseadas nunha estimación das primeiras derivadas direccionais do nivel de gris, calculadas para cada punto da imaxe nunha certa veciñanza. Como combinación das distintas compoñentes calculadas obtense unha magnitude que nos informa do brusca que é a variación da cor en cada punto, que xeralmente será unha aproximación ó gradiente no dominio discreto. A continuación hai que definir un criterio que, a partir desta información, determine que puntos pertencen á fronteira dun obxecto. Un criterio é definir un umbral, de forma que tódolos puntos que teñan valor de gradiente maior ca el considérase que pertencen a unha fronteira. A imaxe que resulta de tomar esta opción pode presentar superficies de fronteira de grosor maior que a unidade. Unha solución mellor, que foi a empregada neste traballo, é determinar os máximos locais do gradiente de nivel de gris. Por outro lado están as técnicas baseadas no cálculo das derivadas segundas direccionais, que permiten realizar unha aproximación á laplaciana.
14 CAPÍTULO 2. RECONSTRUCCIÓN DE IMAXES 3D Os puntos de fronteira son nesta ocasión aqueles nos que a función laplaciana cruza polo valor cero. A segunda fase consiste no etiquetado das superficies resultantes da detección de bordes, mediante o cómputo das propiedades que as caracterizan. Estas etiquetas servirán logo para realiza-la selección de superficies pertencentes ó obxecto, baseándonos en coñecemento (figura 2.1, 3º paso). En función do tipo de obxecto a segmentar, definirase unha función de avaliación de cada parche de superficie que indica en que medida se parece cada parcela de superficie ó prototipo da clase de obxecto. Esta función dependerá de propiedades xeométricas, como a curvatura, que caracteriza a forma dun modelo determinada, e dos niveis de gris e de gradiente de nivel de gris que, en imaxes médicas, poden axudar a determinar o tipo de tecido do órgano. Os parches que tomen un valor da función maior son bos candidatos a pertencer á superficie do obxecto a modelar Á final desta etapa o que temos é un conxunto de parches de superficie con alta probabilidade de pertencer ó obxecto, pero que non describen correctamente a súa forma, pois aparecerán rexións inconexas e outras ramificadas pola unión con parches pertencentes a outros obxectos. Sen embargo, estes datos son bos dabondo para realizala posta en correspondencia co modelo, que será o que realmente utilizaremos como descrición da superficie, e non os parches de fronteira, resolvendo así o problema da carencia de información estructural.
2.5 Recoñecemento de baixo nivel 15 Figura 2.1. Secuencia de procesado para a inicialización dun modelo deformable. No paso 1 defínese un modelo de forma a priori a partir da superficie promedio do obxecto extraída de imaxes de adestramento. Os pasos 2 e 3 constitúen a fase de preprocesado de baixo nivel, da que se obteñen os puntos de fronteira para a posta en correspondencia co modelo a priori, no paso 4. IMAXE ORIXINAL 2. EXTRACCIÓN DE PUNTOS DE FRONTEIRA CONXUNTO DE PARCHES DE SUPERFICIE 3. ETIQUETADO E SELECCIÓN DE PARCHES PARCHES DO OBXECTO A MODELAR 4. POSTA EN CORRESPONDENCIA MODELO A POSTERIORI MODELO A PRIORI 1. MODELADO MEDIANTE SUPERCUÁDRICAS PROTOTIPO DE SUPERFICIE DO OBXECTO
3 MODELOS DE SUPERFICIE 3.1 Esquemas de representación Á hora de elixir un esquema de representación de superficie debemos ter en conta varios aspectos. En primeiro lugar hai que considerar se é apropiada para o tipo de imaxes que imos tratar. Neste sentido haberá que estudia-las seguintes propiedades: • potencia; qué variedade de formas á capaz de adoptar e se estas son suficientes no dominio da aplicación. • completitude; buscarase un método que non presente ambigüidades, é dicir, que dous obxectos diferentes non teñan a mesma representación. • unicidade; que a un obxecto só lle corresponda unha única representación. A definición formal destes conceptos pódese encontrar no traballo de Flynn e Jain [24]. En segundo lugar hai que ter en conta outras cuestións máis prácticas, relacionadas co custe computacional da obtención do modelo, cos requirimentos de memoria e coa facilidade de implementación dos métodos de reconstrucción. A continuación preséntase unha revisión dalgunhas das técnicas de representación máis comúns. Para cada unha delas indícanse as súas propiedades, que quedan resumidas na táboa 3.1 á final deste apartado. Tamén se estudia unha característica moi importante en aplicacións de modelado, que é a invariancia ante transformacións ríxidas. 3.1.1 Distribución de puntos de fronteira Describir unha superficie ou un volume por medio dunha nube de puntos ten a vantaxe de que a súa obtención a partir dunha imaxe é bastante directa. Sen embargo, a información aportada soamente polas coordenadas dos puntos é insuficiente á hora de visualizar o resultado ou de aplicar un postprocesado. Haberá que dotar á representación de estructura, incluíndo información acerca da conectividade entre os puntos. Por este motivo, os datos organízanse xeralmente en forma de grafo [24]. Segundo o tipo de aplicación será interesante organizar a información dun xeito ou doutro. Pódese estructurar de forma xerárquica, por exemplo nunha árbore, de tal xeito que os nodos folla sexan os puntos, e capas superiores da xerarquía emparellarían nodos para formar arestas, agruparían arestas para formar polígonos, polígonos para construír poliedros e incluso se podería operar sobre poliedros simples para obter figuras máis complexas. Haberá ocasións nas que a complexidade destas representacións será innecesaria. Para estas situación hai técnicas máis sinxelas, como é a coñecida como wireframes (figura
24 CAPÍTULO 3. MODELOS DE SUPERFICIE determina-lo plano óptimo. Por outra parte, este método non soluciona o problema da non invariancia fronte a transformacións ríxidas, co agravante de que, a diferencia dos octrees, a árbore non é única para unha imaxe dada. Esto é moi fácil de ver. Se temos unha imaxe negra cunha figura branca, teremos infinitos planos que cortan a imaxe en dúas rexións e que deixan a figura totalmente contida nunha delas, quedando a outra baleira, mentres que na árbore de octantes as divisións están prefixadas. Marching Cubes A técnica de Marching Cubes, máis que un esquema de representación, é un método de reconstrucción de superficies trianguladas, polo que se podería incluír esta técnica no primeiro desta serie de apartados, xunto coa triangulación de Delaunay. Sen embargo, ten a particularidade de que fai uso dunha voxelización da imaxe para dividir o espacio en rexións, a partir do cal extrae a triangulación. O método consiste na extracción dunha isosuperficie a partir dunha imaxe volumétrica. A isosuperficie é o análogo a unha curva de nivel, pero estendida a 3D. Neste caso vén definida por un nivel de gris, que marca a fronteira do obxecto, así que, dalgún xeito esta é unha técnica de detección de puntos de fronteira baseada en rexións, xa que require dunha umbralización. Tamén a árbore de octantes se pode estender deste xeito, pois as imaxes definidas en niveis de gris hai que transformalas a binarias. O mecanismo é o seguinte. Supoñamos, simplificando un pouco as cousas, que temos un volume onde o espacio baleiro aparece coloreado cun nivel de gris en torno a i, e que ese volume contén un obxecto de cor media j. Na zona de transición haberá unha variación do nivel de gris de i a j que será máis ou menos suave. Entre eses dous valores podemos definir un umbral de nivel de gris que establece unha división entre puntos internos e externos ó obxecto. Este é un proceso de baixo nivel, con tódalas (a) (b) (c) Figura 3.6. Quadtree
3.1 Esquemas de representación 25 limitación que esto supón. Do que se trata agora e de obter unha representación da isosuperficie en forma de triangulación a partir desa clasificación dos puntos. A idea é que, se dous puntos adxacentes están situados en rexións distintas, é dicir, se un pertence ó interior do obxecto e o outro é externo, entón a isosuperficie debe pasar por entre os dous. O que fai o algoritmo de Marching Cubes é, en primeiro lugar, agrupar os puntos en cubos (voxels). Para cada cubo estudiarase a que rexión pertence cada un dos seus 8 vértices. A continuación haberá que determinar cal é o plano ou planos que dividen o cubo de forma que os puntos internos e externos queden situados en semiespacios separados. Se por exemplo temos un único vértice pertencente ó volume do obxecto, entón o plano de división será un que corte ás arestas do cubo que comparten ese vértice. A intersección deste plano e as arestas do cubo prodúcese por tres puntos, definindo un triángulo da superficie do obxecto. Este é só un dos 15 posible casos que se poden dar, contando co caso trivial, onde tódolos vértices están na mesma rexión. A triangulación correspondente a cada caso pódese ver na figura 3.7. Agora só resta decidir por que punto da aresta vai a cruzar cada plano de división. A solución máis sinxela é cortar sempre polo punto medio da aresta. Esta aproximación moi poucas posibilidades para as orientacións dos triángulos, polo que a superficie que se obtén resulta moi pouco suave, especialmente se a imaxe se toma cunha frecuencia de mostraxe baixa. Outra aproximación, máis usada que a anterior, é realizar unha interpolación lineal dos niveis de gris sobre os puntos das arestas do cubos. O punto de corte será o que tome o nivel de gris umbral que se definiu para a clasificación das rexións. Deste xeito obtéñense superficies máis suaves. 3.1.6 Transformación do eixe medio Como comentabamos, a representación dun obxecto mediante voxelización, indicando para cada punto da imaxe se pertence ó volume do obxecto ou non, é moi custoso no almacenamento. Algunhas das técnicas que vimos conseguen comprimi-la representación, substituíndo a voxelización punto a punto por outra máis compacta, onde os voxels son rexións do espacio de tamaño variable. A técnica de transformación do eixe medio (MAT, Medial Axis Transform) ,[48 ,49], permite obter resultados aínda mellores, pois un volume vai quedar reducido ó seu esqueleto. O esqueleto dunha superficie é unha liña,ou conxunto de liñas que recolle a súa estructura. A esta liña é que se lle chama eixo medio. Se o que temos é un volume, o esqueleto será unha superficie, polo que ,en ocasións, farase referencia ó obxecto medio (medial object), en lugar do eixo medio, para englobar ámbalas dúas definicións. Para determina-lo esqueleto dun obxecto hai que definir unha métrica que proporcione, para cada punto do interior do volume, unha medida da distancia ó punto do contorno da superficie ou superficie do volume máis próximo. Unha vez determinado o mapa de distancias, extráese a liña que une os máximos e ese será o esqueleto. Intuitivamente vese que esta liña vai equidistar do contorno da superficie, de aí o nome transformada de eixo medio. A representación MAT consiste no esqueleto e no valor da distancia asociado a cada un dos seus puntos. Por tanto, a métrica empregada vai ser moi importante para a posterior reconstrucción da figura. Algunhas das medidas que se soen utilizar son a distancia euclídea, a distancia Manhattan, tamén chamada city block, que representa a suma das compoñentes do vector diferencia, e a distancia tomada como o
26 CAPÍTULO 3. MODELOS DE SUPERFICIE máximo de entre estas compoñentes. O esqueleto calculado a partir da distancia euclídea pódese interpretar como o lugar xeométrico dos centros das circunferencias, no caso 2D, inscritas no contorno, que con radio máximo, percorren o interior da superficie. En 3D pódese definir igual só que agora son esferas. Esto é moi fácil de ver mediante un exemplo, coma o da figura 3.8, que mostra o esqueleto tanto para dúas como para tres dimensións. As dúas últimas métricas son aproximacións que, en certas aplicacións poden ser suficientes. Supoñen un aforro en tempo de cálculo. En calquera dos tres casos a relación entre unha figura e a súa representación MAT é sempre válida e biunívoca, ou sexa, que a representación é completa e única. Ó mesmo tempo, a potencia de representación é moi alta, pois permite definir formas arbitrarias con gran resolución, e todo esto nunha representación compacta que require moi pouco espacio para o seu almacenamento. 3.1.7 Representacións cualitativas As representacións cualitativas manexan un tipo de información totalmente distinto ó dos esquemas anteriores. Non se parametrizan magnitudes características dunha figura xeométrica, senón que se sobe no nivel de abstracción para intentar facer unha descrición cualitativa do aspecto que presenta. Aínda que o tratamento cuantitativo é necesario en moitas das áreas da visión artificial, no caso concreto da realización de tarefas de análise ou recoñecemento pode ser suficiente establecer unha distinción entre clases de obxectos empregando un conxunto de descriptores das súas propiedades, relacionados coa percepción que temos deles. A idea foi proposta por Biederman, quen elaborou unha clasificación do que el chamou geons (GEometrical iONS), diferenciándose cada un deles polos valores dun conxunto de propiedades cualitativas. En tarefas de análise esta idea pode ser interesante, sempre que o número de clases de obxectos definidos sexa suficiente no campo de aplicación. En caso contrario, a clasificación resulta ambigua, pois nunha mesma clase estamos englobando obxectos distintos que nos interesa diferenciar. En xeral este tipo de representación soe ser incompleta. Figura 3.7 As 15 posibles divisións do cubo
3.1 Esquemas de representación 27 3.1.8 Xeometría constructiva: CSG CGS son as siglas de Constructive Solid Geometry. Esta representación fai uso dalgunhas das que vimos anteriormente, a modo de primitivas de volume. Obxectos máis complexos confecciónanse a partir das primitivas aplicando sobre elas transformacións xeométricas e/ou combinando varias figuras mediante operacións de unión, intersección ou diferencia. Para crear un obxecto pode ser necesaria toda unha xerarquía de operacións sobre as primitivas, que se estructuran en forma de árbore binaria (figura 3.9). As primitivas de forma que se empregan soen vir definidas mediante ecuacións implícitas, incluídos algúns volumes de revolución como cilindros e conos, e tamén poliedros, pero estes últimos, en lugar de representárense mediante unha distribución de vértices, veñen dados polo volume encerrado por un conxunto de planos. Deste xeito, todas estas primitivas se poden descompoñer en subárbores, nos que os nodos folla son expresións do tipo ( ) 0,, ≥zyxf ás que se lles chama semiespacios, pois definen unha superficie que divide o espacio en dúas rexións. Se f é lineal teremos un poliedro, e se é de segundo grao, teremos as cuádricas. Substituíndo na árbore CSG a subárbore correspondente a cada primitiva, poderemos obter formas complexas como combinación de semiespacios. Tamén existen sistemas deste tipo que, ademais das primitivas anteriores, empregan superficies spline pechadas e incluso NURBS. Este esquema de representación está pensado para aproveitar a potencia de representación das distintas primitivas que emprega. Nese sentido a solución é mellor que calquera das anteriores, pero a costa dunha serie de inconvenientes moi importantes. Un é que o emprego de semiespacios e de operacións booleanas fan desta unha representación moi descritiva, volvéndose dificultoso o cálculo de propiedades xeométricas. O outro é que a árbore CSG non é única para un obxecto determinado. Este sistema é completo e válido na medida en que o son as primitivas de superficie empregadas. O feito de combinar estas primitivas non impón ningún tipo de ambigüidade adicional nin da lugar representacións erróneas. Figura 3.8. Obxecto medio en 2D e 3D
28 CAPÍTULO 3. MODELOS DE SUPERFICIE 3.1.9 Esquemas de representación comparados Botando unha ollada ás características que presentan os distintos esquemas de representación descritos nesta revisión podemos construír a táboa 3.1, onde se especifican para cada esquema as propiedades que se mencionaban ó principio desta sección. A modo de resume preséntase a continuación unha comparativa entre os distintos esquemas, facendo referencia ás propiedades de potencia de representación, completitude e unicidade. Sobre a potencia de representación, pódese dicir que os esquemas baseados en nubes de puntos teñen unha potencia de representación tan alta como se desexe, en función da densidade superficial de puntos que se elixa. Sen embargo, esta representación estará limitada a figuras poliédricas. O mesmo problema presentan os esquemas de subdivisión do espacio. Dende este punto de vista, as ecuacións explícitas son máis potentes, pois é posible escoller o grao das funcións polinómicas, o que permite definir superficies con curvaturas tan suaves como se crea conveniente. Tamén a representación MAT ten unha potencia de representación alta. No extremo oposto están as ecuacións implícitas, que soen ter un número moi limitado de graos de liberdade. Respecto ós GEONS, o simple feito de ser unha representación cualitativa xa supón unha forte restricción en canto á capacidade descritiva, pero eso dependerá de como se definan as clases de obxectos. Finalmente, temos os sistema CSG que, como herdan as propiedades das primitivas que emprega para construír obxectos, é o esquema máis flexible. Da completitude diciamos que se cumpre cando unha determinada configuración dos parámetros do esquema de representación lle corresponde unha única superficie. Esto verifícase en xeral para tódolos esquemas que acabamos de ver, coa excepción dos GEONS, que nunha mesma clase poden agrupar a distintas superficies. De entre os esquemas de representación que acabamos de revisar, só a transformación de eixe medio e os geons cumpren a propiedade de unicidade. Para o resto das técnicas existirán múltiples representacións dun mesmo obxecto. Os esquemas de distribución Figura 3.9. Árbore CSG
3.1 Esquemas de representación 29 de puntos teñen o problema de que existen innumerables formas de determina-la conectividade. Incluso cando se escolle un método sistemático de poligonalización, do tipo da triangulación de Delaunay ou de Marching Cubes, a representación non é única. Por exemplo, se se quere triangular unha superficie plana con forma de cadrado, existen dúas posibles configuracións da conectividade, segundo a diagonal que se elixa para forma-los triángulos. Coas ecuacións implícitas, a non unicidade é debida en xeral ás simetrías de rotación, pero tamén en casos particulares existen outro tipo de ambigüidades, como sucede no caso dos superelipsoides cos parámetros de cuadratura. No caso das representacións mediante ecuacións explícitas, sucede que unha mesma superficie se pode definir empregando distintos conxuntos de puntos de control. Tampouco a construcción da árbore CSG é única para unha superficie determinada, pois un obxecto pódese construír mediante diversas composicións de primitivas simples. A elección dun ou doutro esquema debe ter en conta estes tres aspectos, pero tamén outro de tipo máis práctico, relativos á eficiencia en tempo de cálculo e tamén ó compacta que sexa a representación. Neste sentido, calquera representación que involucre unha descrición en forma de distribución de puntos, de elementos finitos, ou de voxelización do espacio permite realizar certo tipo de tarefas, como por exemplo a visualización, de forma máis rápida. Noutro tipo de representacións haberá que levar a cabo unha discretización da superficie. Sen embargo, non todo son vantaxes. Para obter unha boa resolución na definición das superficies é preciso dispor dun volume de datos moi elevado, co conseguinte uso de espacio en memoria. Ademais, nalgunhas métodos de reconstrucción, o feito de non dispor dunha expresión analítica que defina a superficie pode ser un problema. É o caso das técnicas de rexistro, onde o cálculo das distancias entre superficies e puntos de imaxe non se pode facer de forma directa. Dado que no presente traballo se vai realizar o modelado precisamente mediante un proceso de posta en correspondencia, pensamos que é máis interesante empregar un modelo descrito por unha ecuación implícita. A perda potencia de representación que esto supón, en comparación cos esquemas de nubes de puntos, non é aquí un problema, senón máis ben unha vantaxe, pois a información que se manexa é de máis alto nivel, o que, como veremos, aporta robustez ó proceso. O axuste a pequena escala xa se introduce posteriormente na etapa de deformación local e o cambio de representación que debe mediar entrambas fases efectúase con rapidez por medio dunha simple discretización.
30 CAPÍTULO 3. MODELOS DE SUPERFICIE 3.2 Métodos de reconstrucción Unha vez seleccionado o esquema de representación que se vai empregar para defini-lo modelo, hai que estudiar a estratexia de obtención dos parámetros correspondentes a partir da imaxe. Para a reconstrucción dun modelo existen multitude de técnicas. Nesta sección vanse presentar algunhas das máis comúns, que se clasifican en dous grupos: técnicas de posta en correspondencia con modelos a priori mediante transformacións ríxidas, e modelos de deformacións. En realidade, ambas estratexias son complementarias, como xa se comentou no capítulo introductorio. Para rematar esta revisión presentaranse unha serie de algoritmos de optimización que poden ser aplicados para implementar estes métodos de reconstrucción. 3.2.1 Posta en correspondencia mediante transformacións ríxidas Cando se ten un coñecemento exacto da morfoloxía do obxecto que se desexa modelar, é posible evita-la reconstrucción do modelo para cada imaxe. Partindo dunha conxunto de imaxes de adestramento pódese obter un modelo a priori do obxecto, co que a análise dunha nova imaxe queda reducida á busca da transformación que superpón modelo e obxecto. Estas transformacións deben manter a forma orixinal do modelo a priori, polo que se aplican transformacións ríxidas: de translación, de rotación e de Esquema de Representación Potencia Completitude Unicidade Wireframes alta si non Poligonalización alta si non Triangulación de Delaunay alta si non Cuádricas moi baixa si non Supercuádricas baixa si non Cilindros xeneralizados baixa si non Splines alta si non B-splines alta si non NURBS moi alta si non Árbores de octantes alta si non Marching cubes alta si non MAT alta si si Geons moi baixa non si CGS moi alta - non Táboa 3.1 Características de varios esquemas de representación
3.2 Métodos de reconstrucción 31 escala. A forma de determina-los parámetros da transformación varía segundo o tipo de representación elixido. Besl e McKey [6] presentaron unha revisión completa dos distintos métodos de rexistro de imaxes. A continuación explícanse os métodos utilizados para modelos de superficies paramétricas, definidos mediante ecuacións implícitas, e modelos de distribucións espaciais de puntos. Funcións implícitas Para a reconstrucción de funcións implícitas sóense aplicar técnicas que implican un proceso de optimización dunha función que reflexa a calidade do axuste do modelo ós datos. Unha das aproximacións máis usadas consiste en tratar de minimizar unha función que represente unha medida relacionada coa distancia da superficie ós puntos da imaxe [6, 9]. En xeral, calcular a distancia exacta entre un punto do espacio e unha superficie pode ser complexo e custoso computacionalmente, o que leva a utilización doutras métricas. A propia función implícita, dada por f(x,y,z), pode ser empregada para estimar en que medida se desvían os puntos da imaxe do modelo, sempre e cando a función sexa monotonamente crecente ou decrecente coa distancia exacta. Na revisión de Bolle e Vemuri [9], e tamén no traballo de Whaite de Ferrie [58] preséntanse algunhas aproximacións para os casos de cuádricas e supercuádricas. Tamén con este esquema de representación se pode aplicar un método probabilístico, intentando maximiza-la probabilidade de que a superficie sexa a correcta. Para superficies simples, como cuádricas é sinxelo definir unha función de probabilidade. Por exemplo, para o caso da esfera, onde haberá que determina-los valores óptimos do seu centro c r e do radio R r , preséntase en [9] a seguinte función de probabilidade: () −−−∝ ∑ = N i iRcxp 1 2 2 1 exp r rr σ Para outro tipo de funcións implícitas máis complexas, como conos ou supercuádricas, é máis complicado definir unha función de probabilidade axeitada. Nubes de puntos Cando a modelo non vén dado pon unha función analítica, senón que se representa por unha colección de puntos contidos na superficie, o cálculo das distancias entre puntos do modelo e da imaxe é máis complicado. O máis frecuente é que a correspondencia entre puntos de obxecto e imaxe non sexa coñecida, polo que non se pode aplicar unha técnica de minimización de distancias de forma directa. Haberá que establecer un criterio para emparellar os puntos da superficie do modelo cos do obxecto. Esto é importante non só na posta en correspondencia do modelo co obxecto, senón tamén na definición dos modelos a priori, que xeralmente se obteñen a partir de varias imaxes de proba, que hai que relacionar entre si. Para salvar este problema Krahnstöver [29] e Feldmar e Ayache [22] utilizan unha técnica baseada en marcas (landmarks), que son un conxunto reducido de puntos característicos da superficie. O primeiro dos autores permite a definición destas marcas por parte do usuario. Definidas estas marcas tanto no modelo coma no obxecto, temos unha relación un a un entre elas e podemos calcular distancias entre os pares de puntos. Entón xa se pode aplicar o algoritmo de rexistro de ambas imaxes. Se o modelo se vai
32 CAPÍTULO 3. MODELOS DE SUPERFICIE empregar como superficie inicial dun algoritmo do tipo dos anteriores, basta con determinar a transformación ríxida que hai que aplicar ó modelo para que as marcas queden o máis próximas posible. Se se desexa un maior axuste dos puntos ó modelo, pódese permitir transformacións afíns doutro tipo, ademais das ríxidas, que supoñen unha deformación global do modelo, e incluso outras transformacións máis complexas. Pola súa banda, Feldmar presenta un algoritmo totalmente automático, no que a correspondencia entre puntos de ambas figuras se realiza en función de propiedades xeométricas. Con este obxectivo propón unha técnica na que se van aplicando de forma sucesiva unha serie de transformacións sobre o modelo. En primeiro lugar realízase unha transformación ríxida, onde un conxunto de puntos da superficie do modelo se emparellan con aqueles outros pertencentes ó obxecto que teñan curvaturas similares. Esto vai seguido dunha deformación global do modelo para acadar un maior axuste, aplicando unha transformación afín non ríxida. Á hora de implementa-la busca da transformación afín global, utiliza un criterio de emparellamento co máis próximo, só que a medida da distancia que se utiliza está tomada nun espacio 8D, onde se mide a proximidade entre as coordenadas espaciais, pero tamén entre as curvaturas principais e as compoñentes do vector normal en cada punto. Concretamente utilízase unha medida da distancia euclídea d, onde os eixes das 8 coordenadas están escalados con cadanseu factor αi, de modo que cada termo teña o mesmo peso no cálculo da distancia. A expresión é a seguinte. () ( ) ( ) ( ) ( ) () ()()() 2/1 2 228 2 117 2 6 2 5 2 4 2 3 2 2 2 1 }'''' ''''{', kkkknnnn nnzzyyxxPPd zzyy xx −+−+−+− +−+−+−+−= αααα αααα Finalmente, temos a etapa de deformación local. A diferencia dos modelos de forma global, onde hai que efectuar un cambio de representación, discretizando a superficie, na representación mediante distribución de puntos xa se pode aplicar unha fase de deformación local. Esta podería ser algunha das que máis adiante se exporán nesta sección. No algoritmo de Feldmar óptase polas transformacións afíns, pero esta vez serán independentes para cada punto da superficie. 3.2.2 Métodos deformables Os métodos deformables, sexan de deformacións locais ou globais, precisan da definición dun modelo de partida, a partir da cal se deforman para axustarse ós datos da imaxe imaxe. A inicialización da superficie deformable pode realizarse mediante a interacción co usuario pero, dependendo do modelo de deformación empregado, pode ser necesaria unha maior aproximación inicial. Cando esto é así, é habitual utiliza-lo modelo de deformación como un refinamento dun modelo a priori, que previamente se puxo en correspondencia coa imaxe mediante transformacións ríxidas. Neste apartado preséntanse algunhas técnicas que utilizan deformacións para a reconstrucción detallada dos obxectos. Snakes. A técnica coñecida como modelo de contornos activos -superficies activas no caso 3Dou snake é un método de segmentación moi coñecido pola súa robustez. Entre os
3.2 Métodos de reconstrucción 33 moitos autores que investigaron sobre os modelos de snakes, tal vez os máis destacados sexan Terzopoulos e McInerney [28, 33, 35, 36, 37, 39, 55]. Trátase de simular o comportamento da superficie como un corpo sometido ás leis da física. Para levar esta idea a cabo hai que definir un conxunto de forzas de interacción dos puntos da superficie entre eles e cos puntos da imaxe. A cada par de puntos do mallado da superficie asóciaselle unha forza de atracción definida en función dunha serie de propiedades dos materiais, como elasticidade ou rixidez. Estas interaccións, comunmente chamadas forzas externas ou de imaxe, van selas que controlen as propiedades xeométricas intrínsecas da superficie, regulando a curvatura e a continuidade. As forzas externas, tamén chamadas de imaxe, deben producir unha atracción cara os puntos da imaxe. A superficie final será aquela para a que o total das forzas estean en equilibrio. Se en lugar de empregar a formulación Newtoniana, aplicamos as ecuacións de EulerLagrange, o problema tradúcese na minimización da integral dunha función de densidade de enerxía ou, no dominio discreto, a suma dunha función de enerxía local estendida a tódolos nodos da superficie. A solución de mínima enerxía correspóndese coa posición de equilibrio. As forzas antes citadas derivaríanse desta enerxía: ( ) ( ) ( ) [ ] ∫∫ += S extsnake dudvfEfEE ..int onde () vuff ,= é a superficie, expresada en función dos parámetros naturais [] 1,0, ∈vu . Expliquemos agora máis en detalle os termos que contribúen á enerxía dun nodo. Por unha parte teremos a enerxía responsable da atracción cara a fronteira dos obxectos. Pódense definir de moitas formas, en función das características das imaxes en xogo, pero habitualmente provocan a atracción da superficie cara ós máximos do gradiente, onde se supón que debe haber unha fronteira entre dous obxectos, como por exemplo ( ) ( ) 2 .fIfEext ∇−= Para regula-la continuidade ou a curvatura da superficie introdúcense termos de enerxía interna relacionados coa elasticidade e a flexibilidade da superficie, como se se tratase dunha lámina en tensión ou unha membrana. Unha posible definición de enerxía interna é a seguinte () ( ) ( ) ( ) 22222 2vvuvuuvu fffffffE ++++= βα A variable α controla a elasticidade, a resistencia ó estiramento. Cando α toma o valor cero, unha forza suficientemente grande pode produci-la rotura da membrana. A variable β é a responsable da rixidez e regulará a facilidade coa que se dobra a lámina e ata que punto se pode flexionar. Se β é nulo, estamos permitindo que haxa descontinuidades na orientación da superficie, é dicir, que se formen esquinas. Os parámetros de elasticidade e de rixidez poden ser constantes ou funcións. Os valores de α e de β escóllense en función do tipo de obxectos a reconstruír e da calidade das
40 CAPÍTULO 3. MODELOS DE SUPERFICIE onde o sistema evoluciona a pasos constantes cara ó extremo. En cada iteración do algoritmo de gradiente o vector de parámetros a optimizar x r desprázase unha distancia k, sendo k un valor constante, na dirección do vector de gradiente normalizado, do seguinte xeito ( ) () i i ii xf xf kxx r r rr ∇ ∇ += +1 Esta aproximación presenta un problema que pode ser bastante grave. Se o valor da constante k é grande, pódense producir oscilacións en torno ó extremo sen chegar a acadalo nunca con precisión suficiente. Para eso habería que reduci-lo valor do paso, pero faría o algoritmo moi lento. Unha posible solución a este problema sería permitir que k fose diminuíndo nas sucesivas iteracións do algoritmo. En cada iteración pódese calcular o valor óptimo de k. O valor de k fai a que conduce ó menor valor posible de f na seguinte iteración, determínase pola condición de mínimo: () ( ) () 0 1 r r r rr = ∇ ∇ + ∂ ∂ = ∂ ∂ + i i iii xf xf kxf k xf k Calculado o valor de ki, substitúese na expresión de 1+i x r . Esto diminúe a taxa de converxencia, pero o tempo de cálculo en cada iteración increméntase. Máis eficiente é o algoritmo de Newton-Raphson, que emprega o cálculo de derivadas segundas para acadar un maior axuste do extremo. Derívase empregando unha aproximación en serie de Taylor de primeira orde () ( ) ( ) ( ) iiiii xxxfxfxf r r r r r − ∇ + ≈ ++ 11 e da condición de extremo () 0 0 = ∇xf r. Se o sistema está inicializado nas proximidades do extremo podemos realiza-la seguinte aproximación () ( ) ( ) ( ) ⇒−∇∇+∇==∇ ++ iiiii xxxfxfxf r r r r r r 11 0 ( ) ( ) iiii xfxHxx r r r r ∇−= − + 1 1 Este algoritmo ten o inconveniente de que hai que calcula-la matriz hessiana H en cada iteración. De novo podemos recorrer á suposición de estar preto dun extremo local,o que nos serve para poder realizar unha aproximación de segunda orde da función. co cal a matriz hessiana é constante e non será necesario recalculala de cada vez. Deste xeito, a aplicación deste algoritmo resulta moi interesante para a optimización de funcións non lineais. O principal problema será conducir as variables do sistema ás proximidades do extremo local. Se esto non é posible hai que recorrer á técnica de gradiente descendente, que pola súa parte, non ofrece unha boa aproximación do mínimo. O ideal sería poder empregar a técnica de gradiente descendente nunha etapa inicial do algoritmo e finalmente o algoritmo de Newton-Raphson para un mellor axuste do extremo.
3.2 Métodos de reconstrucción 41 O algoritmo que consegue levar a cabo este propósito é o de Levenberg-Marquardt. É unha variante do algoritmo de Newton-Raphson que introduce un parámetro que regula o a evolución. O que se fai é substituír a matriz hessiana por outra ( ) iiii HA λ + = 1 ijij HA = j i ≠ Inicialmente o parámetro λ toma un valor elevado, de forma que na matriz A dominan os elementos da diagonal, dando lugar a un comportamento similar ó de gradiente descendente. Esto conduce ó sistema ás proximidades do extremo máis próximo. Se o longo da execución do algoritmo se fai decrecer o parámetro, A vai tender a H, e entón teremos de novo o algoritmo de Newton-Raphson, acadando o extremo con precisión. Temperado simulado O temperado simulado (Simulated Annealing, SA) é un bo algoritmo para a obtención de extremos globais, especialmente útil cando a función posúe ademais moitos extremos locais. É apropiada tanto para funcións continuas definidas en espacios multidimensionais, como é o caso que nos interesa, como para funcións avaliadas sobre un espacio de configuracións discreto pero extenso, onde, por un lado, non se poden aplicar algoritmos de gradiente, e por outro, realizar unha busca exhaustiva resulta inabordable. O principio de funcionamento está extraído da termodinámica, concretamente da cristalización das moléculas cando son sometidas a un arrefriamento lento. Dunha forma moi simplificada, pódese describir este fenómeno deste xeito: a perda de enerxía das moléculas ó descender a temperatura produce unha reducción da mobilidade térmica. Se esto sucede a un ritmo suficientemente lento, estas partículas teñen a oportunidade de “bailar” entre distintos estados ata atopar o de mínima enerxía, estabilizándose e formando a estructura cristalina. Para un determinado valor da temperatura, a probabilidade que ten unha partícula de estar nun estado de enerxía E vén dada pola fórmula de distribución de probabilidade de Boltzmann. () − −∝ kT EE Ep min exp Se T decrece lentamente ata cero chégase á situación en que a probabilidade de estar no estado de mínima enerxía será a unidade e o resto serán inaccesibles. Pero para calquera valor de T>0 sempre existe unha certa probabilidade de que a partícula evolucione cara un estado enerxeticamente peor. Vexamos como se aplica esta idea. Se () xfE r = é unha función que queremos optimizar –vaise aplicar ó caso da busca dun mínimo global-, o que faremos será ir xerando valores de x r aleatoriamente en cada iteración do algoritmo. Para cada nova configuración determínase a súa correspondente enerxía e avalíase a probabilidade do sistema de evolucionar cara o novo estado mediante
42 CAPÍTULO 3. MODELOS DE SUPERFICIE () ( ) > − − ≤ = − − − 1 1 1 exp 1 ii ii ii EE kT xExE EE prr En función desa probabilidade resólvese se o salto se fai efectivo ou non. En caso afirmativo substitúese o vector de parámetros anterior polo novo. Observando a función vemos que efectivamente, para un valor de T maior que cero, e a diferencia do que sucede nun algoritmo de gradiente, a probabilidade de pasar a un estado con enerxía maior non é nula, o que vai permitir saír de mínimos locais. Esto é moi útil nas primeiras iteracións do algoritmo, onde se require unha exploración do espacio de configuracións na procura do mínimo global, pero a medida que o sistema evoluciona, esto supón un inconveniente, pois introduce imprecisión á hora de situar con exactitude o extremo. Aquí é onde intervén o parámetro T, que simula o efecto do descenso da temperatura. Inicialmente tómase un valor alto de T, co cal asegúrase que o sistema poida cambiar tanto para mellor como para peor. Despois dunha exploración do espacio onde se estima en que rexión pode estar situado o mínimo global, a reducción paulatina da temperatura fai que diminúa a probabilidade de que a enerxía suba. Finalmente case tódalas iteracións se resolverán a favor do descenso, cambiando o comportamento do sistema. En conclusión, este método consegue detectar un estremo global cunha precisión boa, aínda que non se pode asegurar ó cen por cen que o extremo sexa o óptimo. Canto máis lento sexa o decrecemento da temperatura maior seguridade se poderá ter no resultado. Por iso, a pesar de ser un sistema moi efectivo, é pouco eficiente computacionalmente. Algoritmos xenéticos A aplicación dos principios da xenética á busca adaptativa e a optimización de funcións deu lugar ó desenvolvemento do que se coñece como algoritmos xenéticos (AX). Os fundamentos sobre os que se constrúe un algoritmo xenético, a saber, os mecanismos de evolución dos cromosomas e o principio de supervivencia do mellor adaptado de Darwin, son simulados artificialmente sobre unha poboación de posibles solucións a un problema determinado. Esto da lugar a un eficaz sistema para aborda-la busca de extremos globais, válido para todo tipo de funcións, por complexas que sexan estas. Algúns estudios interesantes sobre este tipo de sistemas recóllense no traballos de Goldberg [25] e Holland [26]. Os AX son, coma o temperado simulado, técnicas probabilísticas de exploración do espacio de solucións dun problema, pero a diferencia doutras técnicas desta natureza, non traballan cunha única solución en cada iteración, que de aquí en diante imos chamar xeración, senón que se xera toda unha poboación de individuos, representados mediante unha cadea de parámetros a modo de xenes dun cromosoma. As operacións que conducen dunha xeración á seguinte son principalmente a reproducción, que é a transcrición sen alteracións dun cromosoma, o intercambio de xenes, chamado cruzamento (crossover), e a mutación. O papel que xoga a función a optimizar na evolución da poboación, é o de outorgar a un individuo maior ou menor probabilidade de intervir na creación da nova xeración, avaliando a bondade de cada solución, que se
3.2 Métodos de reconstrucción 43 correspondería co grao de adaptación ó medio. Pos eso á función a optimizar se lle chama función de adaptación (fitness function). Polo xeral, para a implementación das operacións de variación xenética, soe ser de gran utilidade a codificación dos vectores de parámetros, aínda que esto non é imprescindible. Habitualmente empréganse codificación de tipo binario. Así, unha representación moi usual dunha solución da función é unha cadea de bits formada pola concatenación das codificacións asociadas a cada un dos parámetros. Esto facilita tanto a operación de mutación, que consistiría na selección aleatoria dun dos bits da cadea e a negación do seu valor, como o cruzamento, para o que bastaría seleccionar unha ou máis posicións de truncamento da cadea e o intercambio dos anacos resultantes. Aínda que existen múltiples variantes do AX, a liña xeral de funcionamento é a que describe o seguinte pseudocódigo. codificación. inicialización aleatoria da poboación a N individuos. avaliación da poboación. repetir repetir seleccionar 1 individuo e reproducilo. seleccionar 2 individuos e recombinalos. seleccionar 1 individuo e mutalo. ata xerar N novos individuos. avalia-la nova poboación. ata que se acade o criterio de terminación. A terminación da execución do algoritmo pode fixarse cando se acada un determinado grao de aproximación da solución, pero para iso hai que dispor dun certo coñecemento sobre a función que se está a optimizar. Outra solución é deter o proceso de evolución cando se deixen de observar variacións significativas na composición da poboación. Por se esta situación tampouco chega a darse, tamén se pode limitar o número de xeracións. O por que do funcionamento dos AX está implícito na forma de representar e manipular as variables. A hipótese dos bloques constituíntes (building blocks) xustifica a converxencia do algoritmo baseándose no estudio dos esquemas (schemata) que compoñen a cadea de xenes. En [5, 25, 27] preséntase unha introducción a esta teoría. Un esquema é unha secuencia de valores determinados dun subconxunto de xenes pertencentes ó cromosoma –na codificación binaria soen representarse polos valores correspondentes, marcando o resto dos bits como variables, como por exemplo 10***, ou 11*0*-. A lonxitude dun esquema mídese como o número de bits de separación entre o de maior e o de menor orde. Os esquemas representan rexións no espacio de solucións. Para que o AX converxa, os seus mecanismos de evolución deben conseguir que dentro da poboación abonden os individuos con esquemas correspondentes a
44 CAPÍTULO 3. MODELOS DE SUPERFICIE rexións próximas á solución óptima, ou polo menos a solucións boas. Mediante o principio de supervivencia do máis adaptado permítese unha maior participación dos individuos con esquemas bos nas operacións de cruzamento, permitindo que esquemas curtos se combinen para formar esquemas longos, é dicir máis preto da solución óptima. Son os esquemas curtos os que xogan un papel máis importante, pois son os que maior probabilidade de sobrevivir xeración tras xeración, dado que canto máis longo é o esquema maior é a probabilidade de ser fragmentado. Por este motivo se tentará deseñar algoritmos que permitan aproveitar ó máximo a información aportada por estes esquemas curtos asociados a altos valores da función de adaptación, que se soen chamar bloques constituíntes. Pola súa banda, a mutación serve para introducir esquemas novos na poboación. Deste xeito, cruzamento e mutación constitúen un método de exploración excelente na procura do extremo global da función. Tódalas consideracións anteriores marcan unha serie de diferencias co resto de técnicas de optimización de propósito xeral, que a continuación se resumen. O AX traballa con varias solucións simultaneamente. Digamos que ten un paralelismo implícito. Ademais, pódese dicir que as solucións existentes nunha determinada xeración teñen memoria das xeracións anteriores, a diferencia de outras técnicas como o SA, onde a información de estados anteriores non se aproveita. Estas dúas características permiten unha visión máis global do espacio de busca e polo tanto maior probabilidade de dar co resultado correcto. O AX non traballa coas variables do problema, senón que polo miúdo o fará unha cunha codificación das mesmas. É dicir, as operacións de evolución xenética non precisan dispor de ningún tipo de coñecemento sobre o problema en particular, nin sobre a natureza das funcións a optimizar. Non sucede así cos algoritmos de gradiente, que a parte de precisar dunha inicialización próxima á solución desexada, necesitan realizar cálculos do gradiente e da matriz hessiana da función. Se a función é demasiado complexa, estes cálculos serán aproximados e introducirán erros. O tratamento das restriccións nos valores dos parámetros é moi sinxelo. Só é necesario engadir un termo de penalización á función cando se estea a avaliar unha solución non valida. Como inconveniente, hai que dicir que os AX aseguran a converxencia ó estremo global pero nun tempo ilimitado, o que quere dicir que nos podemos cansar de esperar e que o criterio de terminación se cumpra antes de que se acade unha boa estimación da solución óptima. Sen embargo, este problema pódese paliar en gran medida realizando unha implementación paralela, o que resulta sinxelo dadas as características do algoritmo.
45 4 SUPERCUÁDRICAS As funcións implícitas chamadas supercuádricas ou superelipsoides son unha familia de superficies paramétricas que se caracterizan pola flexibilidade que posúen en contraposición coa simplicidade da súa representación. Moitos son os autores que nos últimos anos se están interesando polo uso desta representación, ben sexa de modo único ou combinada con outros modelos. Proba desto é a multitude de traballos publicados sobre aplicacións que fan uso de superelipsoides, dos que cabe destacar os de Terzopoulos e Metaxas [39, 55], Solina [30, 51], Bardinet [2, 3, 4], e Chen [15] entre outros . Os superelipsoides son a extensión a 3D das figuras 2D coñecidas como superelipses, ou supercuádricas. Estas obtéñense a súa vez como unha extensión das cuádricas, substituíndo o expoñente 2 das ecuacións dunha curva cónica por unha variable real. Nunha parametrización natural, a superelipse ten a seguinte expresión θ πθπθ ε ε sen - ,cos by ax = <≤= onde o parámetro ε é un número real positivo e a e b son os radios sobre os eixes x e y respectivamente. Como as funcións trigonométricas coseno e seno poden tomar valores negativos, hai que substituí-la operación de potenciación, de forma que en todo o rango de definición de θ estea tamén definida a superelipse no dominio dos números reais. Entón, na expresión da superelipse substitúese a operación de potenciación dunha variable u a un número real positivo ε pola seguinte ( ) ε ε uusignou =: O parámetro ε coñécese como parámetro de cuadratura, pois é o responsable da forma da curva. Cando ε tende a cero, a superelipse aproxímase a unha figura rectangular, sendo a e b as dimensións dos lados, para ε =1 tense a ecuación dunha elipse e para ε =2, a superelipse é un rombo con a e b como semidiagonais. Para valores maiores os lados dos rombo vanse curvando cara ó interior, formando unha figura cóncava con aparencia de estrela. Realizando o producto esférico da superelipse consigo mesma, obtéñense as ecuacións do superelipsoide. ( ) ( ) () () () θ πφπφθ πθπφθ ε εε εε 1 3 21 2 21 1 sen - ,sencos 2/2/- ,coscos az ay ax = <≤= <≤=
46 CAPÍTULO 4. SUPERCUÁDRICAS A superficie que describen estas ecuacións pode tomar formas moi diversas en función dos valores concretos dos seus parámetros. O parámetro ε1 determina a cuadratura nos planos que conteñen ó eixo z e ε2 nos planos normais ó eixo z. Se z, ε1 e ε2 son moi pequenos a superficie aproximarase a un paralelepípedo, para ε1<<1 e ε2=2 temos un cilindro e con ε1=ε2=2 as ecuacións anteriores corresponden ás dun elipsoide. Para valores destes parámetros maiores que 2 a superficie presenta concavidades. En calquera caso, as supercuádricas sempre son simétricas respecto ós eixes principais. A ecuación implícita correspondente é a seguinte: () 1, 1/2 3 1/2 2/2 2 2/2 1 = + + = ε εε εε a z a y a x qrf rr sendo () zyxr ,,= r o vector de posición dun punto sobre a superficie e () 21321 ,,,, ε ε aaaq = r o vector de parámetros do superelipsoide. 4.1 Medidas de erro Para poder empregar os superelipsoides como esquema de representación do modelo de forma global, é preciso definir unha función de medida do axuste da superficie ós puntos da imaxe. Como se comentou anteriormente, os dous métodos máis usados para a optimización de funcións implícitas son o que define unha función que estima a probabilidade de que haxa correspondencia e a que emprega funcións relacionadas coa distancia entre modelo e imaxe. Esta última é a opción escolleita no presente estudio para a optimización dos superelipsoides. Se se dispón dun conxunto de N puntos da imaxe pertencentes á fronteira do obxecto, a función E que mide o erro cometido no axuste dos parámetros da superficie, será a suma das contribucións D de cada un dos puntos. () ( ) ∑ = =N i iqrDqE 1 2,rrr O cálculo exacto da distancia dun punto a unha superficie pode resultar moi custoso computacionalmente, polo que se soe utilizar unha aproximación para o cálculo de D. Para o caso concreto das supercuádricas a propia ecuación da superficie pode ofrecer unha medida de erro, posto que a función ( ) qrf r r ,, para valores positivos de qi,, é monotonamente crecente coa distancia radial se o punto pertence á rexión do espacio que queda fóra da superficie, e monotonamente decrecente no interior. • () qrf rr,=1 o punto está sobre a superficie • () qrf rr,>1 o punto está no exterior da superficie • () qrf rr,<1 o punto está no interior da superficie Por esta razón ( ) qrf r r , é coñecida como a función interior-exterior. Unha análise en profundidade desta función e de como se pode utilizar para obter medidas da función de
4.1 Medidas de erro 47 erro pódese ver mo traballo de Whaite e Ferrie [58]. A función de erro que se deriva directamente desta propiedade é a seguinte: ( ) 1, 1 − = qrfD r r Sen embargo esta non é unha boa función de erro, pois asigna valores menores a superelipsoides cun valor de ε1 maior. Por esta razón, as superficies resultantes dun axuste con esta función de erro tenderán a ter unha sección transversal máis redondeada que a superficie óptima. Unha función moi empregada que reduce este problema é D2, defínese a partir da anterior, de xeito que compensa este comportamento. ( ) 1, 1 2−= qxfD r r ε En ocasións vai ser necesario realiza-lo axuste sobre imaxes incompletas, nas que non se dispón de información sobre a forma dunha parte da superficie. Pode darse o caso de que, por exemplo, nunha zona da imaxe haxa unha perda de contraste, o que provoca que non se detecte o borde do obxecto. Noutros casos é preciso facer unha reconstrucción de obxectos 3D a partir dunha única vista 2D, polo que hai unha oclusión da parte posterior do obxecto. En ambas situacións, a utilización de superelipsoides, soluciona este problema de falta de información baixo a hipótese de simetría axial. A pesar desto, a ambigüidade non queda totalmente resolta, xa que para unha imaxe descrita parcialmente pode haber varias superficies igualmente válidas. A medida de erro D2, sen embargo, vai avaliar estas superficies de distinto modo, asignando valores menores ás superficies que encerran un volume maior. Esto débese a que puntos con igual distancia radial á superficie contribúen máis se están no exterior que se están no interior, polo que o tamaño tenderá a aumentar para aproximarse ós puntos do exterior e compensar esta diferencia. Dado que o habitual é realizar hipóteses de volume mínimo é frecuente encontrar que a medida de erro D2 se modifica, engadindo un factor g(a1,a2,a3) relacionado coas dimensións do superelipsoide, que penaliza ás superficies que encerran un volume maior. No artigo de Solina e Bajcsy [51] realízase un estudio empregando distintas expresións para o factor g, para chegar a conclusión de que a función de erro que perceptualmente mellores resultados da é ( ) ( ) { } 1, 1 2/1 3213 −= qxfaaaD r r ε Ningunha das funcións presentadas ata agora é invariante ante transformacións de escala. Para que se cumprise esta condición necesitariamos unha medida de erro que fose unha métrica euclídea, pero, como se dixo antes, o cálculo da distancia exacta sería moi custoso, debido á gran non linearidade da función. Unha aproximación máis sinxela de calcular, que si é euclídea, é a distancia radial δ do punto á superficie (ver figura 4.1), estudiada por Whaite e Ferrie [58] e por Bardinet e Col. [2]. Se trazamos unha recta que una o punto en cuestión co centro do superelipsoide, a distancia radial mídese entre os puntos r e r0, no que a recta interseca á superficie. Esta distancia pódese calcular facilmente coñecendo a relación que garda coa función f. Da distancia radial sabemos que verifica 0 rr r r += δ . Como non coñecemos a expresión de 0 r r, buscaremos unha segunda ecuación que a relacione cos termos
48 CAPÍTULO 4. SUPERCUÁDRICAS coñecidos. Na figura 4.1 pódese ver que, dado que os puntos o, r e r0 son colineares, débese cumprir que rr rr µ = 0, sendo µ >0. É fácil comprobar que a función () () rfrF rr 2/1 ε = é homoxénea de orde 1, polo que se ten que F(µx, µy, µz)= µF(x,y,z). Tendo en conta que F(0 r r)=1, obtemos µ=1/F(x,y,z). Operando con estas expresións tense 1/2 0 1 ε δ += r fr Desta expresión derívase unha ecuación para δ, da que se pode eliminar 0 r r, tendo en conta que 0 rr rr += δ () −== qrf rD rr r , 1 11 4 ε δ Hai que recalcar que esta última función de erro sobreestima a medida exacta da distancia, excepto naqueles puntos onde a dirección da recta que une o punto co centro do superelipsoide coincida coa normal á superficie en 0 r r . Con todo, os autores están de acordo en que esta última é a mellor estimación do erro do axuste, a pesar de que en casos especiais poida ser conveniente empregar algunha das outras. Por último, é posible facer unha aproximación de primeira orde á distancia euclídea exacta, d, como apunta Bardinet en [2], que corresponde a distancia medida sobre unha recta normal a superficie. ( ) rFrFD r r ∇−= /1)( 5 Esta aproximación é boa tan só para puntos próximos á súa proxección na superficie, ou sexa, para valores pequenos de d. Ademais require o cálculo do gradiente de F. 4.2 Xeneralización para posición e orientación arbitrarias Calquera que sexa a medida de erro que se escolla, vai ser sempre preciso o cálculo da función interior-exterior () qrf rr,. Tal e como a presentamos, esta función describe unha figura centrada na orixe do sistema de referencia e cos eixes principais ó longo dos eixes coordenados, mentres que o obxecto a modelar pode estar situado en calquera parte dentro do volume da imaxe e cunha orientación calquera. Se unha superficie descrita por unha ecuación implícita non ten o centro na orixe de coordenadas, e os eixes principais orientados na dirección dos eixes coordenados, na ecuación da superficie van aparecer novos parámetros, que ata agora non se tiveron en conta. A expresión xeral dos puntos da superficie obtense aplicando as transformacións de translación e rotación necesarias para situala primitiva nunha posición e orientación xenéricas.
4.3 Ambigüidade na representación mediante superelipsoides 49 Para determina-la transformación T que hai que aplicar ó superelipsoide para que se superpoña co obxecto na imaxe, ou sexa, para pasar o superelipsoide do sistema de referencia da imaxe a outro centrado no obxecto, será necesario introducir os parámetros correspondentes no proceso de optimización, polo que, ós cinco parámetros de forma do superelipsoide, haberá que engadir os tres de translación, t1, t2, e t3, máis os tres ángulos de Euler, α, β e γ que describen a rotación. () γ β α ε ε ,,,,,,,,,,' 32132121 tttaaaq = r Transformar o superelipsoide respecto á imaxe para superpoñe-la co obxecto, equivale á realiza-la transformación inversa sobre a imaxe para superpoñe-lo obxecto co superelipsoide. Se r’ é o punto da imaxe dado no sistema de referencia centrado no obxecto, entón hai que substituír ( ) qrf r r , na función de erro seleccionada por esta outra expresión () ( ) ( ) ( ) ( ) qrfqrTfqrfqrf r r r r r r r r ,',',, ˆ1=== − Deste xeito pódense incorporar novas transformacións á definición da superficie, como veremos máis adiante, cando apliquemos deformacións globais ó superelipsoide. 4.3 Ambigüidade na representación mediante superelipsoides Na sección dedicada ós distintos esquemas de representación adiantabamos xa que os superelipsoides xeneralizados non cumpren a propiedade de unicidade. As simetrías de rotación que caracterizan tanto a cuádricas como a supercuádricas, que provocan unha ambigüidade na determinación dos ángulos de Euler. O caso máis claro é o da esfera, para a cal os ángulos de Euler están totalmente indeterminados. Coa figura do cubo, calquera rotación de 90º respecto ós eixes principais da lugar a un novo valor válido do ángulo de Euler. En xeral, calquera outra forma do superelipsoide presentará polo menos simetría ante calquera rotación de 180º en torno ós eixes principais. Pero as simetrías non son as únicas fontes de ambigüidade da representación. Tamén é posible obter formas idénticas variando os valores dos parámetros de cuadratura. A orixe deste comportamento está en que os valores ε<<1 e ε=2 na ecuación da superelipse, da que se deriva o superelipsoide, correspóndense respectivamente co Figura 4.1 Distancias radial, δ, e normal, d, dun punto a unha superficie.
56 CAPÍTULO 4. SUPERCUÁDRICAS ( ) ϕ θ γ β α ε ε ,,,,,,,,,,,,,,, 32121321 kCCtttaaaq yx = r Para unha imaxe dada, a obtención do modelo particular farase empregando algunha das funcións de erro descritas nesta sección para o axuste ós puntos de fronteira. O método de axuste da función elixida que se vai utilizar nesta ocasión é un algoritmo xenético que, como se demostra en [31], resulta moi axeitado para este tipo de funcións. No capítulo 3 fíxose unha introducción ós AX, da que se concluíu que é un algoritmo de optimización global moi robusto, pero que o seu tempo de converxencia pode ser excesivo. Como a aplicación a desenvolver non está pensada en principio para sistemas dinámicos, a rapidez non é aquí unha prioridade. Sen embargo, existe un serie de variacións sobre o AX básico, que de seguido se van a especificar, que permiten axilizar o proceso sen afectar substancialmente os seus principios de funcionamento, mantendo tódalas prestacións do AX máis sinxelo. Utilizando algoritmos xenéticos temos unha gran liberdade á hora de especificar cales e como van se-los mecanismos evolutivos. Xa só do operador de recombinación pódense implementar varias versións, que serán máis ou menos axeitadas segundo o tipo de problema a tratar. Ademais, non existe unha norma que nos indique como fixar certos parámetros do algoritmo, como o número de individuos da poboación, probabilidades de cruzamento e mutación, porcentaxe de individuos que se reproducen, etc., polo que a elección destas magnitudes se fará de forma empírica. O esquema de creación dunha nova xeración de individuos inclúe a aplicación de operadores de reproducción, mutación e cruzamento. No traballo que aquí se presenta, experiméntase con tres variantes do cruzamento, que se diferencian no número de puntos polos que se secciona un par de cromosomas para a recombinación: por un punto, por dous puntos e uniforme, onde cada xene do cromosoma se pode intercambiar independentemente dos demais. Este último operador ten o inconveniente de que non favorece a formación de bloques constituíntes, dada a elevada fragmentación dos cromosomas, o que dificulta a herdanza da información contida nos cromosomas dunhas xeracións a outras. 4.6.1 Codificación Decididos os operadores a empregar, o primeiro e posiblemente máis importante aspecto a considerar é como codificar o vector de parámetros do modelo. Para dar coa codificación máis axeitada haberá que ter en conta o seguintes factores: • como afecta o tipo de codificación, é dicir, que sexa binaria, real, etc., á eficiencia dos operadores xenéticos. Esta eficiencia estará relacionada coa capacidade para xerar esquemas útiles; • como inflúe na eficiencia a ordenación dos parámetros na cadea; Posto que os algoritmos xenéticos constitúen un método de optimización moi robusto, o tipo de codificación que se manexe, non afectará á efectividade do algoritmo, pero si a eficiencia. Por exemplo, coas codificacións real e binaria os operadores xenéticos impleméntanse de forma distinta. Se unha variable está representada por un número
4.6 Optimización con algoritmos xenéticos 57 real, ese número constituirá por si só un xene do cromosoma, mestres que na codificación binaria, cada un dos bits equivale a un xene. Esto fai que a implementación dos operadores xenéticos para a codificación binaria sexa máis sinxelo, posto que a representación como unha restra de bits permite que se modifique un único bit ou que se seccione en dous anacos o parámetro. Esto non é posible coa codificación real, para o que hai que definir operadores máis complexos, como os que propón Kalyanmoy en [27]. Unha opción mellor incluso que a codificación binaria é a utilización do código Gray [25, 27]. A vantaxe respecto ó código binario normal radica en que dúas cadeas de bits que representan números consecutivos diferéncianse en un e só un bit. Desta forma, a modificación dun bit da cadea que representa o cromosoma, ten sempre a mesma repercusión, independentemente de en que posición da cadea estea situado dito bit nin, dos valores que nese momento teñan o resto de bits. Non sucede así co binario normal, onde para pasar do valor 7 ó 8, hai que modificar os 3 primeiros bits e sen embargo pasar de 1000 a 1100 implica variar o valor da cadea en 8 unidades decimais. Esto provoca que haxa interacción entre os xenes. A interacción, en ocasións chamada epistase, entre os xenes consiste en que a contribución dun xene á función de adaptación depende do valor de outros xenes. No que respecta á ordenación dos parámetros na cadea, é importante que se faga de modo que as variables que están relacionadas entre si se sitúen en localizacións contiguas, pois esta é a maneira de que os esquemas correspondentes a estes xenes sobrevivan ó cruzamento. Se se atopasen moi separados na cadea correrían maior risco de seren separados. O problema é que moitas veces estas relacións entre os parámetros non son coñecidas. Outras veces existen demasiadas interrelacións entre as variables como para poder xuntalas a todas adecuadamente. Entón o obxectivo será buscar a colocación óptima. 4.6.2 Rendemento do AX. Outro aspecto a ter en conta no deseño do algoritmo xenético é a elección de parámetros coma o tamaño da poboación, que é o número de individuos por xeración, ou a probabilidade de que, unha vez seleccionados os individuos que van intervir nunha operación de mutación ou cruzamento, esta sexa efectiva. Na escolla destas cantidades haberá que estudiar o seu efecto tanto na eficacia como no rendemento do algoritmo, que, coma sempre, serán obxectivos en conflicto. No caso do tamaño da poboación sucede que, canto maior é esta, maior probabilidade existirá de dar coa solución correcta, posto que haberá un maior número de esquemas en xogo, o que se traduce nunha exploración máis ampla do espacio de busca. Pero, por outra parte, a poboación tardará máis en converxer cara a rexión de función de adaptación alta, pois haberá que ir eliminando un gran número de esquemas non adaptados. Coa probabilidade de mutación sucede algo parecido. Se é moi baixa, realizarase unha exploración precaria do espacio. Lembremos que a mutación encárgase de xerar novos esquemas dentro da poboación, e se esto non se produce en suficiente medida, a poboación estancarase, non haberá evolución. Pola contra, se a mutación se produce con demasiada frecuencia, pode suceder que esquemas moi adaptados desaparezan temporalmente ou tarden en dominar na poboación, o ralentizaría o proceso.
58 CAPÍTULO 4. SUPERCUÁDRICAS No estudio realizado por De Jong [16], no que se testaron diversas configuracións do AX na optimización dun conxunto de funcións de características diversas, puido comprobarse que os mellores resultados se obteñen para tamaños da poboación moderados, probabilidade de cruzamento elevada e probabilidade de mutación baixa, da orde da inversa do número de individuos da poboación, o que correspondería a unha media de unha mutación por xeración. Mellores estimacións destes parámetros se poden establecer en base a cálculos estatísticos para problemas concretos, como fixo Kalyanmoy [27] para o cálculo da poboación inicial. 4.6.3 Mecanismos de selección A converxencia do AX está supeditada á correcta simulación do principio de supervivencia dos máis adaptados. A maneira máis común de simular este comportamento é a de permitir unha maior intervención dos individuos con maior función de adaptación nos procesos responsables da creación das sucesivas xeracións. O modo máis simple se levar esto a cabo é asignar a cada cromosoma unha probabilidade de ser seleccionado como pai dun novo cromosoma que sexa directamente proporcional ó seu grao de adaptación, no caso do problema da maximización da función de erro –un problema de minimización pode transformarse facilmente nun de maximización substituíndo a función pola súa inversa -. Este o que se chama esquema de selección de ruleta (roulette wheel) ou tamén mostraxe estocástica con reposición. Pero os AX que utilizan este método adoecen dun descenso prematuro da variedade de esquemas na poboación. As causas deste comportamento son, por un lado, o feito de que non se pode medir o valor promedio da función de adaptación para un esquema determinado, e por outro, que o número de veces que un individuo intervén nos procesos evolutivos presenta unha gran desviación respecto ó número que se estima que é o correcto, debido a que en cada selección non se ten en conta cales foron as seleccións anteriores. Existen múltiples esquemas de selección que tentan resolver estes problemas. No libro de Goldberg [25] realízase unha comparación de varios destes esquemas. Un método de selección que tenta paliar este último problema é o modelo de valor esperado, tamén chamado mostraxe estocástica sen reposición, para diferencialo do esquema de ruleta, pois agora, a medida que un individuo vai sendo seleccionado para o cruzamento, vai perdendo posibilidades de volver a selo. Esto lévase a cabo calculando para cada individuo o número de descendentes que debe ter como a súa probabilidade de selección pi=fi/fmedia, e cada vez que sae elixido para cruzarse, réstaselle 0.5 se o cruzamento é efectivo e 1.0 se non o é. Este método reduce de considerablemente a perda precoz de esquemas. Outro método de selección moi popular é o torneo (tournament selection). Consiste en tomar toda a poboación e dividila en grupos de n individuos. De cada grupo escóllese o máis adaptado para participar nos procesos de recombinación. Habitualmente os individuos compiten por parellas. Deste xeito asegúrase que o máis adaptado da poboación participe na formación da seguinte xeración. Unha variante de esta técnica, chamada selección de torneo probabilístico (probabilistic tournament selection), permite que nas primeiras xeracións exista unha probabilidade de que o menos
4.7 Posta en correspondencia con modelos a priori 59 adaptado dos dous individuos saia elixido. Esta probabilidade vai descendendo a medida que evoluciona o sistema, de forma similar a como o fai o temperado simulado. 4.6.4 Función de adaptación A probabilidade que se lle asigna a un individuo de formar parte na creación dunha nova xeración, vai depender unicamente da súa función de adaptación. Se empregamos unha técnica de selección proporcional ó valor da función de adaptación e na poboación inicial aparecen individuos cunha adaptación moi alta, aínda que non estean cerca do óptimo, de contado dominarán a poboación, contribuíndo de novo á perda de esquemas e polo tanto reducindo o grao de exploración do espacio de busca. Por esta razón é práctica habitual realizar un transformación das mediadas de erro. A forma máis sinxela é realizar unha ordenación por adaptación dos individuos, chamada ranking. Sendo f’ a función de adaptación transformada, f’max e f’min as cotas do intervalo ó que se mapea f, e rank(i) a función que da a posición do individuo i nunha ordenación descendente con f, entón: () min minmax 'rank 1 '' 'fi N ff f− − − = Tamén se soe facer unha normalización linear ( ) ( ) iNf rank' − ⋅ = σ en función da desviación estándar σ dos valores de f, sendo N o número de individuos da poboación. 4.7 Posta en correspondencia con modelos a priori O procedemento que se acaba de describir para a obtención dos parámetros da ecuación implícita é igualmente válido para a reconstrucción dun modelo e para a posta en correspondencia entre o obxecto nunha imaxe e un modelo a priori. A única diferencia está en que no segundo caso o número de parámetros a determinar é menor, pois toda a información asociada á forma da superficie pode ser “reciclada”. Entón a optimización vese reducida ó cálculo dos parámetros de transformación ríxida: translación, rotación e escalado e un novo parámetro de escalado, a0, que se introduce nas ecuacións do superelipsoide do seguinte xeito. ( ) ( ) () () () θ πφπφθ πθπφθ ε εε εε 1 30 21 20 21 10 sen - ,sencos 2/2/- ,coscos aaz aay aax = <≤= <≤= O obxectivo que nos propuxemos neste traballo, a inicialización de superficies deformables, pode abordarse mediante a segunda das dúas tarefas. Esto require da definición do modelo a priori, pero a reconstrucción deste modelo farase só unha vez. Para que o modelo sexa representativo da clase de obxectos cos que vai traballar, é
60 CAPÍTULO 4. SUPERCUÁDRICAS preciso utilizar algún tipo de técnicas estatísticas para a súa obtención. É usual utilizar para a optimización da superficie a priori unha imaxe promedio do obxecto, calculada a partir dun conxunto de imaxes de proba. Tamén é interesante a utilización dun modelo PCA, o que requiriría, ademais da optimización dos parámetros de transformación ríxida, os pesos da combinación linear de modos de variación principais. Esta técnica pode ser especialmente útil cando se considere que o modelo promedio non é suficientemente representativo do conxunto de imaxes de proba, é dicir, cando os diversos exemplares do obxecto presenten unha gran varianza. Outra opción é empregar como modelo a priori o modelo medio, pero permitindo na posta en correspondencia a variación dos parámetros de forma, ou sexa, transformacións globais non ríxidas, aínda que aquí se faría de forma moi limitada. En ambos casos, para a extracción da colección de puntos sobre os que se axustará o modelo a priori, é preciso poñer en correspondencia as imaxes do conxunto de proba. Calquera dos métodos comentados na sección 3.2.1 para posta en correspondencia entre distribucións espaciais de puntos pode ser utilizada nesta ocasión.
61 5 EXTRACCIÓN DE PARCHES DE SUPERFICIE Ata agora fíxose unha descrición bastante detallada da reconstrucción dos superelipsoides a partir dos puntos da imaxe, pero en ningún momento se fixo mención a como se efectúa a selección de puntos que interveñen na optimización. Neste capítulo explícanse os procesos de baixo nivel empregados, os cales ofrecerán como resultado unha distribución de puntos de fronteira, que serán os que se empreguen para o cálculo das distancias á superficie do modelo. Unha das solucións máis comunmente usadas para a detección de puntos de fronteira en 2D é o filtrado espacial para a estimación das derivadas direccionais dos niveis de gris das imaxes. Nesta ocasión débese aborda-la tarefa do filtrado en imaxes volumétricas, e entón debemos preguntarnos ata que punto a aproximación 2D se pode estender ó tipo de imaxes 3D que imos manexar garantindo a calidade dos resultados. O problema que aparece é que as imaxes médicas, como por exemplo as provenientes de CT, non son verdadeiramente 3D, senón un conxunto apiado de imaxes 2D. A repercusión que este feito ten sobre a calidade dos datos é evidente. Todo dependerá da resolución na dirección axial que se ofreza, pero xeralmente observarase que a frecuencia de mostraxe nesa dirección será inferior que nas direccións contidas no plano das imaxes 2D, dando lugar a descontinuidades nos valores de intensidade ó pasar dun corte ó seguinte. Cando a distancia entre cortes non é excesiva pódese efectuar unha interpolación que xere imaxes en planos intermedios, o que daría maior continuidade á secuencia. A pesar deste inconveniente, o método que se aplica para o filtrado 3D das imaxes neste traballo é unha extensión directa do filtrado 2D. Para solucionar o problema que acabamos de citar, considérase que un suavizado previo é suficiente. O motivo polo que, a pesar de todo, se elixe este tipo de detección, é que a extensión de filtros 1D a filtros 3D separables poden ser implementados de forma recursiva. Deste xeito, o tempo de computación vese reducido de forma importante en relación co correspondente á convolución cunha máscara no dominio espacial, e dado o volume de datos manexados en imaxes 3D, esta idea resulta moi atractiva. A continuación farase unha descrición dos filtros recursivos. Estes filtros son empregados para o suavizado e a estimación das primeiras e segundas derivadas direccionais de cada punto da imaxe. As derivadas primeiras utilízanse para o cálculo da posición dos máximos do gradiente, nos que se localizarán os puntos de fronteira. Unha vez localizados os puntos de fronteira, procederase a agrupalos en parcelas de superficie. Como van a existir múltiples estructuras na imaxe, e como ademais a superficie do obxecto de interese pode aparecer fragmentada, vai ser necesaria unha etapa de selección dos parches que intervirán no axuste do modelo. Aquí é onde intervén o cálculo das derivadas segundas direccionais, que van servir para a estimación das curvaturas extremas locais de cada parche de superficie. A partir destes
62 CAPÍTULO 5. EXTRACCIÓN DE PARCHES DE SUPERFICIE cálculos é posible extraer algunha medida que poida servir como descriptor da forma dos parches. Con estas e outras características que se poden obter para cada parche, establécese un criterio de selección. A extracción e selección de parches de superficie detállase na sección 5.4. 5.1 Filtros multiescala recursivos Nesta sección descríbese en detalle a implementación utilizada para o filtrado encargado de estimar as primeiras e segundas derivadas direccionais das imaxes. Estudiaranse dúas aproximacións diferentes, para as que se levará a cabo unha comparación das prestacións que ofrece cada unha no capítulo de resultados. Ambas aproximacións representan funcións de resposta de impulso dependentes de parámetros asociados á escala do filtro, de modo que poidan ser utilizados con imaxes de distinto tipo e distintas características de ruído. En ambos casos farase unha implementación recursiva da función co obxectivo de reducir o tempo de cálculo das imaxes filtradas. 5.1.1 Filtros multiescala As características que se lle esixen a un filtro de detección de puntos de fronteira dependen do uso que se lle queira dar. Un filtrado axeitado para a detección de estructuras nunha imaxe será con frecuencia pouco apropiado para a boa localización e detección de puntos de fronteira a nivel local. No primeiro caso precisase un filtro que poida captar información relevante na imaxe, mentres que no segundo farán falta operadores que preserven as pequenas estructuras sen alterar a súa forma. É dicir, que distintos niveis de abstracción requiren filtrados a diferentes escalas. Ademais, un obxecto real soe manifestar características que que se manifestan a escalas diferentes, polo que a escala ideal pode variar ó longo do contorno. En calquera caso, nas aplicacións máis usuais, como a que aquí se propón, a descrición mediante unha única escala é suficiente. Neste traballo aplicáronse, tanto para o suavizado das imaxes como para o cálculo das derivadas primeiras e segundas, dúas familias de filtros moi utilizadas en filtrado multiescala polas boas prestacións que ofrecen para este tipo de tarefas. Por unha parte, están o filtro de gaussiana e as súas derivadas [20]. O filtro de derivada de gaussiana, que para o caso unidimensional con desviación estándar σ vén dado pola expresión seguinte, permite unha alta supresión de ruído ademais da boa localización dos puntos de fronteira. Se n é unha variable discreta () 22 2/ 2 2 1 σ πσ n eng − = Por outra parte, utilizáronse outra familia de filtros que se derivan dun filtro de detección de puntos de fronteira deseñado utilizando o criterio de Canny [11] estendido a filtros infinitos, que ofrece unha expresión analítica para un filtro xenérico unidimensional definido no dominio continuo, f(x), dadas unha serie de especificacións en relación ós seguintes criterios:
5.1 Filtros multiescala recursivos 63 • boa detección, asegurando unha baixa probabilidade de detección dun falso punto de fronteira e ó mesmo tempo unha baixa probabilidade de fallo na detección dun punto de fronteira verdadeiro. Dado que ambas probabilidades son monotonamente decrecentes coa relación sinal-ruído, SNR (Signal to Noise Ratio), o criterio de boa detección satisfarase coa maximización desa magnitude. Se g(x) é un sinal unidimensional calquera cun punto de borde situado en x = 0, corrompida por ruído gaussiano branco de amplitude media por unidade de lonxitude n0, a SNR virá dada por: ( ) ( ) () () 2/1 2 0∫ ∫− =dxxfn dxxfxg SNR • boa localización, de xeito que a posición na que se sitúe a fronteira detectada do obxecto estea o máis próxima posible ó centro do borde, mediante a maximización dun factor de localización, Λ. ( ) ( ) () () 2/1 2 0' '' ∫ ∫− =Λ dxxfn dxxfxg onde o apóstrofe indica que se está tomando a derivada primeira da a función correspondente respecto da variable independente x. • resposta única, de maneira que a resposta do filtro ante un sinal ruidoso non teña demasiados máximos arredor do punto de borde. Para eso o que se fai é limita-la distancia media entre dous picos da resposta do filtro ante un sinal de ruído gaussiano, ou o que é o mesmo, entre dous puntos de cruce por cero da derivada segunda () ( ) () () 2/1 2 2/1 2 0'' ' ∫ ∫ =dxxf dxxf x π sendo f’’ a derivada segunda de f respecto de x. Aplicando estes criterios para a optimización dun filtro de derivada primeira f(x)=h’(x) de segunda orde, Deriche [18, 19] deseñou unha familia de filtros dende o filtro de suavizado, h(x) ata o de derivada terceira. As expresións correspondentes ó filtro de suavizado e ás derivadas primeira e segunda, nas súas versións discretizadas, son as que seguen. () () ( ) ( ) ( ) nnn enkkxhnekxhenkxh ααα αα −−− −==+= 2310 1'' ' 1 onde os parámetros k0, k1, k2 e k3 son factores de normalización. k0 vén dado pola condición
64 CAPÍTULO 5. EXTRACCIÓN DE PARCHES DE SUPERFICIE () ∑ ∞ −∞= = n nh 1 As expresións de k1, k2 e k3 escóllense de modo que h’ e h’’ teñan resposta nula ante unha entrada constante, esixindo ademais que () ∑ ∞ −∞= = n nh 0'' As expresións resultantes son: () ( ) 1 2 1 1 21 1 3 2 2 2 1 2 2 0= − = − −= −+ − =− − − − −− − k e e k e e k ee e k α α α α αα α αα Estes filtros veñen dados en función dun parámetro α relacionado coa escala á que opera o filtro, polo que a este filtro se lle chama filtro alfa. Se o filtro de suavizado de gaussiana e o alfa se deseñan de modo que teñan a mesma enerxía total, verifícase a relación α·σ=5/2√π. Valores típicos de α están entre 0.5 e 1.5. 5.1.2 Implementación recursiva dos filtros Tanto para a familia de filtros gaussianos [20], coma para os filtros alfa [17], realizouse unha implementación recursiva con obxecto de reduci-lo número de operacións por punto da imaxe. A expresión que relaciona o valor dun elemento de saída do filtro cos elementos de entrada, dado pola seguinte convolución, () ( )( ) ∑ − = −= 1 0 N k kixkhiy non é recursivo. N é o número de elementos da secuencia de resposta de impulso. Este tamaño de máscara do filtro pode resultar excesivo se se desexa obter boas prestacións en canto a localización e detección. Para a implementación, por exemplo, do filtro de gaussiana, débese tomar un valor de N que sexa alomenos da orde de 4σ, sendo σ a desviación estándar da gaussiana, para que ofreza bos resultados. A idea do filtrado recursivo é deseñar un filtro que teña un comportamento equivalente ó filtro de convolución, ou polo menos que se aproxime o suficiente, pero que veña descrito pola seguinte ecuación de diferencias () ( ) )( 1 1 0 kiyakixbiy n k k m k k−⋅−−⋅= ∑∑ = − = Con esta expresión redúcese o número de operacións necesarias para obter un valor da secuencia resultante de N a n+m, polo que o tempo de computación decrece de forma importante, xa que basta tomar valores moi pequenos de n e m para obter boas aproximacións. Haberá que determinar os valores de k a e k bda función de transferencia do sistema recursivo
5.1 Filtros multiescala recursivos 65 () () () ∑ ∑ ∑ = − − = −− ∞ = −− ⋅+ ⋅ == n k k k m k k k n n RR za zb znhzH 1 1 0 1 0 1 1 de maneira que o seu comportamento sexa igual ó do filtro non recursivo, ou, polo menos, se aproxime o suficiente de acordo con algunha medida de erro. Aquí aplícase un criterio de minimización do erro ó cadrado, onde o erro vén dado por: () ()() 2 0 2∑ ∞ = −= k Rkhkh ε sendo () khR a resposta de impulso do sistema recursivo. Filtros alfa No caso da aproximación recursiva dos filtros alfa, sucede que a función de transferencia asociada pódese implementar de forma exacta mediante a ecuación de diferencias antes presentada. As ecuacións de diferencias para un elemento de saída y(n) do filtro de suavizado son () ( )() ( ) [ ] ( ) ( ) () ( )( ) ( ) [] () ( ) () () () nynyny nyenyenxenxekny nyenyenxenxkny 21 2 2 2 2 12 1 2 101 212211 21211 += +−+++−++= −−−+−−+= −−−− −−− αααα ααα α α Desta forma é posible implementar un filtro de suavizado realizando 8 multiplicacións e 7 sumas, independentemente do valor de α, mentres que cun filtro non recursivo sería preciso un número de operacións proporcional ó tamaño do filtro. A modo de exemplo, un filtro de 16 bits de tamaño de máscara con α=0.5 necesitaría 57 operacións, mentres que con 0.25 serían 113. As ecuacións do filtro de derivada primeira, pola súa parte, precisan tan só de 5 operacións tanto de producto como de suma por elemento de saída. () ( ) ( ) ( ) ()() () ( ) () () () [] nynyekny nyenyenxny nyenyenxny 211 2 2 22 1 2 11 2121 2121 −= +−+++= −−−+−= − −− −− α αα αα As ecuacións de diferencias utilizadas para a estimación da derivada seguinte implican 7 operacións de suma e 7 de multiplicación por elemento de saída, e son as que seguen. () ( )( ) ( ) ( ) ( ) () ( )() ( ) () ( ) () () () nynyny nyenyenxenxkeny nyenyenxkenxny 21 2 2 2 2 2 1 2 11 212211 21211 += +−+++−+−= −−−+−+−= −−−− −−− αααα ααα α α Filtros de gaussiana
72 CAPÍTULO 5. EXTRACCIÓN DE PARCHES DE SUPERFICIE A función recursiva busca_puntos_adxacentes fai o seguinte: Entrada: último_punto, superficie Para cada punto adxacente a último_punto facer Se punto foi percorrido continua Marca punto como percorrido Se gris(punto) < gris_umbral_menor continua Engade punto á lista de puntos de superficie Aumenta a área da superficie nunha unidade último_punto ← punto Busca_puntos_adxacentes(último_punto, superficie) Fin Fin Unha vez rematado este proceso dispoñemos xa dunha organización da información de máis alto nivel que a aportada pola imaxe de niveis de gradiente de gris. Pero de momento as propiedades que se describen, tanto o modulo de gradiente de nivel de gris como as curvaturas, seguen sendo de carácter local, pois van asociadas a puntos individuais da imaxe, polo que non é posible facer unha discriminación de parcelas de superficie. Do que se trata agora é de operar sobre a nova descrición en forma de estructuras de superficie para determinar propiedades relativas, non a cada punto por separado, senón ó conxunto de todos eles. Desta forma seguimos subindo no nivel de abstracción, chegando a unha descrición que permitirá discriminar os puntos do obxecto a modelar. 5.4.2 Selección de parches de superficie As etiquetas que se asignan a cada parche de superficie, veñen dadas como unha función das súas características. Neste traballo estudiamos un número reducido delas pero o método é extensible a calquera tipo de etiquetado. As características que aquí se manexan son as seguintes: Área. Considérase que empregando filtros que operen a escala adecuada, as parcelas de superficie máis pequenas corresponden a ruído, e poden ser descartadas. É posible que con este criterio se estean eliminando anacos de pequeno tamaño da superficie do obxecto de interese, pero que non teñen peso no resultado final. Nivel de gradiente de gris promedio. Informa da magnitude da variación do nivel de gris. Esta variación pode ser baixa cando é moi gradual, pero tamén cando se produce un salto moi brusco entre rexións de niveis de gris moi próximos. De maneira que descartar parches con nivel de gris promedio baixo, podemos estar eliminando as superficies de interese. Esto, xunto co feito de que soe haber máis dunha estructura presente na imaxe, fan que a a medida do gradiente promedio non sexa suficiente para discriminar parches dun determinado obxecto. Curvaturas media e gaussiana. Aportan información sobre a forma local dunha superficie. Defínense a partir das curvaturas principais que na sección 5.3 foron calculadas. Sendo k+ e kos valores das curvaturas principais nun punto da superficie, a curvatura media H e a curvatura gaussiana K veñen dadas por:
5.4 Extracción de parches de superficie 73 −+ −+ ⋅= + =kkK kk H 2 Os valores destes descriptores en cada punto da superficie dan unha idea aproximada da morfoloxía local. A curvatura gaussiana utilízase para saber se as curvaturas principais teñen o mesmo signo. Aproximando a superficie localmente a unha función cuadrática pódese establecer a seguinte clasificación • K > 0, elipsoide • K < 0, hiperboloide • K = 0 e H ≠ 0, cilindro • H = K = 0, plano Pola súa parte a curvatura media informa sobre a convexidade da superficie. Asumindo que as rexións de nivel de gris alto corresponden ó interior dos obxectos e os niveis de gris baixos ó exterior, entón a curvatura local nunha dirección determinada será positiva cando a superficie sexa cóncava nesa dirección. Tendo esto en conta, o criterio para determinar a convexidade da superficie é o que segue • H > 0, superficie cóncava • H < 0, superficie convexa • H = 0 e K ≠ 0, punto de sela • H = K = 0, plano Tanto a curvatura media como a gaussiana permiten distinguir entre varias clases de superficies, indicando como é a súa curvatura, pero non canto. Para aportar unha clasificación cualitativa das superficies debemos introducir unha nova medida, que a continuación se define 22 −+ += kkKm Con estas cinco magnitudes que se acaba de definir xa podemos dispor dunha idea aproximada das propiedades da superficie en cada punto. Sen embargo, non son as propiedades locais as que interesan, pois quérese etiqueta-los parches de superficies, non os puntos de forma individual. As medidas que realmente se utilizan para clasificar unha superficie segundo a súa forma, estarán un promediadas a tódolos puntos de fronteira dunha parcela. A utilidade das medidas de curvatura promediadas estará limitada a parches de pequeno tamaño ou grandes superficies moi regulares. Cando, como resultado do filtrado da imaxe se obtén toda a superficie dos obxectos, se estes teñen unha forma moi complexa, ó promediar os descriptores de forma perderán o seu significado. Nestes casos basta unha clasificación por área e gradiente para a selección dos obxectos. En función de todas estas medidas pódese definir unha etiqueta que represente a clase de parche que temos. No traballo de Pardo e Cabello [44] suxírese unha expresión para o cálculo da etiqueta, aplicada o caso 2D, baseada nos valores do gradiente promedio e da lonxitude dos segmentos de borde. A expresión, que asigna valores maiores a
74 CAPÍTULO 5. EXTRACCIÓN DE PARCHES DE SUPERFICIE segmentos de curva S con valores do promedio do módulo do gradiente e da lonxitude l maiores, é a seguinte: () () 2/1 2 1lpI l SE i Sk ki ∇= ∑ ∈ sendo pk os puntos de contorno e I o gradiente. A expresión que se utilizou neste estudio é unha extensión da anterior a 3D, á que se lle engaden factores que teñen conta das medidas de curvatura anteriormente citadas. Dependendo do tipo de obxecto que haxa que segmentar haberá que incorporar á fórmula factores relacionados coa curvatura gaussiana e media dun ou doutro xeito. En canto a Km pediráselle usualmente que sexa o menor posible, pois as superficies orixinadas polo ruído soen ser moi rugosas. Se se dispón de coñecemento sobre os valores de Km esperados pódense introducir na expresión da etiqueta. A expresión xenérica da etiqueta dun parche P vén dada por: () () caso noutro , 2/12 max − = = CK HKf ag CKE PE m m onde a notación utilizada representa : () () () () () () () () ∑ ∑ ∑ ∑ ∈ ∈ ∈ ∈ = = = ∇= Pk kmm Pk k Pk k Pk k pK a PK pK a PK pH a pH pI a Pg 1 1 1 1 e sendo a a área do parche e C o valor esperado de Km. Se no se dispón de coñecemento a cerca de C, haberá que substituír o seu valor por cero. Segundo o tipo de superficie que se desexe segmentar defínese o factor f de distinta maneira. A continuación preséntanse algúns exemplos. Para superficies con caras planas, introdúcese un factor que aumente cando as curvaturas media como gaussiana tenden a cero en promedio, como pode se-lo seguinte ( ) ( ) HKHKf ⋅= /1, Se a superficie é cilíndrica haberá que distinguir entre o caso cóncavo e o convexo. Por exemplo, no caso cóncavo, f ten a seguinte expresión ( ) ( ) KHsignoHKf /, =
5.5 Observacións 75 Para o caso caso convexo, f cambia de signo. Cando se desexa segmentar superficies do tipo elipsoide ou hiperboloide hai que ter en conta o signo tanto de H como de K de forma separada. En concreto, a función utilizada para a segmentación da superficie externa da tibia a partir das imaxes de CT, deberá ter en conta que a súa forma é hiperbólica e cóncava, resultando unha expresión coma a seguinte. () caso noutro 1 0 e 0 ,0 0 2/12 max − ≠>< = = CK ag KKH CKE PE m m m Para segmentar a superficie interna esixirase que H sexa negativa. O problema con estas magnitudes é que son moi susceptibles de verse afectadas polo ruído. Se as curvaturas se calculan utilizando unha veciñanza de puntos de superficie moi pequena, a influencia do ruído será moi grande, e pode resultar que en lugar de estar medindo a curvatura da superficie esteamos medindo a rugosidade ocasionada polo ruído. Aplicando un filtro de suavizado atenúase este efecto. Tamén é conveniente utilizar filtros de escala grande para o cálculo das derivadas de nivel de gris. Por outro lado, é moi frecuente que aparezan parches duns poucos puntos en zonas de nivel de gris en principio uniformes debidos ó ruído. Estes parches poden ter curvaturas nulas, polo que o seu valor da etiqueta sería máximo. Para evitar esta situación recórrese a unha fase previa de umbralización por área, cun nivel umbral moi baixo. Unha vez calculada a etiqueta de cada superficie, a estratexia a seguir para a selección dos parches de superficie, pode ser unha simple umbralización do valor da etiqueta. A elección dun umbral sempre supón un problema. O habitual é determinar un umbral distinto para cada tipo de imaxes, segundo a técnica de obtención da imaxe e tamén por contido, sobre un conxunto de imaxes de proba. 5.5 Observacións A pesar da mellora que supón utilización de umbralización con histérese, que atenúa en certa medida o efecto de fragmentación das superficies, o proceso de extracción de parches de superficie non ofrece resultados totalmente satisfactorios. Ben é certo que o modelado global vai compensar a perdida de información producida no filtrado. Sen embargo, hai que ter en conta que o axuste vai ser mellor canto máis fiables sexan os datos extraídos da segmentación de baixo nivel. Se na umbralización se define o menor dos dous umbrais moi alto, as rexións da superficie de baixo contraste poden ser eliminadas. Dado que a forma do obxecto vén definida a priori, a desaparición dunha parte da superficie no filtrado non vai afectar ós parámetros da ecuación implícita nin das deformacións globais, pero si á localización, de xeito que o modelo tenderá a desprazarse cara a zonas onde a densidade de puntos de fronteira sexa maior. Se pola
76 CAPÍTULO 5. EXTRACCIÓN DE PARCHES DE SUPERFICIE contra o umbral menor é baixo de máis, a superficie do obxecto pode fusionarse con artefactos da imaxe, que tamén poden contribuír de forma importante ó erro do axuste. Pensando sobre todo na segmentación de imaxes de moi baixa calidade, pode ser interesante a utilización doutras técnicas de filtrado, que ofrezan mellores prestacións en canto a supresión de ruído e detección de estructuras. Referímonos concretamente ás técnicas baseadas en filtros de difusión non lineares e anisótropos. Trátase de suavizar os niveis de gris dunha imaxe de forma análoga a un proceso físico de difusión, coma o do calor. Establécese entón un fluxo que determina a evolución da configuración de niveis de gris co tempo. A imaxe filtrada corresponde ó estado en que se atopa o sistema despois difundirse durante un tempo t determinado, relacionado coa escala do filtro. O operador de suavizado obtense como resultado de resolver a ecuación de difusión () ()() ()() = ∇⋅= ∂ ∂ zyxIzyxu tzyxuDdiv t tzyxu ,,0,,, ,,, ,,, onde () zyxI ,, é o sinal de niveis de gris e D é unha matriz simétrica, chamada tensor de difusión. O filtro de suavizado gaussiano é un exemplo de filtrado que simula un proceso de difusión, pero neste caso séguese un modelo linear e isotrópico, é dicir, que a difusión se produce ó mesmo ritmo en tódolos puntos da imaxe e en tódalas direccións. Neste caso o tensor de difusión redúcese a unha constante, e a solución u á ecuación da difusión correspóndese co producto do sinal de niveis de gris I por unha función gaussiana que ten por desviación estándar σ=√2t. A non linearidade permite definir fluxos variables segundo a posición na imaxe, de modo que sexan menos intensos en rexións de fronteira. Esto da lugar a filtros que operan a escala alta en zonas uniformes, eliminando máis ruído, mentres que teñen un efecto máis local cerca dos contornos, preservando as estructuras presentes na imaxe. Outra vantaxe que ofrece a non linearidade é que atenúa o efecto de deslocalización dos contornos que ten lugar con filtrados coma o gaussiano. Un exemplo de filtro non linear isotrópico é o deseñado por Perona e Malik [46], onde o tensor de difusión é unha función escalar definida mediante a seguinte expresión 1 2 2 1 − ∇ += λ u D sendo λ>0 un factor constante. Introducindo este factor na ecuación da difusión, o fluxo tu ∂∂ / diminúe en puntos onde o gradiente aumenta. A maior desvantaxe deste filtro é que as zonas de fronteira se ven máis afectadas polo ruído. Pola súa parte a anisotropía dota ó filtrado da posibilidade de definir fluxos diferentes para distintas direccións. A difusividade en cada dirección virá determinada por factores estructurais. Un esquema que combine anisotropía e non linearidade permitirá establecer procesos de difusión lentos nas direccións normais en cada punto das
5.5 Observacións 77 superficies de fronteira, mentres que nas direccións tanxentes a estas será posible levar a cabo difusións máis rápidas. Por tanto, mantense unha boa localización e detección das zonas de transición entre rexións mentres que se ofrece un suavizado maior nas direccións tanxentes a esas superficies de fronteira, eliminando ruído. Como exemplo desta técnica téñense os filtros de Weicket [57]. Este último presenta unha aproximación na que se define o tensor de difusión en función de descriptores de forma relacionados con cálculos da curvatura, análogos ós presentados na sección anterior deste mesmo traballo. Esta é unha das vías de mellora desta metodoloxía que nestes momentos estamos considerando. O perfeccionamento das técnicas de preprocesado pode supoñer unha mellora importante na calidade dos resultados en imaxes moi ruidosas, co que a aplicación do modelado con superelipsoides pode estenderse a outro tipo de imaxes médicas de máis baixa calidade que as de CT, para as que este procesado demostrou ofrecer bos resultados.
79 6 RESULTADOS Acabamos de presentar unha metodoloxía orientada á resolución dalgúns problemas comúns dentro da reconstrucción de imaxes en xeral, e especialmente da imaxe médica, onde estes problemas se fan máis patentes. Definiuse un modelo capaz describir a estructura global dun obxecto, o que permite dotar dunha inicialización robusta ós modelos de deformación pola súa capacidade para evitar ou compensar a perda de información que se produce nos procesos de baixo nivel. O método deseñado benefíciase tanto do coñecemento dispoñible sobre a forma global, por medio do uso de modelos a priori, como da información a nivel local extraída da propia imaxe, a partir da cal se efectúa a identificación do obxecto na imaxe. As prestacións do sistema desenvolto dependen de tres aspectos: a capacidade dos modelos empregados para representar adecuadamente os obxectos no ámbito de aplicación do sistema, a eficiencia dos métodos de optimización, e a eficacia do preprocesado na extracción dos puntos de fronteira. Neste capítulo mostraremos algúns resultados que pretenden ilustra-las prestacións que nestes tres aspectos proporciona a técnica desenvolta. 6.1 Optimización da función implícita Comezaremos pola definición dos modelos a priori. Para a elección dunha das medidas de erro propostas na sección 4.1 realizáronse probas con imaxes de tibia. As función de erro que se estudian aquí son as D2, D3, e D4. ( ) 1, 1 2−= qxfD r r ε ( ) ( ) { } 1, 1 2/1 3213 −= qxfaaaD r r ε () −== qrf rD rr r , 1 11 4 ε δ A configuración do algoritmo xenético empregado para a optimización fíxose partindo dos valores dos parámetros da simulación aconsellados na obra de Goldberg [25], que posteriormente se someteron a diversas probas para un maior axuste. Na configuración final a función de erro sométese a un mapeado de ranking e emprégase selección de torneo probabilístico, que ofrece resultados considerablemente mellores que os de ruleta ou torneo sinxelo. A poboación está composta por un total de 80 individuos. O operador de recombinación é o de cruzamento por dous puntos, tomando a probabilidade de cruzamento como pc=0.8. En canto a mutación tomouse unha probabilidade bastante
80 CAPÍTULO 6. RESULTADOS elevada, pm=0.1, en vista da prematura estabilización da poboación obtida para valores máis baixos. Os parámetros codificáronse empregando código Gray de 16 bits para cada un. Aínda que sería posible aumentar o número de bits por parámetro, considérase que dado o tipo de aplicación do modelo non se precisa de máis resolución. Baixo estas condicións fíxose evolucionar a poboación ata observar que esta se estabilizaba ou se chegaba ás 2000 iteracións. Nunha primeira proba aplicouse o algoritmo así configurado a unha distribución de puntos de superficie dunha sección de tibia. As imaxes da figura 6.1 mostran distintas vistas desa superficie, representados mediante un esquema de triangulación para a súa mellor visualización. As imaxes das figuras 6.2 e 6.3 representan os resultados obtidos para as funcións D4 e D3 respectivamente. Visualmente pode apreciarse a diferencia de aspecto dos correspondentes modelos. Pero aínda que perceptualmente a función D4 ofrece mellores resultados, habería que realizar unha comparación cuantitativa para cercionarse de que efectivamente esta é a medida de erro máis axeitada. Para estimar a calidade dos distintos axustes utilizouse unha imaxe sintética dun superelipsoide, que se pode ver na figura 6.4. Esta imaxe presenta tódolos tipos de deformación global que se permite o noso modelo. Os valores concretos dos parámetros do superelipsoide pódense ver na táboa 6.1. A magnitude que nos serve para comparar numericamente os modelos obtidos defínese pola función () ',qqd , que representa a distancia entre o vector q dos parámetros coñecidos da superficie sintética e o vector 'q dos parámetros calculados utilizando unha determinada función de erro, dada polo módulo do vector diferencia q-'q. Para que o axuste dun parámetro en particular non teña máis peso no cálculo da medida de calidade que os demais, cada compoñente está escalada cun factor ai igual á orde de magnitude de qi, resultando a expresión de d como segue. () ( ) ∑−= i iii qqaqqd 22 '', Os valores dos factores ai aparecen tamén na táboa 6.1 ó lado dos correspondentes parámetros. En lugar de realizar unha soa proba con cada función de erro, leváronse a cabo varios axustes, para obter cadanseu vector 'q de parámetros calculados promedio. Os valores de d que se obtiveron son as que se contemplan na táboa 6.2. Pódese observar que en xeral, os resultados son peores para a función D2 que para as outras dúas. Entre D3 e D4 é difícil determinar cal das dúas é mellor. Se calculamos por separado a función d para os parámetros de forma e os de transformación ríxida, comprobamos que a entre D3 e D4 pequena diferencia observada débese principalmente ó axuste dos parámetros de transformación ríxida.
6.1 Preprocesado 81 q i a i q i’ con D 2 q i’ con D 3 q i’ con D 4 ε 1 0,2 0,1 2,174e-01 2,132e-01 2,151e-01 ε 2 0,8 0,1 8,329e-01 8,032e-01 8,054e-01 a 1 20 10 2,003e+01 2,002e+01 2,001e+01 a 2 30 10 3,007e+01 3,002e+01 3,008e+01 a 3 50 10 4,979e+01 4,979e+01 4,978e+01 c x 0 1 2,151e-04 1,040e-03 2,808e-04 c y 0,2 0,1 2,092e-01 1,990e-01 1,990e-01 θ 0,003 0,001 2,361e-03 3,026e-03 2,999e-03 k 0,001 0,001 7,735e-04 8,053e-04 1,133e-03 φ 0 1 6,729e-04 -3,097e-02 1,040e-02 α 0 1 -1,846e-03 -3,072e-02 4,473e-02 β 0 1 7,661e-05 8,999e-06 -1,979e-05 γ 0 1 -8,281e-03 3,180e-02 -4,480e-02 t 1 0 1 1,085e-01 8,941e-02 2,312e-01 t 2 0 1 -2,959e-03 2,875e-03 -7,610e-03 t 3 0 1 -2,020e-02 -4,960e-02 -1,781e-02 Táboa 6.1 Valores dos parámetros q i do superelipsoide da figura 6.4 cos correspondentes factores de peso a i. Á dereita, o promedio dos valores dos parámetros q i’ calculados para as funcións de erro D 2, D 3 e D 4 nun total de 12 probas. Función de erro d d forma d t. ríxida D2 0,787 0,779 0,111 D3 0,266 0,242 0,111 D4 0,319 0,210 0,240 Táboa 6.2 Estimación da calidade do axuste coas distintas medidas de erro.
88 CAPÍTULO 6. RESULTADOS z = 30 z = 50 z = 69 x = 50 y = 50 Figura 6.6 Imaxe de proba con ruído gaussiano de desviación estándar σ = 20. z = 30 z = 50 z = 69 x = 50 y = 50 Figura 6.7 Imaxe de proba con ruído gaussiano de desviación estándar σ = 30. z = 30 z = 50 z = 69 x = 50 y = 50 Figura 6.8 Imaxe de proba con ruído gaussiano de desviación estándar σ = 40.
6.3 Posta en correspondencia 89 z = 10 z = 20 z = 30 z = 41 z = 50 Figura 6.9 Imaxe de proba σruído = 20 tras un filtrado gaussiano con σ = 3 z = 10 z = 20 z = 30 z = 41 z = 50 Figura 6.10 Imaxe de proba σruído = 20 tras un filtrado alfa con α = 0.5 z = 10 z = 20 z = 30 z = 41 z = 50 Figura 6.11 Imaxe de proba σruído = 30 tras un filtrado gaussiano con σ = 3 z = 10 z = 20 z = 30 z = 41 z = 50 Figura 6.12 Imaxe de proba σruído = 40 tras un filtrado alfa con σ = 3
90 CAPÍTULO 6. RESULTADOS Figura 6.13 Cortes z=50 en negativo con filtrados gaussiano σ=3 (esquerda) e alfa α=0.5 (dereita) Figura 6.14 Cortes z=5 en negativo con filtrados gaussiano σ=3 (esquerda) e alfa α=0.5 (dereita) u1 = u2 = 50 u1 = u2 = 20 u1 = 50, u2 = 20 Figura 6.15 Agrupamento de parcelas de superficie con distintos valores dos umbrais maior e menor.
6.3 Posta en correspondencia 91 z = 20 z = 35 z = 50 Figura 6.16 Selección de superficies con K < 0 z = 20 z = 35 z = 50 Figura 6.17 Selección de superficies con K > 0 z = 20 z = 35 z = 50 Figura 6.18 Selección de superficies con K = 0 z = 20 z = 35 z = 50 Figura 6.19 Selección de superficies con H > 0
92 CAPÍTULO 6. RESULTADOS z = 20 z = 35 z = 50 Figura 6.20 Selección de superficies con H < 0 z = 50 x = 50 y = 50 Figura 6.21 Selección de superficies cilíndricas na imaxe de proba con ruído σ = 20. z = 50 x = 50 y = 50 Figura 6.22 Selección de superficies cilíndricas na imaxe de proba con ruído σ = 30. z = 50 x = 50 y = 50 Figura 6.23 Selección de superficies cilíndricas na imaxe de proba con ruído σ = 40.
6.3 Posta en correspondencia 93 z = 10 z = 50 z = 100 z = 140 z = 165 Figura 6.24 Imaxes CT de perna dereita. z = 10 z = 50 z = 100 z = 140 z = 165 Figura 6.25 Imaxes CT de perna dereita filtradas con filtro gaussiano. z = 10 z = 50 z = 100 z = 140 z = 165 Figura 6.26 Agrupamento en parches de superficie
94 CAPÍTULO 6. RESULTADOS z = 1, K > 0 z = 99, K > 0 z = 1, K < 0 z = 99, K < 0 Figura 6.27 Selección de superficies por curvatura gaussiana. z = 1, H > 0 z = 99, H > 0 z = 1, H < 0 z = 99, H < 0 Figura 6.28 Selección de superficies por curvatura media.
6.3 Posta en correspondencia 95 z = 10 z = 50 z = 100 z = 140 z = 165 Figura 6.29 Selección con umbral menor u2=40. Figura 6.30 Reconstrucción triangulada da superficie da Figura 6.29 z = 10 z = 50 z = 100 z = 140 z = 165 Figura 6.31 Selección con umbral menor u2=30. Figura 6.32 Reconstrucción triangulada da superficie da figura 6.30
96 CAPÍTULO 6. RESULTADOS z = 30 z = 50 z = 69 x = 50 y = 50 Figura 6.33 Posta en correspondencia da imaxe de proba de ruído σ = 20 cun modelo de cilindro. z = 30 z = 50 z = 69 x = 50 y = 50 Figura 6.34 Posta en correspondencia da imaxe de proba de ruído σ = 30 cun modelo de cilindro z = 30 z = 50 z = 69 x = 50 y = 50 Figura 6.35 Posta en correspondencia da imaxe de proba de ruído σ = 40 cun modelo de cilindro
6.3 Posta en correspondencia 97 z = 20 z = 35 z = 50 z = 165 x = tx = 98 x = ty = 106 Figura 6.36 Posta en correspondencia entre imaxe e modelo calculado coa función de erro radial z = 20 z = 35 z = 50 z = 165 x = tx = 98 x = ty = 106 Figura 6.37 Posta en correspondencia entre imaxe e modelo calculado coa función de erro D3.
104 BIBLIOGRAFÍA [49]. H. Samet. “ The design and analysis of spatial data structures”. Addison-Wesley, Reading, Massachusetts, 1990. [50]. S. Sandor, R. Leahy. “Surface-based labeling of cortical anatomy using a deformable atlas”. Medical Image, vol. 16, nº 1, 1997, pp 41-54. [51]. F. Solina, R. Bajcsy. “Recovery of parametric models from range images: the case for superquedrics with global and local deformations”. Trans. on PAMI, vol. 12, nº 2, Febreiro 1990, pp 131-146. [52]. L. Staib, A. Chakraborty, J.S. Duncan. “An integrated approach for locating neuroanatomical structure from MRI”. International Journal of Pattern Recognition and Artificial Intelligence, vol. 11, nº 8, 1997, pp 1247-1269. [53]. L. Staib, J.S. Duncan. “ Boundary finding with parametrically deformable models” IEEE Trans. on PAMI, vol. 14, nº 11, 1992, pp 1061-1075. [54]. G. Storvik. “A bayesian approach to dynamic contours through stochastic sampling and simulated annealing”. IEEE Trans. on PAMI, vol. 16, nº 10, Outrubro 1994. [55]. D. Terzopoulos, D. Metaxas. “Dynamic 3D models with local an global deformations: deformable superquadrics”. IEEE Trans. on PAMI, vol. 13, nº 7, Xullo 1991, pp 703-714. [56]. D. Terzopoulos, A. Witkin, M. Kass. “Constraints on deformable models: recovering 3D shape and nonrigid motion” Artificial Intelligence, vol. 36, 1988, pp 91-123. [57]. J. Weickert. “Multiscale texture enhancement”. CAIP95, pp 230-237. [58]. P. Whaite, F.P. Ferrie. “From uncertainty to visual exploration”. IEEE. Trans on PAMI., vol. 13, nº 10, Outubro 1991.