UNIVERSIDADE ESTADUAL PAULISTA FACULDADE DE CIÊNCIAS E TECNOLOGIA Programa de Pós-Graduação em Ciências Cartográficas PRESIDENTE PRUDENTE 2019 MICHELLE SAYURI YANO YWATA EXTRAÇÃO AUTOMÁTICA DE CONTORNOS DE TELHADOS DE EDIFÍCIO NO ESPAÇO-OBJETO INTEGRANDO UM ESTÉREO PAR DE IMAGENS AÉREAS DE ALTA RESOLUÇÃO E MODELOS 3D DE TELHADO MICHELLE SAYURI YANO YWATA EXTRAÇÃO AUTOMÁTICA DE CONTORNOS DE TELHADOS DE EDIFÍCIO NO ESPAÇO-OBJETO INTEGRANDO UM ESTÉREO PAR DE IMAGENS AÉREAS DE ALTA RESOLUÇÃO E MODELOS 3D DE TELHADO Tese de doutorado apresentada ao Programa de Pós-Graduação em Ciências Cartográficas da Universidade Estadual Paulista “Júlio de Mesquita Filho” – Faculdade de Ciências e Tecnologia – Campus de Presidente Prudente, como parte dos requisitos para a obtenção do título de Doutor em Ciências Cartográficas. Orientador: Prof. Dr. Aluir P. Dal Poz Coorientador: Prof. Dr. Milton H. Shimabukuro Y99e Ywata, Michelle Sayuri Yano Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado / Michelle Sayuri Yano Ywata. -- Presidente Prudente, 2019 122 p. Tese (doutorado) - Universidade Estadual Paulista (Unesp), Faculdade de Ciências e Tecnologia, Presidente Prudente Orientador: Aluir Porfírio Dal Poz Coorientador: Milton Hirokazu Shimabukuro 1. Extração de contornos de telhado. 2. Imagem aérea. 3. Dados LASER. 4. Snakes. 5. Programação dinâmica. I. Título. Sistema de geração automática de fichas catalográficas da Unesp. Biblioteca da Faculdade de Ciências e Tecnologia, Presidente Prudente. Dados fornecidos pelo autor(a). Essa ficha não pode ser modificada. DEDICATÓRIA À minha batiam, Mitsuko Tomiyoshi, por todos os ensinamentos e boas lembranças deixadas em meu coração. ありがとうございました。 AGRADECIMENTOS Ao meu grande amor, Fuyu, por todo o apoio, ajuda, ensinamentos e companheirismo durante todos esses anos. Aos meu pais, Celso e Beth, pela educação e oportunidades a mim proporcionadas durante toda a minha vida. Ao meu irmão, Gabriel, pelo apoio e confiança de um irmão caçula. Ao meu orientador Prof. Dr. Aluir Dal Poz, pela orientação, ensinamentos e confiança depositada em mim desde a graduação. Ao meu coorientador Prof. Dr. Milton Shimabukuro, pelos ensinamentos e contribuições dadas para o desenvolvimento deste trabalho. Ao Prof. Dr. Ayman Habib, pela recepção na Purdue University e por todo apoio durante o meu doutorado sanduíche. Aos professores do Departamento de Cartografia e do PPGCC, que contribuíram para minha formação acadêmica, desde a graduação até o doutorado. Aos amigos do PPGCC, principalmente aqueles que estão comigo desde o primeiro ano da graduação e aqueles que fizeram parte do mesmo grupo de pesquisa. Ao meu amigo Prof. Dr. Érico Fernando de Oliveira Martins pelo fornecimento do código de Programação Dinâmica utilizado como base no desenvolvimento do programa implementado neste trabalho. Ao meu amigo Prof. Dr. Henrique Candido de Oliveira pela geração dos mapas de visilibidade utilizados neste trabalho. À FCT/UNESP, ao Programa de Pós-Graduação em Ciências Cartográficas (PPGCC) e à Purdue University por proporcionar todos os meios para o desenvolvimento deste trabalho. O presente trabalho foi realizado com apoio da Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES) – Código de Financiamento 001. 為せば成る 為さねば成らぬ何事も 成らぬは人の為さぬなりけり “Se você tentar, você pode ter sucesso. Se você não tentar, não terá sucesso. Isso é verdade para todas as coisas. Não ter sucesso é o resultado de não tentar.” (Uesugi Harunori) RESUMO Neste trabalho foi proposta uma metodologia para a extração de contornos de telhados de edifícios no espaço-objeto, a partir da integração de um estéreo par de imagens aéreas de alta resolução e modelos 3D aproximados de telhado obtidos a partir de dados de varredura a LASER. Um modelo matemático considerando as propriedades radiométricas e geométricas dos telhados foi formulado a fim de representar o contorno do telhado no espaço- imagem, tendo como base o modelo de contorno ativo Snake. Esse modelo foi então adaptado para descrever os contornos no espaço-objeto considerando um estéreo par de imagens aéreas. Finalmente, o polígono ótimo que representa um dado contorno do telhado foi determinado a partir da otimização, via Programação Dinâmica, da função de energia criada. A solução obtida é uma representação mais acurada para o correspondente contorno do modelo 3D do telhado. O método desenvolvido apresenta também mecanismos para realizar a compensação automática de três tipos de problemas comuns em ambientes urbanos e que podem prejudicar a extração automática de telhados: obstruções perspectivas causadas por edifícios elevados, obstruções diretas causadas por vegetação que se eleva acima do telhado e sombras adjacentes aos telhados, as quais podem ser confundidas com as bordas do telhado. Os experimentos foram realizados utilizando imagens aéreas com GSD ≈ 0,10 m e nuvem de pontos LASER com densidade média de 6 pontos/m2. Os resultados mostraram que o método funciona adequadamente na extração dos contornos de telhados em áreas contendo obstruções e sombras, apresentando valores médios de completeza e correção acima de 90%, e valores médios de RMSE abaixo de 0,5 m para as componentes E e N, e abaixo de 1 m para a componente H. Palavras-chave: Extração de telhados, dados LASER, imagens aéreas, Snake, Programação Dinânima. ABSTRACT In this work a methodology was proposed for extracting building roof contours in the object-space, by integration of a high-resolution aerial images stereo pair and 3D roof models reconstructed from LASER scanning data. A mathematical model considering the radiometric and geometric properties of roofs was developed in order to represent the roof contour in the image-space, based on the Snake active contour model. Then, the model was adapted to represent the contours in the object space considering a stereo pair of aerial images. Finally, the optimal polygon representing a selected roof contour was obtained by optimizing the proposed energy function using Dynamic Programming algorithm. The solution obtained, i.e., a polygon representing each 3D roof contour, will be a higher accurate representation for the correspondent contour of the 3D roof model. The proposed method also presents mechanisms to perform the compensation of three types of common problems in urban environment and which can disturb the automatic roof extraction: perspective occlusions caused by high buildings, occlusions caused by vegetation that covers the roof and shadows that are adjacent to the roofs which can be misinterpreted as roof edges. The experiments were performed using aerial images with GSD ≈ 0,10 m and LASER point cloud with average density of 6 points/m2. The results showed that the proposed method works properly in contour extraction of roofs with occlusion and shadows areas, presenting completeness and correctness average values above 90%, RMSE average values below 0,5 m for E and N components, and below 1 m for H component. Keywords: roof extraction, LiDAR data, aerial images, Snake, Dynamic Programming. LISTA DE FIGURAS Figura 1 – Sistema fotogramétrico. .......................................................................................... 22 Figura 2 – Sistema fotográfico de coordenadas do negativo. ................................................... 22 Figura 3 – Sistema de imagem. ................................................................................................ 23 Figura 4 – Condição de colinearidade entre o CP, o ponto imagem p e o ponto no espaço- objeto P. .................................................................................................................................... 25 Figura 5 – Distorção radial simétrica. ...................................................................................... 26 Figura 6 – Efeitos da distorção radial simétrica. ...................................................................... 26 Figura 7 – Distorção descentrada. ............................................................................................ 27 Figura 8 – Perfil dos tipos de bordas existentes em uma imagem. ........................................... 29 Figura 9 – Detecção de bordas por operadores de derivação. .................................................. 30 Figura 10 – Máscara de gradiente. (a) Região de tamanho 3x3 em uma imagem; (b) Operadores de Sobel; (c) Operadores de Prewitt. ............................................................... 32 Figura 11 – Linha do tempo da criação dos principais operadores de interesse. ..................... 33 Figura 12 – Janelas 3×3 e os vetores nas direções 0o, 45o, 90o e 135o. .................................... 34 Figura 13 – (a) Visada em perspectiva central. (b) Visada ortogonal. ..................................... 37 Figura 14 – Ilustração do efeito de duplo mapeamento............................................................ 38 Figura 15 – (a) Imagem original em projeção perspectiva. (b) Efeito de duplo mapeamento. (c) Área de oclusão (região em vermelho). .............................................................................. 39 Figura 16 – Situações em que ocorrem múltiplos retornos. ..................................................... 41 Figura 17 – Componentes de um sistema de VLA e suas relações geométricas. ..................... 43 Figura 18 – Perfil gerado pela conexão de pontos LASER sequenciais numa região de vegetação. ................................................................................................................................. 49 Figura 19 – Exemplos de elementos estruturantes. .................................................................. 50 Figura 20 – Efeitos da erosão utilizando o EE octógono de tamanho 21 pixels. ..................... 52 Figura 21 – Efeitos da dilatação utilizando o EE octógono de tamanho 21 pixels. ................. 54 Figura 22 – Exemplo da transformação top-hat por abertura em um sinal unidimensional. ... 56 Figura 23 – Exemplo da transformação top-hat por fechamento em um sinal unidimensional. .................................................................................................................................................. 56 Figura 24 – Representação geométrica de uma curva Snake.................................................... 61 Figura 25 – Problema clássico de Programação Dinâmica. ..................................................... 66 Figura 26 – Fluxograma do método proposto. ......................................................................... 72 Figura 27 – (a) Amostra do mapa de visibilidade, em azul a região oclusa. (b) Imagem aérea correspondente. ......................................................................................................................... 74 Figura 28 – Etapas da classificação da nuvem de pontos LASER. .......................................... 75 Figura 29 – Pontos de vegetação. ............................................................................................. 76 Figura 30 – Ilustração do elemento estruturante octógono de tamanho 3. ............................... 77 Figura 31 – Exemplo de imagem de vegetação gerada pelo método proposto. ....................... 77 Figura 32 – Exemplo de imagem de sombra gerada a partir de uma imagem sem suavização. (a) Detalhes de uma imagem aérea de alta resolução. (b) Sombras detectadas (em branco). .. 78 Figura 33 – Imagem de sombra gerada pelo método proposto................................................. 79 Figura 34 – Fator de compensação de sombra. ........................................................................ 82 Figura 35 – Ideia utilizada na criação do termo para o cálculo da resposta de quina. ........ 84 Figura 36 – Princípio da ambiguidade. ..................................................................................... 87 Figura 37 – Criação do espaço de busca 1D no espaço 3D. ..................................................... 89 Figura 38 – Criação do espaço de busca 2D no espaço 3D. ..................................................... 90 Figura 39 – Criação do espaço de solução no espaço 3D. ........................................................ 92 Figura 40 – Aplicação da PD. ................................................................................................... 95 Figura 41 – Gráficos do número total de combinações obtidas em cada algoritmo a partir de diferentes quantidades de pontos no espaço de solução. ........................................................ 101 Figura 42 – Áreas utilizadas no Experimento 1 contendo telhados cobertos por vegetação. . 102 Figura 43 – Resultado visual do Experimento 1. ................................................................... 104 Figura 44 – Área utilizada no Experimento 2 contendo sombras adjacente aos telhados. ..... 106 Figura 45 – Resultado visual do Experimento 2. ................................................................... 108 Figura 46 – Áreas utilizadas no Experimento 3 contendo obstruções perspectivas. .............. 110 Figura 47 – Resultado visual do Experimento 3. ................................................................... 112 Figura 48 – Detalhe da mancha branca no telhado que atrapalhou o delineamento correto da borda, contorno detectado (verde) e contorno correto (rosa tracejado). ................................. 113 LISTA DE TABELAS Tabela 1 – Tabela da matriz de “custos parciais”. .................................................................... 67 Tabela 2 – Tabela da matriz de "custos parciais" g2(e2,e3). ...................................................... 68 Tabela 3 – Tabela da matriz de "custos máximos" f1(e2).......................................................... 68 Tabela 4 – Tabela da matriz de "custos acumulados" f1(e2) + g2(e2,e3). ................................... 68 Tabela 5 – Tabela da matriz de "custos máximos". .................................................................. 69 Tabela 6 – Última tabela da matriz de "custos máximos". ....................................................... 69 Tabela 7 – Caminho inverso com a superposição de 1 nó........................................................ 69 Tabela 8 – Parâmetros e limiares utilizados. ............................................................................ 98 Tabela 9 – Testes de comparação dos algoritmos. ................................................................. 100 Tabela 10 – Comparação entre a quantidade de pontos no espaço de solução e número de combinações resultantes em cada algoritmo. ......................................................................... 101 Tabela 11 – Parâmetros de qualidade dos contornos iniciais do Experimento 1. .................. 105 Tabela 12 – Parâmetros de qualidade dos contornos extraídos pelo Experimento 1. ............ 105 Tabela 13 – Parâmetros de qualidade dos contornos iniciais do Experimento 2. .................. 109 Tabela 14 – Parâmetros de qualidade dos contornos extraídos pelo Experimento 2. ............ 109 Tabela 15 – Parâmetros de qualidade dos contornos iniciais do Experimento 3. .................. 113 Tabela 16 – Parâmetros de qualidade dos contornos extraídos pelo Experimento 3. ............ 113 LISTA DE SIGLAS 1D Unidimensional 2D Bidimensional 3D Tridimensional BTH Black Top-hat CCD Charge Couple Device CP Centro perspectivo EE Elemento estruturante EI Espaço-imagem RMSE Raiz do erro médio quadrático EO Espaço-objeto FCT Faculdade de Ciências e Tecnologia GNSS Global Navigation Satellite System GSD Ground Sample Distance IFOV Instantaneous Field of View IMU Inertial Measurement Unit INS Inertial Navigation System ISPRS International Society for Photogrammetry and Remote Sensing LASER Light Amplification by Stimulated Emission of Radiation LPS Leica Photogrammetry Suite MDS Modelo digital de superfície MDT Modelo digital de terreno MM Morfologia matemática PD Programação Dinâmica POE Parâmetros de orientação exterior POI Parâmetros de orientação interior PP Ponto principal SGBM Surface-Gradient-Based Method TIM Triangulated Irregular Network TM Transverso de Mercator UNESP Universidade Estadual Paulista UTM Universal Transverso de Mercator VLA Varredua a LASER aerotransportado WGS84 World Geodetic System 1984 WTH White Top-hat SUMÁRIO 1. INTRODUÇÃO ............................................................................................................... 14 1.1. Considerações iniciais ................................................................................................ 14 1.2. Hipótese e Objetivos .................................................................................................. 17 1.3. Justificativa ................................................................................................................ 18 1.4. Estrutura do trabalho .................................................................................................. 19 2. FUNDAMENTAÇÃO TEÓRICA ................................................................................. 21 2.1. Sistemas de coordenadas ........................................................................................... 21 2.2. Conceito de borda ...................................................................................................... 28 2.3. Conceito de quina ...................................................................................................... 33 2.4. Conceito de sombra ................................................................................................... 35 2.5. Áreas de oclusão e mapa de visibilidade ................................................................... 36 2.6. Sistema de varredura a LASER aerotransportado (SVLA) ....................................... 40 2.7. Morfologia Matemática ............................................................................................. 50 2.8. Modelo de contorno ativo: Snakes ............................................................................. 57 2.9. Programação Dinâmica .............................................................................................. 65 3. MATERIAL E MÉTODOS ............................................................................................ 70 3.1. Material ...................................................................................................................... 70 3.2. Métodos ..................................................................................................................... 71 3.3. Forma de análise dos resultados ................................................................................ 95 4. RESULTADOS ................................................................................................................ 97 4.1. Introdução .................................................................................................................. 97 4.2. Parâmetros e limiares utilizados ................................................................................ 97 4.3. Resultados e análise ................................................................................................... 98 5. CONCLUSÕES E RECOMENDAÇÕES ................................................................... 114 REFERÊNCIAS ................................................................................................................... 117 Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 14 1. INTRODUÇÃO 1.1. Considerações iniciais Atualmente a humanidade está em uma era na qual as informações espaciais desempenham um importante papel para o público em geral. No entanto, o grande desafio continua sendo fornecer informações de alta qualidade, para permitir o processamento geoespacial avançado e apoiar a solução colaborativa de problemas (CHEN et al., 2016a). Diversas questões urbanas como a modelagem e visualização da poluição, atualizações cadastrais, estimativa de volume predial, planejamento urbano e projeto de redes de comunicação tem utilizado os modelos tridimensionais das cidades como base de estudo (BENCIOLINI et al., 2018). A reconstrução tridimensional dos edifícios presentes nos modelos de cidades envolve as etapas de detecção, extração do contorno e reconstrução do telhado. Embora muitos esforços tenham sido feitos nas últimas décadas para o desenvolvimento de novos métodos de extração de edifícios, a extração exata e automática ainda continua sendo um desafio para as comunidades de Sensoriamento Remoto e Visão Computacional (HUANG et al., 2019). Isso se deve às condições de complexidade e variabilidade das cenas, que torna desafiadora a questão de extração de telhados no âmbito do reconhecimento de objetos. Dessa forma, se torna interessante o desenvolvimento de metodologias que explorem a sinergia de diferentes fontes de dados para realizar essa tarefa, visto que cada fonte possui as suas vantagens e desvantagens. Diversos autores (ROTTENSTEINER et al., 2014; GERKE e XIAO, 2014; ZHOU e ZHOU, 2014; WU et al., 2015; OLIVEIRA, 2016a; LARI et al., 2017) têm apresentado a integração entre imagem e dados de varredura a LASER (Light Amplification by Stimulated Emission of Radiation) aerotransportado (VLA) para a extração de edifícios. Essa combinação é atrativa, uma vez que esses dois tipos de dados possuem informações complementares. Quando se trata da obtenção de planos de telhados e sua orientação, os dados de VLA fornecem melhores resultados, pois eles apresentam informações altimétricas mais precisas. Já com relação à extração de contornos de telhados, as imagens ópticas fornecem melhores resultados se comparados com dados de VLA de baixa ou média densidade, uma vez que elas possuem melhores informações posicionais ao longo das linhas de quebra (DAL POZ et al., 2009). Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 15 No entanto, a questão de como integrar de forma eficiente essas duas fontes de dados para a extração de edifícios, ainda precisa ser resolvida e continua sendo estudada por diversos autores (CHEN et al., 2016a). No trabalho de Gerke e Xiao (2014) uma nova abordagem para o processo de segmentação de áreas urbanas é apresentada, integrando a qualidade geométrica dos dados de VLA e as informações espectrais presentes em imagens, para detectar edifícios, árvores e terreno. Os atributos de geometria e textura da imagem são relacionados a cada ponto da nuvem de pontos LASER, formando um conjunto de voxels. Esse conjunto é então classificado a partir de duas abordagens diferentes (supervisionada e não supervisionada) com o objetivo de detectar edifícios, árvores, terrenos cobertos por vegetação e terrenos descobertos. Zhou e Zhou (2014) apresentam uma metodologia para a extração de edifícios através da elaboração de um grafo de atributos baseado nas características dos edifícios (geometria, estrutura, forma), e construído a partir da combinação das imagens aéreas e dos dados de VLA. Lari et al. (2017) propõem um método para melhorar a interpretabilidade das superfícies planares extraídas de dados de VLA, a partir do uso de informações de textura provenientes de imagens coletadas através de sistemas de mapeamento moderno de baixo custo. Primeiro as superfícies planares são extraídas dos dados de VLA através de um novo processo de segmentação. Em seguida, as partes não oclusas das superfícies planares extraídas são decompostas em segmentos e um processo de renderização é realizado a fim de texturizar esses segmentos usando as imagens disponíveis. Como pode ser observado nos exemplos apresentados anteriormente, as metodologias de extração de edifícios existentes na literatura exploram diversas técnicas para alcançar o objetivo desejado. Embora existam poucas abordagens baseadas na modelagem Snake, essa teoria vem ganhando espaço no campo da extração de feições. Yan et al. (2015) realiza a construção de modelos 3D de edifícios a partir de dados de VLA com base em um novo modelo para o algoritmo Snake, denominado “algoritmo Snake 2D”, o qual minimiza a função de energia associada à topologia planar através da técnica de redução de grafos. Essa técnica consiste em um conjunto de operações gráficas capazes de reduzir uma topologia planar bidimensional a vértices isolados, enquanto mantém a energia mínima da topologia. No trabalho de Chen et al. (2016b), um modelo melhorado de Snake é proposto para refinar os contornos iniciais de edifícios extraídos de dados de VLA. Primeiro, um modelo Snake aprimorado é construído com as restrições de ângulo, gradiente e área. Em seguida, os nós do contorno são movidos para encontrar o melhor resultado utilizando o algoritmo guloso (ZIVIANI, 2011). Griffiths e Boehm (2019) apresentam um método para aumentar a Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 16 qualidade de grandes conjuntos de dados com baixa precisão, para permitir o treinamento de redes de segmentação. Os pontos centrais de cada edifício são calculados e usados como ponto semente para a aplicação do modelo Snake. O modelo é utilizado para determinar as bordas dos edifícios e realizar uma segmentação precisa. A imagem segmentada é então utilizada como dado de treinamento para a aplicação em redes neurais convolucionais. De uma forma geral, os métodos de extração dos contornos de telhados existentes na literatura, tanto aqueles que usam imagem ou dados de VLA quanto aqueles que usam os dois tipos de dados, podem trabalhar de três diferentes formas. No primeiro caso, têm-se os métodos que extraem os contornos de telhados diretamente na imagem aérea (COTE e SAEEDI, 2013). O segundo caso envolve métodos que extraem contornos tridimensionais de telhados utilizando nuvem de pontos em 3D de VLA (SEO et al., 2014; CHEN et al., 2014). E, no terceiro caso, têm-se os métodos que combinam de alguma forma ambas fontes de dados: imagem e nuvem de pontos LASER (GERKE e XIAO, 2014; ZHOU e ZHOU, 2014; GILANI et al., 2015; OLIVEIRA, 2016a). Este trabalho inspira-se inicialmente na ideia proposta em Fazan (2014), o qual realiza a extração dos contornos de telhados diretamente no espaço-objeto (3D), integrando numa única função (função de energia Snake) uma imagem aérea e dados de VLA. Esse método de extração combina as propriedades radiométricas do telhado, obtidas a partir da imagem aérea, e as propriedades geométricas do telhado, obtidas a partir dos dados de VLA. A presente proposta diferencia-se principalmente por usar: 1) um estéreo par de imagens aéreas, em vez de apenas uma imagem; 2) modelos tridimensionais de telhado previamente reconstruídos de dados de VLA, em vez de dados brutos dessa mesma fonte de dados; e 3) informações contextuais da cena, obtidas a partir das imagens e dos dados de VLA. Além desses aspectos diferenciadores, apresenta-se um método totalmente automático, contrapondo-se à estratégia semi-automática proposta em Fazan (2014). Vale ressaltar que embora não sejam utilizados dados de VLA em sua forma bruta, mas sim como modelos 3D de telhados, o método proposto para o delineamento do contorno do telhado é considerado totalmente automático visto que não há necessidade de interferência humana em nenhuma das etapas do método de extração. Uma grande vantagem do método proposto é a possibilidade de tratar obstruções perspectivas, que são dependentes do ponto de vista de tomada das imagens e que não podiam ser tratadas pelo método proposto em Fazan (2014). Outro diferenciador é a incorporação de mecanismos para superar problemas comuns quando se trata de extração de contornos de Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 17 telhados em ambientes urbanos, sendo eles: obstruções causadas por vegetação que se eleva acima dos telhados; e sombras adjacentes aos telhados, as quais podem ser confundidas com as bordas dos telhados. O método proposto visa aproveitar a sinergia entre dados de VLA e imagens ópticas no contexto de extração de edifícios, isto é: planos de telhados são melhor definidos em dados de VLA e contornos de edifícios são melhor definidos em imagens aéreas de alta resolução (Ground Sample Distance (GSD) ≈ 10 cm). Isso será originalmente realizado através de uma função de energia baseada na formulação Snake, cuja solução (polígonos representando contornos refinados de edifícios) será obtida por otimização via Programação Dinâmica (PD). 1.2. Hipótese e Objetivos 1.2.1. Hipótese Antes da apresentação da hipótese que fundamenta este trabalho, é necessário salientar que o alvo de estudo são contornos de telhados formados por segmentos de retas conectados em 90°. A partir dessa premissa, apresenta-se a seguinte hipótese: “A partir da otimização via PD de uma função Snake com parâmetros adaptativos, é possível, de forma robusta, melhorar os contornos 3D de telhados previamente obtidos de dados de VLA, utilizando imagens aéreas de alta resolução e informações contextuais extraídas das imagens e dados de VLA, atingindo melhor acurácia dos contornos de telhado no espaço-objeto.” 1.2.2. Objetivo Geral Esse trabalho tem por objetivo principal desenvolver um método para a extração de contornos de telhados de edifícios no espaço-objeto (3D), a partir do uso integrado de um estéreo par de imagens aéreas de alta resolução (GSD ≈ 10 cm), modelos tridimensionais aproximados de telhado obtidos a partir de dados de VLA, e informações contextuais obtidas previamente das imagens e dados de VLA, utilizando a modelagem Snake com parâmetros adaptativos e solução em 3D do modelo matemático via otimização por PD. Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 18 1.2.3. Objetivos Específicos 1) Desenvolver uma função de energia baseada na teoria de contorno ativo (Snake) para modelar contornos de telhado de edifícios no espaço-objeto (3D), tendo como base os atributos radiométricos e geométricos de telhados derivados de um estéreo par de imagens aéreas, as relações matemáticas entre os espaços imagem e objeto, e as propriedades geométricas presentes em modelos tridimensionais de telhados; 2) Desenvolver um mecanismo para a compensação automática de obstruções perspectiva sobre telhados, causadas por edifícios vizinhos elevados; 3) Desenvolver um mecanismo para a compensação automática de obstruções diretas sobre telhados, causadas por partes da vegetação que se eleva acima dos telhados; 4) Desenvolver um mecanismo para a compensação automática de sombreamento nas adjacências do edifício, causados por ele próprio e que podem ser confundidos com as bordas de telhados; 5) Desenvolver uma solução para eficientemente otimizar por PD a função de energia desenvolvida; 6) Avaliar experimentalmente a metodologia proposta através da comparação entre os resultados obtidos pelo método e os dados obtidos por um operador humano. 1.3. Justificativa A relevância científica e tecnológica deste trabalho é evidenciada pela importância dada ao tema pela International Society for Photogrammetry and Remote Sensing (ISPRS), na Comissão II (Photogrammetry) que possui como um de seus objetivos realizar de forma integrada os processos de correspondência (matching) e delineamento de objetos. Segundo Chen et al. (2016a), essa integração ainda precisa ser explorada com mais profundidade, uma vez que, atualmente, esses processos são desenvolvidos de forma independente uma da outra, e assim, os erros na solução de uma das tarefas se propaga para a próxima. Dessa forma, a solução simultânea, como é proposta neste trabalho, pode ajudar a superar algumas dessas limitações existentes nos métodos atuais. Como também apontam Chen et al. (2016a), a fusão de dados de diferentes sensores, como dados ópticos e não ópticos, que também é o caso desse trabalho, é atualmente um tópico de grande relevância. Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 19 O estudo dos processos de extração automática de objetos em áreas urbanas é de grande relevância tecnológica, visto que, atualmente, a extração ainda depende da participação de um operador humano, o que demanda muito tempo e elevados custos no processo de aquisição dos dados. Sendo assim, este trabalho se insere no contexto de desenvolvimento de novas metodologias para a eficiente aquisição de informações espaciais a partir do uso de dados provenientes de diferentes sensores. O ineditismo deste trabalho está no fato de se desconhecer na literatura relacionada, metodologias de extração de contornos de telhados de edifícios no espaço-objeto (3D) que utilizem um estéreo par de imagens aéreas juntamente com o modelo tridimensional aproximado do telhado. O grande desafio dessa proposta está relacionado principalmente ao Objetivo Específico 1, o qual propõe o desenvolvimento de uma função de energia baseada na modelagem Snake para a representação de contornos de telhados no espaço-objeto, utilizando a integração de um estéreo par de imagens aéreas, modelo 3D do telhado e informações contextuais da cena. Além disso, constitui também um importante desafio (Objetivo Específico 5) o desenvolvimento de uma estratégia para eficientemente otimizar esta função de energia usando PD. Outros desafios de difícil solução estão relacionados com a presença de vários tipos de obstrução (Objetivos Específicos 2, 3 e 4), os quais requerem estratégias especiais para compensação de seus efeitos. Como será descrito na proposta do método, esta compensação fará parte do próprio processo de otimização, caracterizando um método com solução simultânea e não hierárquica, sendo outro aspecto original desta proposta. 1.4. Estrutura do trabalho Este trabalho está divido em cinco capítulos. O Capítulo 1 apresentou uma contextualização sobre o assunto tratado nesta tese, bem como a hipótese, objetivos e justificativa do mesmo. O Capítulo 2 traz a fundamentação teórica necessária para o desenvolvimento da metodologia proposta, onde são apresentados os sistemas de coordenadas do espaço-imagem e espaço-objeto; os conceitos de quina, borda e sombra; o conceito de ortorretificação, áreas de oclusão e mapa de visibilidade; os fundamentos de um sistema de varredura a LASER aerotransportado; os conceitos básicos de morfologia matemática; o modelo de contorno ativo Snake; e o algoritmo de PD. Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 20 No Capítulo 3 são apresentados os dados e recursos de hardware e software utilizados, bem como o método proposto para a extração dos contornos de telhados, o qual foi dividido em três partes. Primeiramente é apresentado o processo de geração das imagens auxiliares contendo as informações contextuais da cena de interesse. Após isso, é apresentada a formulação da função de energia Snake considerando o estéreo par de imagens. Por fim, é descrita a estratégia utilizada para a solução da função de energia via PD. Nesse capítulo é apresentado também a forma de análise dos resultados. O Capítulo 4 apresenta os parâmetros e limiares empregados na implementação do método e os resultados experimentais obtidos pelo método proposto. Os experimentos foram divididos em três grupos, onde cada um teve o objetivo de demonstrar a eficiência do método em tratar os problemas de obstruções e sombreamento. Também foram realizados testes para analisar a eficiência do algoritmo de PD na solução do método proposto. Por fim, o Capítulo 5 apresenta as conclusões do trabalho realizado e algumas recomendações para trabalhos futuros. Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 21 2. FUNDAMENTAÇÃO TEÓRICA O presente capítulo apresenta uma revisão teórica sobre conceitos e técnicas fundamentais para o desenvolvimento do estudo proposto. A Seção 2.1 apresenta os sistemas de coordenadas do espaço-imagem e espaço-objeto, bem como o procedimento para a transformação de coordenadas entre os dois sistemas. Na Seção 2.2 é apresentada a definição de borda e os principais conceitos sobre o cálculo do gradiente em pontos de borda. A Seção 2.3 aborda o conceito de quina e o operador de Moravec. A Seção 2.4 trata sobre as sombras presentes em imagens aéreas e sobre os métodos de detecção existentes. Na Seção 2.5 é apresentado o conceito de ortorretificação, áreas de oclusão e mapas de visibilidade. A Seção 2.6 trata sobre os fundamentos de um sistema de VLA, o seu funcionamento e os processos de filtragem e classificação de dados de VLA, apresentando as características básicas de um algoritmo de filtragem e a ideia fundamental do funcionamento de um algoritmo de classificação baseado em curvatura. A Seção 2.7 traz os conceitos básicos da morfologia matemática e dos operadores de erosão, dilatação, abertura e fechamento. A Seção 2.8 apresenta os fundamentos do modelo de contorno ativo Snake. E, por fim, a Seção 2.9 trata dos fundamentos do algoritmo de otimização de Programação Dinâmica. 2.1. Sistemas de coordenadas 2.1.1. Sistemas de coordenadas do espaço-imagem O espaço-imagem é considerado como sendo o espaço que envolve o nodo posterior do sistema de lentes da câmara e o plano do negativo (LUGNANI, 1987). A seguir são apresentados os principais sistemas de coordenadas do espaço imagem: sistema com origem no centro da imagem, sistema fotogramétrico e sistema de imagem. 1) Sistema com origem no centro da imagem: é um sistema bidimensional definido da seguinte forma: x Origem coincidente com o centro da imagem (intersecção entre as diagonais da imagem), cuja posição é dada por (Cx, Cy); x O eixo xc tem a direção horizontal mais próxima da linha de voo, sendo as coordenadas positivas nesse sentido; Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 22 x O eixo yc tem sentido positivo a 90o (sentido anti-horário) do eixo xc. 2) Sistema Fotogramétrico: é um sistema tridimensional definido como segue (Figura 1): x Origem no Centro Perspectivo (CP) da câmara; x Os eixos x e y são paralelos e de mesmo sentido dos eixos xc e yc do Sistema com origem no centro da imagem; x O eixo z é perpendicular ao plano da imagem, de forma a tornar o sistema dextrogiro. Nesse sistema a coordenada z de qualquer ponto possui o valor de -f. Figura 1 – Sistema fotogramétrico. Fonte: Lugnani (1987) Para o caso do negativo, a origem desse sistema coincide com o CP, os eixos x e y sofrem uma reflexão em relação ao sistema anterior e o eixo z, é tal que, torna o sistema dextrogiro (Figura 2). Neste caso, a coordenada z de qualquer ponto situado no plano do negativo assumirá o valor de f. Figura 2 – Sistema fotográfico de coordenadas do negativo. Fonte: Lugnani (1987) Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 23 3) Sistema de Imagem: é um sistema bidimensional associado à imagem digital (Figura 3). Uma imagem digital é constituída por um conjunto de elementos de imagem, denominado pixels, espacialmente ordenados em forma de matriz. Cada pixel possui um atributo que indica o nível de cinza (g(c,l)), e a sua posição é dada por um par ordenado (coluna, linha) ou (c, l). A definição desse sistema é apresentada a seguir: x Origem no centro do pixel situado no canto superior esquerdo; x O eixo l coincide com a primeira coluna da imagem; x O eixo c coincide com a primeira linha da imagem. Figura 3 – Sistema de imagem. 2.1.2. Sistemas de coordenadas do espaço-objeto – Sistema Universal Transverso de Mercator (UTM) O sistema UTM é obtido a partir da projeção dos pontos do elipsoide de revolução sobre um cilindro secante, que por sua vez é desenvolvido em um plano. Esse sistema faz parte da categoria de projeções conformes, cuja propriedade principal é a preservação da forma de pequenas áreas. Nessa projeção as coordenadas de um ponto são representadas por (E, N) e podem ser calculadas a partir das coordenadas no sistema Transverso de Mercator (TM) através das seguintes equações: Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 24 0 1 0 1 [ ] 0 1 (Hemisfério sul) (2.1) ou 0 1 0 1 [ ] 0 1 (Hemisfério norte) (2.2) Como as coordenadas no Sistema UTM são referenciadas a um sistema cartesiano, é necessária uma terceira componente para se ter a determinação unívoca de um ponto sobre a superfície física. Essa terceira componente pode ser a altitude geométrica (h), gerando a combinação UTM-h, ou a altitude ortométrica (H), gerando a combinação UTM-H. Essas combinações constituem sistemas caracterizados como híbrido, por serem compostos por superfícies de referência de diferentes características. As coordenadas E e N são coordenadas cartesianas sobre o cilindro, que é desenvolvível no plano, e h e H são obtidas a partir das distâncias dos pontos da superfície física a uma superfície não plana. Em alguns procedimentos essas combinações provocam erros sistemáticos, principalmente no caso de extensas áreas. 2.1.3. Transformação de um ponto no espaço-objeto para o espaço-imagem O modelo matemático comumente utilizado para relacionar pontos no espaço-objeto com os respectivos pontos no espaço-imagem são as Equações de Colinearidade. Este modelo consiste basicamente em duas equações (Equação 2.3), as quais consideram que um ponto no espaço-objeto P(X, Y, Z), o seu correspondente no espaço-imagem p(x, y) e o centro perspectivo CP(XCP, YCP, ZCP) são colineares no instante da tomada da imagem (Figura 4). Sendo assim, a transformação das coordenadas do espaço-objeto para o espaço-imagem se inicia transformando o ponto no espaço-objeto para o Sistema Fotogramétrico, a partir da seguinte equação (ANDRADE, 1998): ( ) ( ) ( ) ( ) ( ) ( ) (2.3) ( ) ( ) ( ) ( ) ( ) ( ) Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 25 onde (x, y) são as coordenadas de um ponto p no Sistema Fotogramétrico; f é a distância focal da câmara; , e , são os elementos da matriz de rotação; (XL, YL, ZL) são as coordenadas do ponto P no sistema do espaço-objeto; (XCP, YCP, ZCP) são as coordenadas do CP no sistema do espaço-objeto. Figura 4 – Condição de colinearidade entre o CP, o ponto imagem p e o ponto no espaço- objeto P. Após isso, o ponto p(x, y) no Sistema Fotogramétrico é transformado para o referencial da imagem (c, l), sendo necessário primeiro introduzir as distorções relacionadas ao sistema de lentes. Isso é necessário visto que a imagem digital original, onde geralmente as operações de análise de imagem são realizadas, é afetada por tais distorções. Em se tratando de imagens capturadas por câmaras digitais, as distorções mais relevantes a se considerar são as distorções das lentes, sendo elas: distorção radial simétrica e distorção descentrada. A distorção radial simétrica é provocada pela dificuldade em fabricar superfícies de lentes com paraboloides de revolução perfeitos. Isso provoca um desvio na trajetória do raio de luz emergente no ponto nodal posterior, em relação ao raio incidente no ponto nodal anterior. Mais formalmente, a distorção radial simétrica pode ser encarada como sendo a componente radial indesejável da refração sofrida por um raio de luz ao atravessar uma lente ou um sistema de lentes (ANDRADE, 1998). Dessa forma, um raio de luz ao atravessar um sistema de lentes com um ângulo de incidência α com o eixo óptico, refrata-se com valor α + δα, provocando um deslocamento ∆r no ponto imagem (Figura 5). Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 26 Figura 5 – Distorção radial simétrica. Fonte: adaptado de Andrade (1998). Essa distorção provoca um deslocamento radial e simétrico com relação ao ponto principal, por isso é denominada de distorção radial simétrica. Vale salientar, que esse efeito deforma geometricamente as imagens. A Figura 6 ilustra as deformações que essa distorção causa nas imagens. A Figura 6(a) apresenta um quadrado sem distorções, a Figura 6(b) ilustra a distorção positiva chamada de distorção em almofada e a Figura 6(c) mostra a distorção negativa, chamada de distorção em barrilete. Figura 6 – Efeitos da distorção radial simétrica. (a) Quadrado sem distorções; (b) Distorção em almofada; e (c) Distorção barrilete. A distorção descentrada é causada pela impossibilidade do fabricante em alinhar perfeitamente os eixos ópticos das lentes que compõe o sistema de lentes. Ela provoca o deslocamento de um ponto pertencente à imagem, normal à linha radial passante pelo ponto (Figura 7). Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 27 Como pode ser visto na Figura 7, a imagem d está em sua posição correta tangencialmente. Já o ponto e foi deslocado no sentido horário de sua posição correta em e’, e o ponto f foi deslocado no sentido anti-horário de sua posição correta f’. Como exemplo, se a imagem de uma reta e’f’ passa por O (ponto principal), a distorção fará com que a reta seja uma linha curva, e isto só não ocorrerá quando a linha coincidir com a linha de distorção nula. Figura 7 – Distorção descentrada. O modelo matemático que permite corrigir ambas as distorções é dado por (LUGNANI, 1987): ( ) ( ) (2.4) ( ) ( ) onde k1, k2, k3, P1, P2 são parâmetros conhecidos de calibração; √( ) ( ) ; (x’, y’) são as coordenadas do ponto imagem p afetadas pelas distorções da lente e (x, y) são as coordenadas do ponto p livre das distorções da lente, sendo essas coordenadas obtidas pelas Equações de Colinearidade (Equação 2.3). A Equação 2.4 permite o cálculo das coordenadas não afetadas pelos erros sistemáticos (x, y) em função das coordenadas afetadas pelos erros sistemáticos (x’, y’). No entanto, como esta equação não admite a inversa analítica, é necessário utilizar alguma solução numérica para encontrar os valores de (x’, y’) que satisfaçam a Equação 2.4. Isso pode ser eficientemente resolvido através de um método numérico para encontrar raízes de funções, com por exemplo o Newton-Raphson (SULI e MAYERS, 2003). Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 28 O ponto resultante p(x’, y’) da introdução das distorções das lentes necessita agora ser transformado para o referencial central da imagem (xc, yc). O referencial central da imagem é paralelo ao referencial fotogramétrico bidimensional, deferindo apenas na origem. A origem do sistema fotogramétrico bidimensional é o ponto principal (PP) e o sistema central possui origem no centro da imagem (Ci). As coordenadas do PP(x0, y0) no referencial central da imagem são determinadas em processos de calibração de câmeras aéreas. Portanto, sendo (xc, yc) as coordenadas do ponto p(x’, y’) no referencial central, tem- se: (2.5) Finalmente, o último passo consiste em transformar o ponto p(xc, yc) para seu correspondente ponto p(c, l) no referencial digital da imagem, como mostrado a seguir: (2.6) onde: - ; - ; - DC e DL são as dimensões horizontal e vertical (em milímetros) de um pixel. 2.2. Conceito de borda 2.2.1. Definição de borda As bordas oferecem importantes informações dos objetos imageados, já que as mesmas correspondem a descontinuidades fotométricas e geométricas. Segundo Gonzalez e Woods (2000), a borda é o limite entre duas regiões com propriedades dinstintas de nível de cinza. A Figura 8 apresenta o perfil dos tipos de bordas presentes em uma imagem, sendo elas (JAIN et al., 1995): Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 29 Borda degrau: ocorre quando a intensidade muda abruptamente de um valor em um dos lados da descontinuidade para outro valor no outro lado da descontinuidade. Borda rampa: é um caso especial da borda degrau, onde a mudança de intensidade não é instantânea, ou seja, acontece de forma gradual em uma distância finita. Borda platô: decorre de uma mudança abrupta de intensidade, voltando à intensidade anterior após um uma curta distância. Borda telhado: é um caso especial da borda platô, onde a mudança de intensidade não é simultânea, ocorrendo em uma distância finita. Figura 8 – Perfil dos tipos de bordas existentes em uma imagem. Fonte: Adaptado de JAIN et al., 1995. A partir da Figura 9 é possível observar que as bordas estão localizadas nos máximos ou mínimos da primeira derivada, e nos pontos de cruzamento por zero da segunda derivada. Além disso, nota-se que a segunda derivada é positiva na parte da transição associada ao lado escuro da borda e negativa na parte da transição associada ao lado claro da borda. Sendo assim, a magnitude da primeira derivada pode ser utilizada para detectar a presença de uma borda na imagem, enquanto o sinal da segunda derivada poder ser utilizado para indicar se um pixel de borda se encontra no lado escuro ou claro da mesma (GONZALEZ e WOODS, 2000). Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 30 Figura 9 – Detecção de bordas por operadores de derivação. Fonte: GONZALEZ e WOODS, 2000. 2.2.2. Gradientes de borda Segundo Gonzalez e Woods (2000), a maioria das técnicas de detecção de bordas são baseadas em um operador local diferencial. O método mais comum de diferenciação é o do gradiente. Considerando uma função f(x,y), por definição, sabe-se que o vetor gradiente dessa função irá apontar na direção onde há uma mudança mais rápida de f na posição (x,y). Sendo assim, em regiões de bordas o vetor gradiente irá apontar na direção transversal a essas bordas (GONZALEZ e WOODS, 2000). Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 31 O gradiente de f nas direções (x,y) é dado pelo vetor: =[ ] [ ] (2.7) A magnitude desse vetor, denotada por ‖ ‖, é uma quantidade muito importante na detecção de bordas. Ela equivale à maior taxa de aumento de f(x,y) por unidade de distância na direção de , e é dada por: ‖ ‖ ( ) [ ] ⁄ (2.8) Geralmente aproxima-se a magnitude do vetor gradiente com valores absolutos: ‖ ‖ | | | | (2.9) Outra quantidade importante é a direção do vetor gradiente. Seja α(x,y) o ângulo da direção do vetor na posição (x,y) e medido em relação ao eixo x, tem-se que no primeiro quadrante: ( ) ( ) (2.10) Considerando a Figura 10(a), onde os z’s são os valores dos níveis de cinzas, a Equação 2.8, por exemplo, pode ser escrita no ponto z5 de várias maneiras. Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 32 Figura 10 – Máscara de gradiente. (a) Região de tamanho 3x3 em uma imagem; (b) Operadores de Sobel; (c) Operadores de Prewitt. Fonte: Adaptado de GONZALEZ e WOODS, 2000. Utilizando os operadores de Sobel (Figura 10(b)), as derivadas baseadas nas máscaras são: ( ) ( ) (2.11) ( ) ( ) Então a magnitude do vetor gradiente (Equação 2.9) resulta em: ‖ ‖ |( ) ( )| |( ) ( )| (2.12) Para os operadores de Prewitt (Figura 10(c)) as derivadas baseadas nas máscaras são: ( ) ( ) (2.13) ( ) ( ) A magnitude do vetor gradiente (Equação 2.9) para esse caso resulta em: ‖ ‖ |( ) ( )| |( ) ( )| (2.14) O cálculo da magnitude do vetor gradiente é realizado em cada pixel da imagem para a obtenção da imagem de magnitude do gradiente. Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 33 2.3. Conceito de quina Segundo Mehrotra et al. (1990), uma quina é definida como um ponto de intersecção entre duas ou mais bordas. Uma importante concepção na área de visão computacional é a definição de “ponto de interesse”, o qual é considerado como sendo um ponto que tem posição bem definida e pode ser detectado de forma robusta. Isso significa que esse ponto pode ser uma quina, mas também pode ser, por exemplo, um ponto isolado com intensidade local máxima ou mínima, final de linhas, ou um ponto em uma curva onde a curvatura é localmente máxima. Dessa forma, quando se trata da detecção de quinas, um detector de quinas deve possuir as seguintes características (MEHROTRA et al., 1990): x Todos os pontos de quina devem ser detectados; x O detector de quina deve ser insensível aos ruídos; x Os pontos de falsas quinas não devem ser detectados; x O detector de quina deve fornecer o ângulo e a orientação da quina. A Figura 11 apresenta a linha do tempo com os anos de criação dos principais detectores de pontos de interesse, também conhecidos como operadores de interesse. Figura 11 – Linha do tempo da criação dos principais operadores de interesse. Fonte: Jazayeri e Fraser (2010). O primeiro operador de interesse foi desenvolvido por Hans P. Moravec em 1977 (MORAVEC, 1977), que introduziu também o conceito de ponto de interesse. Desde então o operado de Moravec vem sendo usado em várias aplicações na área de fotogrametria, incluindo a reconstrução 3D de objetos. No entanto, esse operador possui algumas limitações devido à sua baixa taxa de repetibilidade e sensibilidade a ruídos (JAZAYERI e FRASER, Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 34 2010). A baixa repetibilidade é consequência da natureza anisotrópica do algoritmo e qualquer rotação na imagem gera resultados diferentes na detecção de pontos de interesse. O operador de Moravec se baseia na medida de variância direcional, considerando quatro direções. A resposta desse operador considera os valores dos quadrados das diferenças dos tons de cinza nas direções 0o, 45o, 90o e 135o, como mostra a Figura 12 com janelas de dimensão 3×3 e conectividade 8. Figura 12 – Janelas 3×3 e os vetores nas direções 0o, 45o, 90o e 135o. Fonte: Galo e Tozzi (2002). Para o caso de janelas n×n, com n ímpar e origem do sistema de referência no centro da janela, as diferenças quadráticas são calculadas da seguinte maneira (GALO e TOZZI, 2002): ( ) ∑ ∑ ( ( ) ( )) (2.15) ( ) ∑ ∑ ( ( ) ( )) (2.16) ( ) ∑ ∑ ( ( ) ( )) (2.17) ( ) ∑ ∑ ( ( ) ( )) (2.18) onde N = (n-1)/2. O uso de janelas de dimensão ímpar não é restrição do operador de Moravec, mas se deve ao fato de se ter escolhido como origem o pixel central da janela, sendo a janela simétrica em relação a esse ponto. A partir dos resultados obtidos com as Equações (2.15) à (2.18), a resposta do operador de Moravec para o pixel (l, c) é calculada por: ( ) ( ) (2.19) Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 35 A Equação (2.19) é usada em cada ponto da imagem, e a partir dos resultados obtidos, os máximos locais de ( ) podem ser encontrados. Sendo assim, os pixels nos quais os valores de ( ) forem superiores a um limiar preestabelecido, corresponderão aos pontos de quina. 2.4. Conceito de sombra Em imagens aéreas de alta resolução é comum a presença de sombras quando a cena envolvida compreende regiões urbanas complexas, uma vez que o principal motivo para a ocorrência de sombras é o bloqueio direto da luz solar por meio de objetos altos como árvores e edifícios. Consequentemente, as superfícies afetadas são fracamente iluminadas e aparecem escuras nas imagens. Os efeitos negativos das sombras afetam a qualidade visual das imagens, causando perturbações severas na análise dessas imagens. Isso dificulta tarefas como a correspondência de imagens para a extração de elevações, detecção de alterações cadastrais, extração automática de edifícios e da malha viária (LI et al., 2004; MADHAVAN et al., 2004; MASSALABI et al., 2004). As principais características consideradas em análise de imagem para detectar regiões de sombras ou compensar seus efeitos negativos são (MASSALABI et al., 2004): x Sombras são regiões formadas por pixels com baixo valor de intensidade; x A forma de uma sombra é função da forma do objeto que a produz; x Uma ou mais extremidades de uma região sombreada são orientadas na direção do azimute solar (no instante de tomada da imagem); x A área de uma sombra depende da altura solar (no instante de tomada da imagem) e da altura do objeto que a produz; x A sombra não modifica a saturação da cor de um objeto; x Alguns elementos da textura de uma superfície são invariantes à sombra, ou seja, a textura da superfície praticamente não é afetada pelo sombreamento. Em geral, os métodos de detecção regiões de sombras podem ser classificados em dois grupos. No primeiro as regiões de sombras geralmente são obtidas utilizando-se dados de imagens combinados com outras fontes de informação. Em geral são utilizadas imagens combinadas com informações solares e de altura dos objetos. A estratégia adotada é semelhante em praticamente todos os métodos encontrados na literatura: as regiões de Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 36 sombras são determinadas através de uma simulação (predição) empregando dados de elevação, geralmente derivados de um Modelo Digital de Superfície (MDS), e dados de posição solar (azimute e altura), estimados para o instante do imageamento. Fazan e Dal Poz (2008) apresentam um método para a predição de regiões afetadas por sombras de edifícios em uma imagem aérea. Os contornos de telhado de edifícios obtidos a partir de um MDS e as informações sobre a posição solar no instante de aquisição da imagem são utilizados como base para a obtenção das regiões de sombra, que posteriormente são registradas na imagem aérea. O segundo grupo emprega dados de imagem em conjunto com técnicas convencionais de análise de imagem (limiarização, segmentação, etc.), explorando principalmente as propriedades radiométricas das sombras. No trabalho de Azevedo et al. (2015) a detecção das regiões de sombra em imagens é realizada através de operações de morfologia matemática. O primeiro passo consiste em realçar os padrões escuros da imagem através da transformação top-hat por fechamento de área. A Seção 2.7.2.1 apresenta mais detalhes sobre essa transformação. Como as sombras em imagens aéreas não obedecem um padrão específico, a operação de fechamento algébrico por área é utilizado para que não haja a necessidade de estabelecer uma forma particular para o elemento estruturante. O próximo passo consiste em realizar uma limiarização na imagem resultante da transformação top-hat, a fim de detectar as sombras realçadas. O histograma da imagem resultante da transformação top-hat possui um comportamento bimodal, ou seja, os alvos de interesse se destacam dos demais alvos da imagem. Dessa forma, o método de Otsu é utilziado para realizar a limiarização, pois ele consegue determinar automaticamente o limiar mais adequado de acordo com uma medida de separação de classes. O resultado é uma imagem binária onde os pixels de sombra possuem o valor 255 e os demais pixels possuem o valor 0. Por fim, a operação de abertura morfológica por área é aplicada na imagem binária para eliminar regiões de sombra pequenas, as quais podem ser consideradas ruídos. 2.5. Áreas de oclusão e mapa de visibilidade As imagens aéreas apresentam a projeção das feições do mundo real em um plano por meio da projeção via perspectiva central (Figura 13(a)). Nas projeções de perspectiva central os objetos com diferentes altitudes, mas com a mesma posição planimétrica, são projetados em diferentes posições na imagem, causando a não uniformidade da escala da Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 37 imagem. A Figura 13(a) ilustra um exemplo da projeção perspectiva central, onde os pontos A e B estão na mesma posição planimétrica, no entanto, são projetados em diferentes posições na imagem devido a suas diferenças de altitude (OLIVEIRA, 2016b). O deslocamento devido ao relevo presente na projeção de perspectiva central não ocorre na projeção ortogonal (Figura 13(b)), onde a escala é mantida e as feições são representadas em suas posições planimétricas corretas na imagem (MIKHAIL et al., 2001). Figura 13 – (a) Visada em perspectiva central. (b) Visada ortogonal. Fonte: Oliveira (2016b) Segundo Mikhail et al. (2001), a partir do processo denominado ortorretificação é possível transformar uma imagem, que está em projeção de perspectiva central, para a projeção ortogonal. Como resultado tem-se uma imagem onde todos os elementos são representados ortogonalmente em um plano de referência, assim como acontece em um mapa. Essa transformação é realizada por meio de um Modelo Digital de Superfície (MDS) e dos parâmetros de orientação interior e exterior do sensor utilizado na aquisição das imagens. O MDS é uma representação digital das informações de elevação da superfície, o qual inclue tanto o terreno quanto os elementos acima dele, como por exemplo, vegetação e feições antrópicas (WOLF et al., 2014). A utilização do MDS na ortorretificação tem como objetivo eliminar o deslocamento causado pela variação de altura dos objetos acima do terreno, permitindo assim a construção de uma nova imagem de forma ortogonal (NIELSEN, 2004). No entanto, como é apresentado por Nielsen (2004), Habib et al. (2007) e Oliveira (2013), o uso do MDS no processo de ortorretificação causa o efeito denominado duplo Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 38 mapeamento. A Figura 14 ilustra um caso de duplo mapeamento. Nesse exemplo, nota-se que cada um dos pixels a, b e c relaciona-se com mais de um groundel. O pixel a possui relação com os groundels A e D, o pixel b com B e E, o pixel c com C e F. Entretanto, como D, E e F pertencem ao terreno e não são visíveis ao CP, durante o processo de ortorretificação eles são preenchidos com os mesmos valores de tons de cinza de A, B e C, respectivamente. Isso gera uma ambiguidade, causando o efeito denominado duplo mapeamento. Figura 14 – Ilustração do efeito de duplo mapeamento. Fonte: Oliveira (2016b) A Figura 15(a) apresenta a imagem de um edifício elevado tomado em perspectiva central, onde ele parece estar “deitado” sobre a superfície do terreno. Após a ortorretificação o topo do edifício é projetado em sua posição correta, mas é criado o efeito de duplo mapeamento (Figura 15(b)). A Figura 15(c) apresenta em vermelho a região afetada pelo duplo mapeamento, essa região é denominada área de olcusão. Ao se aplicar um processo de detecção de oclusão, os tons de cinza da área duplicada podem ser convertidos para uma cor sólida, como mostra a Figura 15(c), gerando os mapas de visibilidade. Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 39 Figura 15 – (a) Imagem original em projeção perspectiva. (b) Efeito de duplo mapeamento. (c) Área de oclusão (região em vermelho). Fonte: Oliveira (2016b) Segundo Oliveira (2013), os métodos de detecção de oclusão mais utilizados são: método Z-buffer, foi o primeiro a ser utilizado para realizar a detecção de oclusão na geração de ortoimagens verdadeiras e é muito utilizado em computação gráfica (AMHAR et al., 1998); método baseado em ângulos proposto por Habib et al. (2007), o qual foi criado para evitar algumas imprecisões resultantes do método Z-buffer; método proposto por Volotão (2001), que utiliza a projeção das bordas de edifícios para identificar as células de oclusão. Oliveira (2016b) apresenta também um novo método de detecção das áreas de oclusão, denominado Surface-Gradient-Based Method (SGBM), o qual se baseia na análise de gradientes de altura negativos sobre um MDS. Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 40 2.6. Sistema de varredura a LASER aerotransportado (SVLA) Devido à grande necessidade de aquisição rápida e eficaz de dados de elevação, o uso da tecnologia de VLA tem sido foco de pesquisas desde a década de 1990, uma vez que esses sistemas permitem coletar, em um curto intervalo de tempo, uma alta densidade de pontos tridimensionais com precisão e acurácia equivalentes às técnicas tradicionais de levantamento in situ e Fotogrametria. Isso atribui aos dados de VLA uma variedade de aplicações como, reconstrução de superfície, extração de feições, modelagem urbana, mapeamento de zonas costeiras, entre outros. O funcionamento do sistema de VLA é baseado na utilização de pulsos LASER que são disparados em direção à superfície terrestre com o auxílio de um sistema de varredura. Esses pulsos interagem com a superfície e parte do sinal é refletido de volta em direção ao sensor. O sensor mede a intensidade do sinal retornado e usa o tempo decorrido entre a emissão e a recepção do pulso para calcular a dist ncia entre o sensor e a superfície amostrada. ssas medidas de dist ncias, em conjunto com as informaç es de ngulo de varredura, das posiç es e atitudes da plataforma e medidas de calibração entre todos os componentes do sistema de V A, permitem a determinação das coordenadas tridimensionais de diversos pontos no terreno, formando um conjunto de pontos denominado nuvem de pontos LASER. Os sistemas de VLA usados para mapeamento topográfico operam na região do infravermelho, que compreende o intervalo de comprimento de onda entre 1040 nm e 1060 nm (BALTSAVIAS, 1999). A qualidade e a acurácia dos dados podem ser afetadas por alguns fatores, tais como a superfície do material, altura de voo, integração GNSS (Global Navigation Satellite System)/INS (Inertial Navigation System), ângulo de observação, tipo de sensor utilizado, entre outros. Cada pulso LASER cobre uma área aproximadamente circular, determinada pelo IFOV (Instantaneous Field of View – campo de visada instantânea), que é a área da superfície irradiada pelo sensor ativo. Dentro desta área circular (footprint), o pulso pode gerar um retorno ou múltiplos retornos, caso o pulso LASER encontre um ou mais objetos (por exemplo, folhas e galhos de árvores ou extremidades de outras feições elevadas) antes de atingir a superfície do terreno (Figura 16). Em geral, o primeiro retorno apresenta a reflexão Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 41 de objetos mais elevados e o último retorno apresenta a reflexão da superfície do terreno (JENSEN, 2009). Figura 16 – Situações em que ocorrem múltiplos retornos. 2.6.1. Posição e orientação do feixe LASER A determinação das coordenadas dos pontos amostrados pela VLA é realizada combinando as medidas oriundas de cada componente que integra o sistema de VLA e os parâmetros que relacionam esses componentes. A unidade LASER é reponsável por realizar a medida da distância entre o emissor do pulso e um ponto na superfície atingida. As coordenadas desse ponto são matematicamente associadas a um sistema de referência global relacionando o ângulo de varredura do sistema de VLA, a posição (receptores GNSS) e a atitude (INS) da plataforma no instante da emissão do pulso LASER. Esse procedimento é realizado para todos os pontos amostrados durante o levantamento, obtendo-se assim um conjunto de dados tridimensionais da superfície. A Figura 17 apresenta os componentes de um sistema de VLA (BALTSAVIAS, 1999; PETRIE E TOTH, 2008): Unidade LASER contém o conjunto óptico de emissão e recepção do AS , detector de sinal, amplificador (necess rio para emissão da luz), sistema de contagem de Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 42 tempo e outros componentes eletr nicos. Após ser estimulado, o pulso AS é direcionado ao espelho de varredura, e o conjunto óptico orienta o pulso direção desejada. Ao chegar na superfície, o pulso interage com os objetos, retorna ao sistema de V A e é coletado por meio de um receptor. O sinal analógico captado é convertido para digital e filtrado para a eliminação de possíveis ruídos. or fim, esse sinal é utilizado para o c lculo da dist ncia entre o emissor do pulso AS e a superfície amostrada. INS: é composto por computadores, instrumentos eletrônicos de apoio e a IMU (Inertial Measurement Unit – Unidade de Medição Inercial). A IMU é respons vel pela mensuração das aceleraç es lineares e velocidades angulares da plataforma, o que possibilita a determinação das posiç es e da atitude da plataforma ao longo da trajetória do levantamento. ssas observaç es não são coletadas ao mesmo tempo em que são determinadas as posiç es espaciais estimadas pelo sistema G SS, criando a necessidade de um sincronismo entre os dois sistemas. Receptor GNSS respons vel pela coleta das observaç es necess rias para determinação das posiç es espaciais da plataforma em que se encontra o sistema de VLA, em intervalos de tempo pré-determinados. Na Figura 17 é considerado que o centro de fase da antena do receptor GNSS coincide com a origem do sistema de coordenadas do INS. Essa simplificação pode ser considerada desde que o vetor ligando a origem do sistema de coordenadas do INS e do receptor G S seja conhecido ou determinado a priori com alta precisão. staç es G SS de refer ncia: são utilizadas para realizar um pós-processamento e obter posiç es acuradas para a trajetória da plataforma. Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 43 Figura 17 – Componentes de um sistema de VLA e suas relações geométricas. Fonte: Adaptado de Habib et al. (2008) Segundo Habib et al. (2008), a posição de um ponto i no terreno ( ⃗⃗⃗⃗ ⃗) pode ser determina por meio da somatória de três vetores ( ⃗⃗⃗⃗ ⃗⃗⃗⃗ ⃗⃗⃗ ) e aplicando-se as rotações necessárias ( ): ⃗⃗⃗⃗ ⃗ ⃗⃗⃗⃗ ( ) ( ) ⃗⃗⃗⃗ ( ) ( ) [ ] (2.20) onde: ⃗⃗⃗⃗ ⃗ – coordenadas de um ponto genérico i no sistema de coordenadas terrestre; ⃗⃗⃗⃗ – coordenadas da origem do sistema de coordenadas do GNSS/INS, no sistema de coordenadas terrestre; ⃗⃗⃗⃗ – vetor translação entre a origem do sistema de coordenadas G SS I S e a origem do sistema de coordenadas da unidade AS . sse vetor translação também é conhecido como lever arm offset; – matriz de rotação entre o sistema de coordenadas G SS I S e o sistema de coordenadas terrestre; Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 44 – matriz de rotação entre o sistema de coordenadas G SS I S e o sistema de coordenadas da unidade AS . Os ngulos de rotação que comp em essa matriz são conhecidos como ngulos de desalinhamento (boresight angles); – matriz de rotação entre o sistema de coordenadas da unidade LASER e o sistema de coordenadas do raio AS . Os ngulos que comp em essa matriz são conhecidos como ngulos de varredura – medida de dist ncia entre a origem do sistema de coordenadas da unidade LASER e o ponto objeto i no terreno; t – instante de recepção do pulso AS . A posição ( ⃗⃗⃗⃗ ) e a orientação ( ) da plataforma durante o levantamento são determinadas pela integração do G SS e INS. Em virtude do problema de deriva, inerente ao INS, as posiç es e velocidades obtidas pelo sistema G SS podem ser usadas como medidas externas para atualizar a informação gerada pelo I S, melhorando sua estabilidade e precisão ao longo do tempo. Em contrapartida, quando o sinal de recepção do GNSS é interrompido e/ou a geometria dos satélites não est dentro dos limites aceit veis, o sistema de navegação inercial pode proporcionar as informaç es necessárias para a navegação (EL-SHEIMY e NIU, 2007). Durante o deslocamento as coordenadas obtidas pelo sistema GNSS e as observaç es do sistema inercial não necessariamente coincidem com os instantes de emissão e recepção dos pulsos LASER. Por essa razão é fundamental que o sincronismo entre estes sistemas (GNSS, INS e emissor do pulso LASER) seja estabelecido, para que todos os pontos amostrados possuam coordenadas e valores de atitudes relativos sua real posição no momento da coleta. Outro fator importante é que os elementos de translação ( ⃗⃗⃗⃗ ) e desalinhamento angular (∆𝜅, ∆𝜑, ∆𝜔) devem estar referenciados entre si, ou seja, eles devem ser determinados a priori por meio de um processo de calibração do sistema (HA I et al., 2008). 2.6.2. Filtragem e classificação dos dados de VLA Ao fim da aquisição, o sistema de VLA fornece um conjunto de pontos (nuvem de pontos) distribuídos irregularmente ao longo da linha de voo, e com coordenadas Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 45 tridimensionais em relação ao sistema de referência WGS 84 (World Geodetic System 1984) determinadas através do processo apresentado anteriormente (JENSEN, 2009). Dependendo da finalidade do trabalho, a nuvem de pontos originada pela varredura a LASER pode passar por um processo de tratamento a partir da filtragem e classificação dos pontos. Esse processo auxilia na diminuição de dados irrelevantes, como por exemplo, a eliminação da vegetação em estudos focados na extração de edifícios e topografia. Esse pré-processamento dos dados LASER consiste em separar o terreno de objetos altos através do procedimento de filtragem e, em seguida, analisar individualmente as regiões contendo os objetos altos a fim de classificá-las como vegetação ou edifício, segundo um algoritmo específico de classificação, podendo então, eliminar as regiões indesejadas (vegetação ou edifício). Vale ressaltar que a etapa de filtragem, para a separação do terreno e dos objetos altos, também é necessária e precede o processo de geração do MDT. A seguir são apresentados alguns conceitos sobre cada um desses processos. 2.6.2.1. Filtragem Segundo Sithole e Vosselman (2004), a filtragem de dados LASER é fundamentada na combinação de diferentes elementos. Alguns deles são apresentados a seguir. 1) Estrutura dos dados: a nuvem de pontos produzida pelo sistema de VLA é composta por pontos tridimensionais irregularmente espaçados. Alguns algoritmos de filtragem (AXELSSON, 1999; PFEIFER et al., 1998; SITHOLE, 2001; SOHN e DOWMAN, 2002; ROGGERO, 2001) trabalham diretamente com a nuvem de pontos LASER original. Entretanto, para aproveitar as ferramentas da área de processamento de imagens, alguns algoritmos de filtragem (BROVELLI et al., 2002; ELMQVIST, 2001; WACK e WIMMER, 2002) reamostram a nuvem de pontos original para uma grade regular de pontos, antes de realizar a filtragem. 2) Definição da vizinhança a ser testada: no que diz respeito à vizinhança utilizada pelo algoritmo de filtragem para classificar os pontos em terreno ou objeto, três possíveis formas de classificação podem ser listadas (SITHOLE e VOSSELMAN, 2004): x Ponto contra ponto: nesses algoritmos (SITHOLE, 2001; ROGGERO, 2001) dois pontos são comparados ao mesmo tempo. Uma função discriminante é resolvida com base na posição desses dois pontos e, se a saída da função discriminante estiver acima de certo limiar, então um dos pontos é assumido Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 46 como pertencente ao objeto. Em outras palavras, somente um ponto é classificado de cada vez. x Ponto contra pontos: nesses algoritmos (AXELSSON, 1999; SOHN e DOWMAN, 2002) os pontos vizinhos a um ponto de interesse são usados para resolver uma função discriminante. Com base na saída da função discriminante o ponto de interesse pode ser classificado. Apenas um ponto é classificado a cada vez. x Pontos contra pontos: nesses algoritmos (ELMQVIST, 2001; PFEIFER et al., 1999, BROVELLI et al., 2002; WACK e WIMMER, 2002) vários pontos são usados para resolver uma função discriminante, que é posteriormente utilizada para classificar os vários pontos. Mais de um ponto é classificado nesse método. 3) Medida de descontinuidade: as descontinuidades dos objetos (por exemplo, edifícios) em relação ao terreno são utilizadas pelos algoritmos de filtragem para separar os pontos do terreno e de objetos elevados. Alguns exemplos de medidas de descontinuidade comumente utilizadas são: diferença de altura, declividade, menor distância em relação às faces do TIN (Triangulate Irregular Network), e menor distância de pontos a uma superfície. 4) Princípios de filtragem: geralmente os algoritmos de filtragem baseiam-se em algum modelo geométrico do terreno em uma vizinhança local. Quatro princípios de filtragem são apresentados por Sithole e Vosselman (2004): x Baseado na declividade: nesses algoritmos (SITHOLE, 2001; ROGGERO, 2001) a declividade ou a diferença de altura entre dois pontos é medida. Se o desnível exceder certo limiar, então o ponto mais alto é assumido como pertencente a um objeto elevado. Neste caso é assumida a hipótese básica de que o terreno varia suavemente. x Baseado no bloco mínimo: como apresentado por Wack e Wimmer (2002), a função discriminante neste princípio visa delimitar uma região em 3D (se assemelharia a um paralelepípedo horizontal ou vertical), tendo por referência um segmento local de plano horizontal a partir do qual, acima ou abaixo, os pontos de terreno devem ser encontrados. x Baseado numa superfície: nesse caso é utilizada uma superfície que modela globalmente a superfície do terreno. Por exemplo, tendo em vista que o terreno é suave, Elmqvist (2001) utiliza contorno ativo para reter apenas pontos do Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 47 terreno. Já Axelsson (1999) refina progressivamente uma representação poliédrica grosseira do relevo. x Baseado em segmentação/agrupamento: estes métodos (BROVELLI et al., 2002) levam em conta que pontos que se agrupam e que estão acima de um outro agrupamento adjacente, pertencem a um objeto. 5) Mecanismo da filtragem: o mecanismo do processo de filtragem pode ser iterativo ou não iterativo. Os algoritmos não iterativos (SITHOLE, 2001; ROGGERO, 2001) realizam a filtragem numa única passagem, enquanto os iterativos (ELMQVIST, 2001; AXELSSON, 1999; BROVELLI et al., 2002; PFEIFER et al., 1998; SOHN e DOWMAN, 2002; WACK e WIMMER, 2002) classificam pontos através de múltiplas passagens. Em geral, os métodos não iterativos são computacionalmente mais atrativos devido a sua velocidade na execução. No entanto, em contrapartida, os iterativos são geralmente mais acurados, baseando-se na justificativa de que a cada passagem mais informações sobre a vizinhança de um ponto é incorporada e, assim, uma classificação mais confiável pode ser obtida. 6) Natureza da filtragem: os algoritmos de filtragem podem remover os pontos filtrados do conjunto de dados, ou recolocar os pontos filtrados no conjunto de dados. Os métodos que removem os pontos filtrados (SITHOLE, 2001; ROGGERO, 2001; AXELSSON, 1999; PFEIFER et al., 1998; SOHN e DOWMAN, 2002) normalmente operam sobre a nuvem de pontos LASER original, com os pontos irregularmente espaçados. Já os métodos de recolocação (ELMQVIST, 2001; BROVELLI et al., 2002; WACK e WIMMER, 2002) retornam os pontos filtrados para o conjunto de dados, mas com diferentes alturas, tendo por base alguma estratégia de interpolação a partir de pontos vizinhos. Esses métodos geralmente operam sobre malhas regularmente espaçadas. 2.6.2.2. Classificação Na etapa da classificação das regiões altas, a propriedade do pulso LASER de penetrar entre a vegetação é o elemento chave para a distinção entre vegetação e edifício (DAL POZ, 2013). Ao transpassar a vegetação, o pulso LASER pode atingir o terreno, gerando regiões com superfícies irregulares e rugosas. Já no caso dos edifícios, que possuem telhados planos, o retorno é de apenas um pulso, gerando superfícies regulares. Com base Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 48 nessa diferença, algoritmos podem ser criados para a classificação de vegetação e edifício, como pode ser visto em Axelsson (1999). Os edifícios geralmente possuem telhados planos, o que origina uma rede triangulada de pontos descrita por triângulos conectados com orientações similares em cada face do telhado. Dessa forma, desconsiderando irregularidades geométricas (rugosidade do próprio telhado e detalhes do telhado) e erros grosseiros, é possível modelar um telhado através de segmentos de planos justapostos. Considerando inicialmente o caso unidimensional, ao longo de linhas de varredura, pode-se modelar um perfil de telhado através de sucessivas retas conectadas (DAL POZ, 2013). Para cada segmento de reta, pode-se escrever: ( ) (2.21) onde z é a elevação e x a distância ao longo do perfil. Pela Equação 2.21 tem-se que: (2.22) A Equação 2.22 mostra que idealmente bastaria procurar pontos com derivada segunda nula para detectar os pontos pertencentes ao plano considerado (AXELSSON, 1999). Mas em problemas práticos tem-se que levar em consideração o conceito de planaridade, que está associado a dois fatores básicos: 1) o fato do próprio telhado não ser rigorosamente plano, apresentando uma pequena rugosidade; e 2) o fato dos pontos apresentarem incertezas decorrentes do processo de medida do sistema de VLA. Assim, a derivada segunda precisa ser considerada nula em relação a um certo limiar que reflita essas incertezas, ou seja, um desvio-padrão possível de ser estimado a partir de incertezas associadas com a rugosidade do telhado e da precisão associada a cada ponto da nuvem obtida por posicionamento LASER. Os limites dos planos também podem ser identificados com base no critério de curvatura. Os pontos pertencentes a essas descontinuidades devem apresentar (AXELSSON, 1999): (2.23) Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 49 Os princípios apresentados anteriormente podem ser estendidos para o caso bidimensional, em que a análise é realizada sobre uma representação poliédrica obtida por uma rede triangular de pontos. Neste caso considera-se que as normais às faces triangulares, que compõem o plano em estudo, devem ser aproximadamente paralelas. Assim, a detecção de pontos pertencente a uma região plana definida por pontos LASER estruturados numa rede triangulada, pode ser feita utilizando o seguinte procedimento (DAL POZ, 2013): 1) Para cada ponto P, identificar os triângulos que compartilham o vértice P. 2) Determinar a curvatura da superfície em P através da média (m) dos ângulos entre os vetores normais aos triângulos identificados anteriormente. 3) A curvatura em P pode ser considerada nula (plano) se m for menor que um dado limiar, estimado a partir do conceito de planaridade. Numa sub-região englobando pontos praticamente coplanares, os vetores normais aos triângulos dessa região serão praticamente paralelos, assim, quanto mais plana for a sub- região em análise, mais próximo de zero será m. 4) Por fim, os pontos com curvatura considerada nula podem ser agrupados (como num processo de crescimento de regiões) para formar a região de telhado. No caso de uma região de vegetação, a varredura a LASER pode gerar um perfil do tipo apresentado na Figura 18. Figura 18 – Perfil gerado pela conexão de pontos LASER sequenciais numa região de vegetação. Fonte: Adaptado de Dal Poz (2013). Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 50 Como pode ser notado, o perfil gerado pela conexão dos pontos sequenciais é bastante irregular. Sendo assim, ao longo desse perfil a segunda derivada normalmente é diferente de zero e varia randomicamente (AXELSSON, 1999). Essa variação não é só local, mas se estende por um longo segmento, o que a diferencia das descontinuidades locais típicas de telhado (bordas, cumeeiras e terreno). Na prática, uma região estruturada na forma de uma rede triangular apresentaria os vetores normais às faces dos triângulos variando bastante em direção e de forma randômica. Isso possibilita que pontos com essas características possam ser agrupados para formar regiões de vegetação. 2.7. Morfologia Matemática A Morfologia Matemática (MM) é definida como uma ferramenta para extrair informações relativas à geometria e à topologia de uma imagem através de um subconjunto que possui forma e tamanho pré-definidos, denominado elemento estruturante (EE). O subconjunto conhecido do EE é comparado com o conjunto desconhecido da imagem através de uma transformação, com o intuito de testar e quantificar de que maneira o EE está ou não contido na imagem (FACON, 1996). Sendo assim, a forma e o tamanho do EE devem ser estabelecidos de acordo com a geometria dos objetos de interesse na imagem. A Figura 19 apresenta três exemplos de EE, sendo a origem do EE representada pelo pixel preto. Figura 19 – Exemplos de elementos estruturantes. Fonte: Adaptado de Soille (2004). A linguagem morfológica se baseia na teoria dos conjuntos para representar os objetos presentes na imagem. Em uma imagem binária, por exemplo, o conjunto de todos os pixels brancos pode descrever completamente a imagem, uma vez que os demais pixels só podem ser pretos. Assim, cada elemento do conjunto é um vetor bidimensional (x, y) Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 51 contendo, por convenção, as coordenadas dos pixels brancos da imagem. Já as imagens em níveis de cinza podem ser representadas por conjuntos de vetores tridimensionais. Nesse caso, os vetores têm três elementos, sendo (x, y) as coordenadas do pixel e o terceiro elemento o seu nível de cinza (SOILLE, 2004). Os principais operadores morfológicos são a dilatação e erosão. A partir da combinação desses dois operadores, outras operações mais complexas podem ser realizadas, como por exemplo a abertura e o fechamento. A seguir serão apresentados os conceitos básicos desses operadores. 2.7.1. Erosão e Dilatação Para o caso de imagens binárias, a erosão (ε) do conjunto X pelo elemento estruturante B é dada por (SOILLE, 2004): ( ) * | + (2.24) onde é o elementro estruturante centrado no pixel . Nessas condições, o elemento estruturante B desliza sobre a imagem X, comparando cada pixel com a vizinhança de x. O resultado serão todos os pixels x quando B estiver inteiramente contido em X. Os efeitos da erosão em uma imagem binária são (Figura 20(b)): x Eliminação de objetos com tamanho inferior ao elemento estruturante; x Diminuição dos objetos presentes na imagem; x Aumento dos buracos no interior dos objetos; x Separação de objetos próximos. Para uma imagem h em níveis de cinza, a erosão de h pelo elemento estruturante B, quando a origem de B estiver em x, será o valor mínimo da imagem (h) contida na região definida por B, ou seja: , ( )-( ) ( ) (2.25) Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 52 Os efeitos da erosão em uma imagem em níveis de cinza são (Figura 20(a)): x Escurecimento da imagem; x Conexão de padrões escuros; x Alargamento e aumento de padrões escuros; x Redução e eliminação de padrões claros; x Separação de padrões claros. Figura 20 – Efeitos da erosão utilizando o EE octógono de tamanho 21 pixels. (a) Imagem binária. (b) Imagem em níveis de cinza. Imagem original Imagem erodida (a) Imagem original Imagem erodida (b) O processo de dilatação (δ) do conjunto X pelo elemento estruturante B é dado por (SOILLE, 2004): ( ) * | + (2.26) onde é o elementro estruturante centrado no pixel , e é um conjunto vazio. Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 53 Os efeitos da dilatação em imagens binárias são (Figura 21(a)): x Preenchimento de pequenos buracos; x Aumento do tamanho dos objetos; x Conexão de objetos próximos. Para imagens em níveis de cinza, a dilatação de uma imagem h por elemento estruturante B é dado por (SOILLE, 2004): , ( )-( ) ( ) (2.27) A dilatação em imagens em níveis de cinza causa os seguintes efeitos (Figura 21(b)): x Clareamento da imagem; x Conexão de padrões claros; x Alargamento e aumento de padrões claros; x Redução e eliminação de padrões escuros; x Separação de padrões escuros próximos. Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 54 Figura 21 – Efeitos da dilatação utilizando o EE octógono de tamanho 21 pixels. (a) Imagem binária. (b) Imagem em níveis de cinza. Imagem original Imagem dilatada (a) Imagem original Imagem dilatada (b) 2.7.2. Abertura e Fechamento A abertura (γ) de uma imagem h por um elemento estruturante B é definida pela erosão de h por B, seguida da dilatação usando o elemento estruturante transposto B’ (SOILLE, 2004): ( ) , ( )- (2.28) Os efeitos do processo de abertura em imagens em níveis de cinza são: x Conexão de padrões escuros próximos; x Conservação de padrões escuros distantes; x Eliminação de padrões claros menores que o elemento estruturante; x Separação de padrões claros próximos; x Diminuição dos detalhes da imagem. Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 55 O fechamento (ϕ) de uma imagem h por um elemento estruturante B é definido pela dilatação de h por B, seguida da erosão usando o elemento estruturante transposto B’ (SOILLE, 2004): ( ) , ( )- (2.29) Os efeitos do processo de fechamento em imagens em níveis de cinza são: x Conexão de padrões claros próximos; x Conservação de padrões claros distantes; x Eliminação de padrões escuros menores que o elemento estruturante; x Separação de padrões escuros próximos; x Diminução dos detalhes da imagem. 2.7.2.1. Transformação Top-hat A transformação Top-hat tem como objetivo detectar os elementos relevantes de uma imagem a partir do realce das regiões claras ou escuras. Para isso ela realiza operações aritméticas entre a imagem original e a imagem resultante do processo de abertura ou fechamento morfológico. O processo de abertura ou fechamento de uma imagem por um EE remove as feições que não se ajustam a esse EE. Dessa forma, as feições removidas podem ser recuperadas através da diferença aritmética com a imagem original. Existem dois tipos de transformação top-hat, por abertura ou por fechamento. Top-hat por abertura A transformação top-hat por abertura, denominada White Top-hat (WTH), é obtida a partir da diferença entre uma imagem h e sua abertura γ (SOILLE, 2004): ( ) ( ) (2.30) A Figura 22 apresenta o princípio da transformação top-hat por abertura. É possível notar que a diferença apresentada na equação acima permite ressaltar as informações dos maiores valores de intensidade da imagem, ou seja, permite ressaltar os picos das imagens. Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 56 Figura 22 – Exemplo da transformação top-hat por abertura em um sinal unidimensional. Fonte: Adaptado de Soille (2004) Top-hat por fechamento A transformação top-hat por fechamento, denominada Black Top-hat (BTH), consiste an difenreça entre o fechamento ϕ de uma imagem e a imagem original h (SOILLE, 2004): ( ) ( ) (2.31) A diferença apresentada na equação acima permite detectar os vales de uma imagem, ou seja, ressalta as informações dos menores valores de intensidade da imagem. A Figura 23 ilustra o princípio da transformação top-hat por fechamento, onde todas as estruturas escuras relevantes do fundo foram suprimidas e após a diferença com a imagem original por meio do BTH, elas foram recuperadas. Figura 23 – Exemplo da transformação top-hat por fechamento em um sinal unidimensional. Fonte: Adaptado de Soille (2004) Extração automática de contornos de telhados de edifício no espaço-objeto integrando um estéreo par de imagens aéreas de alta resolução e modelos 3D de telhado YWATA, M. S. Y. 57 2.8. Modelo de contorno ativo: Snakes 2.8.1. Modelos deformáveis generalizados Um modelo deformável generalizado pode ser considerado como sendo uma superfície n-dimensional constituída fisicamente por um material elástico que se deforma sob a ação de forças internas e externas atuando sobre ele