SciELO - Scientific Electronic Library Online

 
vol.21 número3INTEGRAÇÃO DE DADOS DE IMAGENS ORBITAIS DE ALTA RESOLUÇÃO E ALS PARA DETECÇÃO SEMI-AUTOMÁTICA DE EDIFICAÇÕES EM ÁREAS URBANASCOMPATIBILIZAÇÃO DE REFERENCIAIS DE COORDENADAS E VELOCIDADES COM ESTIMATIVA DE PRECISÃO índice de autoresíndice de assuntospesquisa de artigos
Home Pagelista alfabética de periódicos  

Serviços Personalizados

Journal

Artigo

Indicadores

Links relacionados

Compartilhar


Boletim de Ciências Geodésicas

versão impressa ISSN 1413-4853versão On-line ISSN 1982-2170

Bol. Ciênc. Geod. vol.21 no.3 Curitiba jul./set. 2015

http://dx.doi.org/10.1590/S1982-21702015000300032 

Artigos

UM MÉTODO AUTOMÁTICO PARA REGISTRO DE DADOS LASER SCANNING TERRESTRE USANDO SUPERFÍCIES PLANAS

An automatic method for registration of terrestrial laser scanning data using planar surfaces

Nadisson Luis Pavan1 

Daniel Rodrigues dos Santos1 

1Universidade Federal do Paraná Setor de Ciências da Terra Programa de Pós-Graduação em Ciências Geodésicas CEP 81531-Curitiba/PR-Brasil nadissonluisp@gmail.com; danielsantos@ufpr.br

Resumo:

Neste trabalho é apresentado um método automático para registro de pares de nuvens de pontos derivados do sistema LASER scanning terrestre usando superfícies planas. Primeiramente, as superfícies planas são detectadas através do algoritmo RANSAC e, posteriormente, ajustadas pelo método dos Mínimos Quadrados Totais. Em seguida é proposto um algoritmo para estabelecimento automático de correspondências baseado em análises geométricas dos vetores normais aos planos. Finalmente, um modelo matemático baseado em abordagens plano-a-plano é empregado para determinar os parâmetros de transformação. Um experimento foi realizado para avaliar a viabilidade e potencialidade do método proposto. Os resultados obtidos mostraram que esta abordagem é promissora e determina medidas com precisão na ordem do centímetro

Palavras-chave: Registro de pares de nuvem de pontos 3D; LASER scanning Terrestre; Quatérnios, correspondência plano-a-plano

Abstract:

In this work, an automatic method for registration of terrestrial laser scanning data using planar surfaces is presented. Firstly, planar surfaces are detected using the RANSAC algorithm and their normal vectors are adjusted by the Total Least Squares method. After, an automatic algorithm to establish correspondences using the normal vectors and geometric planar analysis is proposed. Finally, the transformation parameters are determined using a plane-to-plane based model. The proposed algorithm solves the transformation parameters without the initial estimates and it applies quaternions to describe spatial rotations, making it faster and better than conventional approaches. Experiment was conducted to verify the feasibility and the effectiveness of the proposed method. The obtained results showed that this approach is promising and can computed positions despite measure errors of a few centimeters

Keywords: Pairwise registration of 3D point clouds; Terrestrial laser scanning; Quaternions; matching plane-to-plane

1. Introdução

Atualmente, sistemas LASER (Light Amplification by Stimulated Emission of Radiation) scanning terrestre têm sido amplamente empregados na reconstrução 3D de ambientes internos e externos, propiciando uma nuvem de pontos 3D de forma rápida, precisa e segura. Para garantir o perfilamento completo da cena usando o sistema LASER scanning terrestre (LST) é necessário instalar o equipamento em diferentes pontos de vista propiciando, desta forma, um conjunto de nuvens de pontos com certa porcentagem de sobreposição.

O desenvolvimento de algoritmos rápidos e robustos para determinação dos parâmetros de transformação entre pares de nuvens de pontos 3D tem marcado uma nova tendência no processamento de dados LST. Geralmente, são encontrados desalinhamentos angulares e deslocamentos lineares entre os pares de nuvens de pontos 3D obtidos no processo de perfilamento. Sendo assim, é essencial que o erro de transformação de corpo rígido entre os pares de nuvens de pontos seja minimizado. A tarefa que consiste em minimizar esses erros é conhecida como registro de pares de nuvens de pontos tridimensionais (3D).

O registro de pares de nuvem de pontos 3D é alvo de intensa investigação científica e pode ser dividido, basicamente, no problema de correspondência e de otimização. O problema de correspondência consiste em encontrar as primitivas na nuvem de pontos de referência com maior similaridade existente na nuvem de pontos de pesquisa. Já o problema de otimização consiste em definir uma transformação que minimiza a soma dos quadrados das distâncias entre dois conjuntos de dados LST. A hipótese conhecida é que se deve alinhar perfeitamente a nuvem de pontos de pesquisa em relação à nuvem de pontos de referência. Desta forma, a transformação escolhida deve permitir a determinação de parâmetros de transformação (parâmetros de rotações e translações) entre os pares de nuvens de pontos 3D. Vale ressaltar que diversas abordagens propostas foram baseadas no emprego da transformação de corpo rígido 3D.

Usualmente, os métodos de registro de pares de nuvem de pontos 3D são divididos em duas categorias, isto é, os métodos analíticos (Horn, 1987; Walker et al., 1991) e os métodos iterativos (Besl e McKay, 1992; Chen e Medioni, 1992). Basicamente, os métodos iterativos exigem valores iniciais aproximados que, frequentemente, introduz erros na estimativa dos parâmetros. Dentro desta categoria, o algoritmo conhecido como ICP (Iterative Closest Point) é o mais popular e foi originalmente desenvolvido por Besl e McKay (1992). Inicialmente, o algoritmo ICP estabelece pseudo correspondências ponto-a-ponto entre o par de nuvens de pontos 3D. O algoritmo admite que os pontos mais próximos sejam correspondentes. Em seguida, um problema de otimização deve ser solucionado. Iterativamente, a função minimiza a distância euclidiana entre os pontos e calcula os parâmetros de transformação. Neste aspecto, o algoritmo é de alto custo computacional implicando em uma das maiores desvantagens da técnica supracitada. De acordo com Rabbani et al. (2007) o ICP, na sua forma original, não prevê qualquer medida sobre a precisão e confiabilidade dos parâmetros estimados. Outra desvantagem do ICP é a necessidade de valores iniciais aproximados, uma vez que o algoritmo tem a tendência de convergir para mínimos locais, caso o desalinhamento inicial seja muito grande (Rabbani et al., 2007). Já as soluções analíticas não necessitam de valores aproximados e evitam linearizações passando a ser alvo de atenção da comunidade fotogramétrica. Dentre os métodos propostos o mais importante é a abordagem apresentada por Horn (1987). O autor empregou quatérnios como forma de representação da matriz de rotação para resolver o problema de registro de pares de nuvens de pontos 3D. Essa solução minimiza a soma dos quadrados da diferença entre as rotações calculadas e as observadas.

Os quatérnios têm um grande potencial de aplicação em Ciências Geodésicas (Galo e Tozzi, 2001). Shen et al. (2006) utilizaram quatérnios em transformações entre dois sistemas referenciais geodésicos. A principal vantagem dessa representação é que a linearização do modelo matemático não é necessária, dependendo de como é dada a solução. Os quatérnios também foram empregados no trabalho de Kim e Golnaraghi (2004) para analisar a orientação de imagens em tempo real, baseado em sinais provenientes de uma unidade de medida inercial de baixo custo. Além disso, os quatérnios preservam a ortogonalidade dos eixos do sistema referencial a serem transformados.

Dirigido pelo o que foi apresentado anteriormente, o presente trabalho tem como proposta apresentar um método automático de registro de pares de nuvens de pontos 3D derivados do sistema LST. O método proposto é baseado na relação plano-a-plano de um conjunto de primitivas detectadas nos pares de nuvens de pontos. A característica mais importante do método proposto é sua simplicidade computacional e a robustez da abordagem de correspondência, além do emprego de quatérnios para estimativa dos parâmetros de transformação entre os pares de nuvem de pontos. A principal característica do método de correspondência proposto é a aplicação de uma hierarquia de restrições que descarta o elevado número de falsas correspondências. De forma geral, as principais vantagens do método proposto são a simplicidade e robustez do algoritmo desenvolvido para estabelecer a correspondência entre os planos e o emprego de quatérnios para solucionar o problema de registro dos pares de nuvem de pontos 3D.

Este texto está organizado em cinco seções. A seção 2 apresenta brevemente os trabalhos relacionados apontando os métodos iterativos e analíticos. A seção 3 introduz o modelo proposto para o estabelecimento de correspondências e a inclusão de quatérnios no problema de otimização. As etapas para solucionar os problemas supracitados são detalhadamente descritas. Na seção 4 é apresentada uma análise comparativa do algoritmo implementado baseada em quatro posições obtidas com o sistema LST e uma estação total. Finalmente, as conclusões e recomendações para trabalhos futuros são apresentados na seção 5.

2. Trabalhos Relacionados

Como descrito anteriormente, o problema de registro de pares de nuvem de pontos 3D pode ser dividido em duas categorias, isto é, métodos iterativos e métodos analíticos. Enquadrado na categoria de métodos iterativos, o ICP sofreu diversas variações desde que foi proposto, tendo como principal característica evitar o estabelecimento de relações ponto-a-ponto. Chen e Medioni (1992) investigaram um algoritmo de registro de pares de nuvem de pontos 3D com base na minimização da soma dos quadrados das distâncias entre pontos e suas superfícies correspondentes. As correspondências ponto-a-plano podem levar a um modelo de estimativa mais robusta que converge rapidamente. Este algoritmo procura o ponto de intersecção entre a superfície de referência e do vetor normal do ponto de origem projetado ortogonalmente. Para encontrar o ponto de origem, os autores usaram uma técnica de pesquisa baseada no método de Newton-Raphson em um sistema de coordenadas ortogonal. Vale ressaltar que este método também não converge sem valores iniciais aproximados.

Fontanneli et al. (2007) resolveram o problema de localização empregando o algoritmo RANSAC (RANdom Sampling Consensus) combinado com o filtro de Kalman estendido para acompanhar a trajetória de um sistema SLT móvel. A principal desvantagem desse método é o alto custo computacional. Bae e Lichti (2008) propuseram dois métodos de registro de nuvem de pontos 3D, conhecidos como ICP modificado e ICP modificado/RANSAC. Os pontos correspondentes são obtidos através da variação de curvatura, do vetor normal da superfície estimada e do ângulo de variação da normalidade prevista. Embora este método de registro de nuvem de pontos 3D seja mais eficiente e preciso na busca de pontos correspondentes, em relação ao que se encontra no estado da arte, também exige valores iniciais aproximados. Durgham et al. (2013)apresentaram um método para o registro de nuvem de pontos através da estratégia de correspondência automática que emprega linhas retas. O processo de registro simula o método RANSAC para determinar os parâmetros de transformação, enquanto o algoritmo ICPP (Iterative Closest Projected Point) é utilizado para determinar a solução mais provável dos parâmetros de transformação dentre as diversas soluções possíveis.

Em relação às abordagens analíticas, a melhor solução possível é obtida em uma única etapa, ou seja, sem a necessidade de iterações e de valores iniciais aproximados. Arun et al. (1987) empregou a decomposição em valores singulares (Singular Valor Decomposition-SVD), e os mínimos quadrados para estimar os parâmetros de rotação. Já Horn (1987) utilizou quatérnios unitários como representação da matriz de rotação para solucionar, de forma analítica, a orientação absoluta de imagens. Vale ressaltar que os métodos analíticos supracitados são divergentes apenas na forma como o problema é formulado.

Outro aspecto importante a ser tratado neste trabalho está relacionado ao tipo de primitiva a ser empregado no processo de registro de pares de nuvem de pontos 3D. Devido à natureza discreta e à densidade de amostragem finita das nuvens de pontos 3D obtidas no processo de varredura do sistema SLT é altamente propenso a incertezas encontrar pontos correspondentes no par de nuvens de pontos (Rabbani et al., 2007). Primitivas lineares e superfícies planas têm sido amplamente utilizadas, uma vez que não exigem pré-sinalização de alvos (Durgham et al., 2013). Tais primitivas podem ser facilmente extraídas e identificadas em áreas urbanas, especialmente em ambientes antrópicos, além de não exigirem pré-sinalização de alvos e podem propiciar parâmetros de transformação entre pares de nuvens de pontos com alta precisão. Jaw e Chuang (2008) investigaram o emprego de pontos, linhas e planos no processo de registro de dados LST. Esta abordagem oferece um alto grau de flexibilidade e propicia melhor precisão. Park e Subbarao (2003) introduziram uma nova técnica de registro de nuvem de pontos 3D combinando a precisão do método baseado em abordagens ponto-a-plano e a velocidade da técnica baseada na projeção do ponto. Nesta técnica, o vetor tangente do ponto de referência é utilizado para encontrar a intersecção da linha de projeção do ponto de referência com o vetor normal do ponto correspondente ao ponto de referência. Para encontrar a intersecção das primitivas supracitadas, o ponto de referência é projetado para a superfície correspondente e reprojetada para o vetor normal que parte do ponto de referência. Rabbani et al. (2007) propuseram dois métodos para registro de pares de nuvem de pontos obtidos em cenas industriais. No primeiro método, os pontos são segmentados como objetos planos, esferas ou cilindros. Os parâmetros de transformação são calculados através do método de ajustamento por mínimos quadrados (MMQ). O segundo método, imita a ideia do ajustamento por feixes de raio de luz perspectivo empregado em Fotogrametria convencional. Estas abordagens combinam o registro e a modelagem da nuvem de pontos e, assim, evitam o acúmulo de erros. Segal et al. (2009) combinaram a versão original do algoritmo ICP e uma de suas variações, baseada em abordagem ponto-a-plano, em uma única estrutura probabilística. Esta abordagem mostrou ser mais robusta no estabelecimento de falsos positivos. Entretanto, os métodos de registro de pares de nuvem de pontos 3D supracitados são lentos e iterativos, além de exigirem uma estimativa aproximada da transformação, em que a convergência desses algoritmos depende da escolha dos parâmetros iniciais.

Brenner e Dold (2007) propuseram o uso de planos correspondentes para obtenção de valores iniciais. Primeiramente, o algoritmo extrai superfícies planas nas nuvens de pontos de referência e de pesquisa. Em seguida, seleciona três do número total de planos extraídos para calcular os parâmetros de transformação entre os pares de nuvens de pontos 3D. Neste trabalho, os autores também investigaram a complexidade combinatória da busca de planos correspondentes. Khoshelham e Gorte (2009) apresentaram um método para registro automático de pares de nuvens de pontos 3D. Os dados de referência eram extraídos de mapas cadastrais, enquanto superfícies planas são detectadas na nuvem de pontos derivada do sistema LST. O método é baseado na suposição de que as edificações são objetos poliédricos. Este método foi testado com um conjunto de dados simulados e reais. Os resultados mostraram um desempenho robusto na presença de planos supérfluos e ausentes. Khoshelham (2010) apresentou um método para localização automática de sensores SLT em ambientes internos. O autor empregou uma estratégia que inicia o processo de correspondência com um mínimo de três ou quatro pares de planos e um processo recursivo para estimativa dos parâmetros é iniciado. À medida que outros planos são detectados e suas correspondências são estabelecidas dentro da região de sobreposição entre os pares de nuvem de pontos, os parâmetros são recalculados e atualizados. Para calcular estas estimativas, o modelo linear de ajustamento por mínimos quadrados foi utilizado. Desta forma, não são necessários aproximações iniciais, propiciando maior eficiência e rapidez no estabelecimento de correspondências. Pathak et al. (2010) abordaram o problema de registro de pares de nuvem de pontos baseado em abordagens plano-a-plano. Esta abordagem não só é mais robusta, quando comparada com abordagens ponto-a-ponto, mas também mais rápida. Wang e Sohn (2010) apresentaram um método para co-registro automático de pontos 3D e plantas cadastrais. Uma limitação óbvia do método proposto é que ele pode produzir resultados equivocados quando duas ou mais fachadas de edificações possuem as mesmas estruturas. Taguchi et al. (2013) apresentou um algoritmo para mapeamento e localização simultâneo (Simultaneous Localization and Mapping - SLAM) de um sensor 3D (Kinect) em tempo real. Neste algoritmo os pares de nuvens de pontos 3D são registrados utilizando pontos e planos como primitivas. O algoritmo RANSAC é empregado como restrição geométrica entre pontos e planos para resolver o problema de correspondência. Esta abordagem mista também permite a determinação dos parâmetros de forma mais rápida e precisa do que usando apenas pontos. Vale ressaltar que todos os testes realizados para validação dos métodos supracitados foram realizados em ambientes internos, cuja presença de diversos planos é facilmente garantida.

3. Método

Basicamente, o método proposto neste trabalho é dividido em quatro etapas, a saber:

  • Superfícies planas são automaticamente extraídas nos pares de nuvem de pontos 3D usando o algoritmo RANSAC. Posteriormente, essas superfícies são ajustadas através do método dos mínimos quadrados totais;

  • A correspondência entre os planos é automaticamente estabelecida através de uma análise geométrica usando os atributos invariantes dos respectivos planos;

  • Pares de planos correspondentes são usados para estimativa dos parâmetros de transformação entre os pares de nuvens de pontos 3D;

  • Finalmente, a localização do sensor é calculada

3.1 Extração e segmentação automática dos planos

A fim de detectar automaticamente características planas existentes a partir de dados do sistema LST é utilizado o algoritmo RANSAC, proposto por Fischler e Bolles (1981). O algoritmo RANSAC empregado segue os seguintes passos (Rusu e Cousins, 2011):

  1. O algoritmo seleciona, aleatoriamente, três pontos não colineares nos dados do sistema LST;

  2. Através desses três pontos é calculado os parâmetros do plano definido geometricamente por esses três pontos;

  3. Conta o número total de inliers desse plano. O critério inlier é baseado em um limiar de distância máxima pré-estabelecida de cada ponto para a superfície plana hipotética;

  4. Se a distância for menor que o limiar pré-estabelecida, o ponto é considerado pertencente ao plano hipotético e agrupado em um conjunto de dados. Os passos (1), (2) e (3) são repetidas para até certo número de iterações também pré-estabelecidas. O conjunto com o maior número de pontos (ou seja, inliers) é classificado como conjunto de pontos contidos em um objeto plano encontrado;

  5. Os pontos contidos no objeto plano são removidos do conjunto destes dados e o algoritmo retorna ao passo (1) para detectar outra superfície plana

Após a segmentação dos planos, seus parâmetros são calculados usando o método dos mínimos quadrados totais proposto por Golub et al. (1980). Inicialmente, cada conjunto de pontos contidos em uma superfície plana é organizado em uma matriz da forma como segue:

sendo, xi, yi e zi as coordenadas do ponto pi para 1≤ in1 i n e n o número de pontos detectados no plano. A Figura 1 mostra um exemplo de superfícies planas detectadas pelo RANSAC.

Figura 1: Extração de superfícies planas. (A) Nuvem de pontos derivados de dados do sistema SLT. (B) Planos detectados pelo algoritmo RANSAC 

O método dos mínimos quadrados totais se baseia na decomposição de valores singulares (SVD). Deste modo, os parâmetros do plano a serem determinados correspondem ao vetor singular à direita da matriz, associado ao menor valor singular.

3.2 Descrição da técnica de registro

Geralmente, um plano π no espaço euclidiano é definido pelo vetor normal unitário;e a distância perpendicular d entre a origem do sistema de coordenadas e o plano (Steinbruch e Winterle, 2006). A condição de que um ponto se encontra em um plano π é expressa pela seguinte equação:

Portanto a, b, c e d são os parâmetros do plano π. Seja ;onde x', y' e z' são as coordenadas do ponto p no sistema de coordenadas da outra nuvem de pontos. Na ausência de erros sistemáticos, a transformação de corpo rígido do ponto p para o ponto p' é dada por:

sendo, M a matriz de rotação no espaço euclidiano tridimensional e t o vetor translação.

Do mesmo modo, seja ;o vetor normal unitário do plano π, onde a', b' e c', são as coordenadas do vetor u' no sistema de coordenadas da nuvem de pontos de pesquisa. Logo, a rotação do vetor u para o vetor u' é dada por:

Como o ponto p' se pertence no plano π, a seguinte equação é dada por:

sendo, d' a distância perpendicular entre o plano e a origem do sistema de coordenadas da nuvem de pontos de pesquisa. Substituindo a Equação 3 na Equação 5 tem-se:

Seguindo as propriedades de transposição de matrizes, e substituindo a Equação 2 na Equação 6 é possível reescrever a Equação 6 da seguinte forma:

Agora, é possível calcular a translação t e a rotação M através das Equações 7 e 4, respectivamente. Como em Brenner e Dold (2007) a determinação da transformação é feita em duas etapas. A primeira etapa calcula a translação, enquanto a segunda etapa calcula a rotação. A translação t é calculada pelo MMQ. Como a Equação 7 é linear em relação aos parâmetros de translação t e observações (d - d'), tem-se uma equação linear para cada par de planos correspondentes. Este conjunto de equações lineares pode ser expresso na forma matricial:

X o vetor de parâmetros de translação t , V o vetor de resíduos, L o vetor das observações, A a matriz formada pela Jacobiana da Equação 7 em relação aos parâmetros, em que cada linha da matriz A é composta pelos coeficientes do vetor normal u' para cada par de planos correspondentes (Dalmolin, 2004). A solução X da Equação 8 é dada por:

Sabendo que o modelo de estimativa é linear não há necessidade de valores iniciais aproximados para a determinação dos parâmetros de translação t. Agora o cálculo da matriz de rotação M é feito utilizando quatérnios como proposto por Horn (1987). A seguir são apresentados os quatérnios e o operador de rotação que será utilizado na seção 3.2.2 para obter a rotação entre os planos correspondentes.

3.2.1 Quatérnios e o Operador de Rotação

O quatérnio unitário pode ser escrito em termos do número real q0 e um vetor q = (q1,q2,q3), como segue:

Deste modo, o operador ρq : R3 → R3 é definido em Gomes e Velho (2008) pelo quatérnio unitário q, a saber:

sendo, q* o conjugado do quatérnio q definido como q* = q0 - q, e (.) representa o produto interno euclidiano (ou produto escalar) entre vetores do espaço vetorial euclidiano R3.

Gomes e Velho (2008) ainda demonstraram que o operador ρq é um operador linear e também um operador ortogonal, ou seja, não altera o comprimento do vetor v. Assim, é admissível determinar a matriz associada ao operador de rotação que é dada por:

sendo, M a matriz de rotação que é ortogonal e com det (M) = 1.

Gomes e Velho (2008) também mostraram que o operador ρq é uma rotação em torno do vetor q, ou seja, M é a matriz de rotação. Uma importante propriedade dos quatérnios p, r e q é apresentada por Horn (1987), como segue:

sendo, < > o produto interno de dois vetores genéricos. Maiores detalhes sobre os quatérnios em Horn (1987), Galo e Tozzi (2001), Gomes e Velho (2008).

3.2.2 Rotação entre os planos correspondentes

Sejam u1', u2' ,...,un' vetores normais dos planos da nuvem de pontos de referência que correspondem aos respectivos vetores normais u1, u2, ..., un da nuvem de pontos de pesquisa, sendo essas correspondências pré-determinadas. Então, pode-se encontrar uma matriz de rotação M ortogonal com det (M) = 1, como solução para a seguinte minimização:

sendo que o vetor corresponde ao respectivo plano com o vetor normal com 1≤ i≤n.

O somatório acima pode ser reescrito nos termos utilizando as propriedades de produto escalar e norma, como segue:

Como todos os vetores normais u'i e ui sendo 1≤ i≤ n são unitários, então, tem-se:

Agora as rotações são representadas utilizando os quatérnios unitários. Este quatérnio q maximiza a equação, da seguinte forma:

Em seguida, utilizando algumas propriedades dos quatérnios sobre o produto interno (Horn, 1987), é possível obter a seguinte igualdade:

para 1≤ i≤ n.

Segundo Gomes e Velho (2008) o produto à direita qui e o produto à esquerda u'i q podem ser reescritos através da multiplicação de uma matriz 4×4 por um vetor coluna 4×1, ou seja

sendo Q = T as componentes do quatérnio q= q 0 + q 1 i + q 2 j + q 3 k para 1≤ i≤ n.

Também pode-se escrever o produto escalar ρq(ui) . u' . como produtos de matriz, e assim para 1≤ i≤ n tem-se:

Logo o somatório acima equivalente ao produto das matrizes:

sendo.

É fácil verificar que a matriz W é uma matriz 4×4 simétrica. Então, W tem os autovalores reais, λ1, λ2, λ3, λ4 onde λ1 ≥ λ2 ≥ λ3 ≥ λ4 com os correspondentes auto-vetores w1, w2, w3, w4, unitários e ortogonais entre eles. Desta forma, pode-se aplicar o teorema apresentado por Anton e Rorres (2001). Portanto, o quatérnio q que maximiza o somatório da Equação 22 é o auto-vetor Q, que corresponde ao maior auto-valor da matriz W. Logo, obtêm-se a rotação para o registro dos dados.

3.3 Método de correspondência de planos

Para estimar os parâmetros de transformação é necessário encontrar pares de planos correspondentes nas nuvens de pontos de referência e de pesquisa. Neste trabalho, o maior problema do registro de nuvens de pontos é estabelecer o maior número de planos correspondentes possíveis. Mas, devido à falta de informação causada pelos efeitos de oclusão e mudança de ponto de vista, alguns planos na nuvem de pontos de referência podem não apresentar correspondentes na nuvem de pontos de pesquisa. Por isso, o algoritmo de correspondência deve ser robusto às situações supracitadas.

Basicamente, como o número de planos é significativamente menor que o número de pontos encontrados na nuvem de pontos 3D, o estabelecimento da correspondência entre os planos pode ser iniciado a partir de todas as combinações possíveis de planos. No entanto, a busca por planos correspondentes é um processo crucial e pode ser de alto custo computacional, uma vez que o espaço de busca deve aumentar conforme o número de planos aumenta (Khoshelham e Gorte, 2009).

Para reduzir o custo computacional envolvido no problema de correspondência Brenner e Dold (2007) sugeriram uma hierarquia de restrições que descarta falsos positivos. A primeira restrição procura classificar os planos detectados em planos de solo e planos da fachada de edificações (parede). Essa restrição é baseada na suposição de que as fachadas e o solo, nos dados LST, têm pequenos desvios de planos verticais e horizontais, respectivamente. Esta é uma suposição real, uma vez que o equipamento SLT é instalado em um tripé e devidamente nivelado utilizando um nível de bolha (Khoshelham, 2010). Para essa tarefa é analisada a terceira componente do vetor normal unitário. Se esta componente estiver próxima do valor unitário (1), ou seja, dentro de um intervalo de 0,9986 até 1,00, este plano é classificado como plano do solo. Por outro lado, se esta componente estiver próxima de 0 (zero), ou seja, dentro de um intervalo de 0,00 até 0,006 este plano é classificado como plano da fachada. Caso esta componente for negativa, todos os parâmetros do plano são multiplicados por -1, para não descartar uma correspondência verdadeira. Esta operação de reflexão é importante, tendo em vista a natureza irregular do perfilhamento a LASER.

Outra informação geométrica invariante, identificada na transformação de corpo rígido de planos, é o ângulo entre dois vetores normais aos respectivos planos. Por exemplo, dado um par de u1 e u2 (da nuvem de pontos de referência) correspondentes ao par de vetores de u1' e u'2 (da nuvem de pontos de pesquisa), os ângulos entre os vetores u1 e u2 devem ser o mesmo que os ângulos entre os vetores u1' e u'2. Devido aos erros de medida presentes nos dados LST, esses ângulos são ligeiramente diferentes; entretanto, são ainda muito próximos. Portanto, o segundo passo deste algoritmo de correspondência entre planos se constitui da seguinte diferença dada pela Equação 23, a saber:

Na Equação 23, se θ estiver próximo de 0°, ou seja, dentro de um intervalo de -3° a 3°, os planos são considerados correspondentes. Caso contrário, a hipótese de correspondência é totalmente rejeitada.

O terceiro passo é aplicado quando o número de pares de planos correspondentes é maior do que três. A razão disso é explicada pelo fato de que a solução da Equação 7(que resulta na translação do registro do par de nuvens de pontos) necessita de no mínimo três pares de planos correspondentes. Uma ressalva importante é que estes planos devem se interceptar em um ponto. Este passo se inicia com a estimação dos parâmetros de transformação, calculados através de cada quatro pares de combinações de planos correspondentes. O cálculo das estimativas é feito através dos métodos apresentados nas subseções 3.2.1. e 3.2.3. A combinação considerada correta é aquela que contém o menor erro de verificação. O erro de verificação utiliza o centroide dos pontos pertencentes ao plano com o vetor normal u i, como segue:

sendo, 1≤ i≤ 4.

No quarto passo, os parâmetros de transformação estimados no terceiro passo são utilizados para encontrar mais pares de planos correspondentes. Para isto, é verificado o erro de cada combinação, como segue.

Se erroplano< erro, então o par de planos é considerado correspondente. Após todas as verificações, os planos correspondentes potenciais são utilizados para refinar os parâmetros de transformação.

4. Experimentos e Análise dos Resultados

O desempenho do método proposto foi testado usando um conjunto de dados reais derivados do sistema LST Leica HDS 3000. A Figura 2 apresenta a nuvem de pontos 3D obtida com o equipamento LST, que foi utilizada nesse trabalho. Para cobrir parte da edificação, o LST foi posicionado em quatro pontos de vista diferentes. Para a avaliação da acurácia absoluta dos resultados de localização, todas as coordenadas planimétricas dos diferentes pontos de vista foram levantadas usando uma estação total Leica TC2002, com uma precisão nominal linear de 1 mm e uma precisão nominal angular de 0,5". Na Tabela 1, é apresentado o número total de pontos perfilados em cada nuvem de pontos, a quantidade de planos extraídos automaticamente, bem como o limiar de distância máxima pré-estabelecido para o funcionamento do algoritmo RANSAC. Também vale ressaltar que em todos os experimentos foram admitidos 1000 repetições para o algoritmo RANSAC.

Tabela 1: Conjunto de experimentos elaborados com dados reais. 

Na Figura 2, cada plano segmentado é representado por uma cor diferente. Pode-se perceber, visualmente, que o algoritmo RANSAC considerou alguns objetos como planos, cuja densidade de pontos é baixa, tais como os beirais das edificações, a vegetação, o latão de lixo etc. Para estes planos não foram estabelecidas nenhuma correspondência. Outra observação que pode ser feita, na Figura 2, é que os pontos contidos no solo foram agrupados em vários conjuntos diferentes. Em outras palavras, o plano do solo foi dividido em vários planos. A principal razão deste efeito é a existência de uma pequena inclinação no nível do solo. Entretanto, para este caso, o algoritmo mostrou um desempenho robusto para o estabelecimento de correspondência entre os planos.

Após os planos serem detectados pelo RANSAC e, devidamente, ajustados pelo MMQ totais foi realizado o processo de correspondência entre os planos. A partir das correspondências estabelecidas, os parâmetros de transformação podem ser estimados para cada par de nuvens de pontos 3D. Consequentemente, pode ser determinada a localização do sensor.

Figura 2: Nuvem de Pontos. (A) Nuvem de Pontos I; (B) Nuvem de Pontos II; (C) Nuvem de Pontos III; (D) Nuvem de Pontos IV 

A Figura 3 mostra a trajetória obtida do SLT empregando o método proposto, assim como as medidas levantadas por meio da estação total supracitada.

Figura 3: Posições do sensor LST calculadas através do método proposto e pela estação total. 

Na Figura 3, a linha em azul representa a trajetória do LST e em vermelho a trajetória da estação total. Para realizar uma análise quantitativa do método proposto foi calculado o erro de verificação. Em outras palavras, a avaliação é feita através da média e desvio padrão da distância entre o centroide dos pontos pertencentes a um determinado plano na nuvem de pontos de pesquisa e seu plano correspondente na nuvem de pontos de referência. Os resultados estatísticos obtidos são apresentados na Tabela 2 e 3.

Tabela 2: Resultados dos registrados de pares de nuvens de pontos. 

Tabela 3: Discrepâncias dos resultados. 

Na Tabela 2, os valores δxe δy são as discrepâncias entre as coordenadas planimétricas das posições de referência (determinadas com estação total) e as obtidas com o método proposto. A raiz quadrada do erro médio (REMQ) da distância e dos resíduos dos vetores normais em cada registro indica a precisão do registro dos dados. Como pode ser verificado na Tabela 2, as diferenças foram calculadas na ordem dos decímetros. Já na Tabela 3 pode ser verificado que registro do par de nuvens de pontos II-I produziu o resultado de pior precisão, uma vez que os REMQ dos resíduos dos parâmetros das distâncias dos planos à origem foram de 0,094 m. Considerando-se a precisão nominal linear do sistema LST, empregado neste trabalho, em torno de 4 mm, estes resíduos são visivelmente grandes. A principal razão para este efeito é o número de pares de planos correspondentes usados no processo de registro do par de nuvens de pontos supracitado, sendo de apenas quatro planos. Os resultados obtidos para os demais registros, isto é, pares de nuvens III-II e IV-III apresentaram valores de REMQ dos parâmetros de distâncias dos planos à origem de 0,00361 m e 0,018 m. Neste caso, as precisões são as mesmas, uma vez que foi empregado, para cada registro, um número superior a cinco planos correspondentes.

Também pode ser observado na Tabela 2 que as discrepâncias entre as coordenadas obtidas com a estação total e as determinadas pelo método proposto apresentam médias de -0,02153 m e -0,06266 m, respectivamente. Isto é, as discrepâncias estão na ordem dos centímetros. Como esperado, o registro da nuvem de pontos II também resultou nos maiores valores de discrepâncias. Este fato está claramente apresentado na Figura 3 (canto superior direito). Isto demonstra a necessidade de no mínimo seis planos correspondentes para obter boa qualidade no registro dos pares de nuvem de pontos.

A Figura 4 apresenta uma vista de aérea de parte da edificação após o registro dos pares de nuvem de pontos. Uma análise qualitativa dos resultados obtidos também pode ser efetuada através de inspeção visual. Esta inspeção visual mostra a qualidade da transformação realizada com os parâmetros obtidos, uma vez que pode ser verificada a sobreposição dos pontos contidos em alguns objetos, como por exemplo, as fachadas.

Figura 4: (A) Vista aérea da edificação completa após o registro dos pares de nuvem de pontos; (B) Vista aérea de uma determinada fachada da edificação após o registro dos pares de nuvem de pontos 

5. Conclusões e Recomendações para Trabalhos Futuros

Este trabalho apresentou um método de registro automático de nuvens de pontos derivados do sistema LST. Para investigar a potencialidade do método proposto, experimentos com dados reais foram realizados.

O processo de correspondência plano-a-plano proposto apresenta, como principal característica, o emprego de restrições geométricas. Os resultados obtidos mostraram que as correspondências entre planos, empregando as restrições geométricas, são obtidas com eficiência e robustez, uma vez que os falsos positivos foram devidamente detectados e removidos do processo. Outra vantagem do método proposto em relação aos métodos apresentados no estado da arte é a busca rápida das correspondências devido ao menor número de planos, quando comparados com métodos baseados em abordagens ponto-a-ponto. No entanto, uma restrição do método de correspondência proposto é a exigência de um número suficiente de quatro pares de planos disponíveis em cada nuvem de pontos e ainda que pelo menos três destes planos se interceptem em um único ponto, de modo a haver uma solução única para a estimativa dos parâmetros de translação do sensor.

A contribuição mais importante do método proposto é a obtenção da solução por meio de uma única iteração, reduzindo o custo computacional. Desta forma o algoritmo de correspondência entre planos torna-se mais simples e eficiente. Isto se dá através da utilização dos quatérnios para determinar o parâmetro de rotação do sensor LST. Além disso, nenhuma carga de trabalho adicional é imposta durante o processo de registro dos pares de nuvens de pontos. A principal vantagem do emprego de quatérnios é o seu baixo custo computacional comparado a outras representações como, por exemplo, os ângulos de Euler. Por outro lado os quatérnios unitários representam apenas rotações. Outra característica do método proposto, neste trabalho, é a incapacidade de utilização direta de objetos mais complexos, uma vez que apenas superfícies planas são empregadas no processo de registro de pares de nuvem de pontos 3D.

Como trabalhos futuros, serão conduzidas investigações para incluir simultaneamente feições pontuais, lineares e planas no processo de registro de pares de nuvens de pontos 3D e o desenvolvimento de técnicas de análise de consistência global para evitar o acúmulo de erros ao longo do processo de registro.

AGRADECIMENTOS

Os autores externam seus agradecimentos a Capes e ao CNPq pelas bolsas de doutorado e produtividade em pesquisa (processo no. 302644/2013-0) concedidas para o desenvolvimento deste trabalho

REFERÊNCIAS BIBLIOGRÁFICAS

Anton, Howard; Rorres, Chris. Álgebra linear com aplicações. Tradução de Claus Ivo Doering.8. ed. Porto Alegre: Bookman, 2001. [ Links ]

Arun, K. S., Huang, T. S., and Blostein, S. D., "Least-squares fitting of two 3-D point sets", IEEE Transactions Pattern Analysis Machine Intelligence., vol. 9, no. 5, (1987) 698-700, Accessed April 17, 2014. doi: 10.1109/TPAMI.1987.4767965 [ Links ]

AL-Durgham, M., Habib, A., and Kwak, E., "A framework for the registration and segmentation of heterogeneous LIDAR data", ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, Turkey, (2013). no. 2, 135 - 145. Accessed March 12, 2014. [ Links ]

Bae, Kwang-Ho David Belton and Lichti, Derek D., "A Closed-Form Expression of the Positional Uncertainty for 3D Point Clouds", IEEE Transactions on Pattern Analysis & Machine Intelligence, vol.31, (2008) no. 4, 577-590, Accessed March 11, 2014. doi:10.1109/TPAMI.2008.116. [ Links ]

Besl, Paul. J., and McKay, Neil. D., "A method for registration of 3-D shapes", IEEE Transactions Pattern Analysis Machine Intelligence, vol. 14, (1992) no. 2, 239-256, Accessed March 07, 2014. doi: 10.1109/34.121791 [ Links ]

Brenner, Claus and Dold, Christoph, "Automatic relative orientation of terrestrial laser scans using planar structures and angle constraints", ISPRS Workshop on Laser Scanning 200, Espoo, Finland, (2007) 84-89, Accessed April 22, 2014. [ Links ]

Chen, Yang; Medioni, Gerard (1991). "Object modelling by registration of multiple range images". Image and Vision Computing, Vol. 10, Newton, MA, USA, (1992) no. 3, 145-155. Accessed November 05, 2013. doi:10.1016/0262-8856(92)90066-C. [ Links ]

Dalmolin, Quintino. Ajustamento por Mínimos Quadrados. Curitiba: Imprensa Universitária (UFPR), 2004. [ Links ]

Fischler, Martin A., e Bolles, Robert C., "Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography," Communications of the ACM, Vol. 24, (1981) no. 6, 381-395. Accessed May 09, 2014. doi:10.1145/358669.358692 [ Links ]

Fontanelli, Daniele., Ricciato, Luigi and Soatto, Stefano, "A fast RANSAC-based registration algorithm for accurate localization in unknown environments using LIDAR measurements", IEEE International Conference on Automation Science and Engineering (CASE), (2007), 597-602. Accessed April 19, 2014. doi: 10.1109/COASE.2007.4341827 [ Links ]

Galo, Maurício and Tozzi, Clésio L., "A representação de matrizes de rotação e o uso de quatérnios em Ciências Geodésicas". Série em Ciências Geodésicas. Curitiba: Editora da UFPR, Vol. 1, (2001), 214-231. Accessed November 21, 2013 [ Links ]

Golub, Gene H. and Loan, Charles. F. Van., "An analysis of the total least squares problem", SIAM Journal of Numerical Analysis, vol. 17, (1980) no. 6, 883-893. Accessed April 19, 2014. doi: 10.1137/0717073 [ Links ]

Gomes, J.; e Velho, L., Fundamentos da Computação Gráfica. Série de Computação e Matemática, IMPA, Rio de Janeiro, 2008. [ Links ]

Horn, Berthold. K. P., "Closed-form solution of absolute orientation using unit quaternions". Journal of the Optical Society of America vol. 4, (1987) 629-642. Accessed March 26, 2014. doi: 10.1364/JOSAA.4.000629 [ Links ]

Jaw, J., and Chuang, T., "Feature-based registration of terrestrial LIDAR point clouds". International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences (ISPRS) Beijing, China, (2008) 303-308. Accessed March 11, 2014. [ Links ]

Kim, A., and Golnaraghi, M. F., "A quaternions-based orientation estimation algorithm using a inertial measurement unit". IEE Position Location Navigation Symp. (2004) 268-272. Accessed March 07, 2014. [ Links ]

Khoshelham, Kourosh, "Automated localization of a laser scanner in indoor environments using planar objects". International Conference on Indoor Positioning and Indoor Navigation (IPIN), Zürich, Switzerland, (2010) 1-7. Accessed March 07, 2014. [ Links ]

Khoshelham, Kourosh and Gorte, Ben. G. H., "Registering point clouds of polyhedral buildings to 2D maps", Proceedings of the 3rd ISPRS International Workshop 3D-ARCH 2009: 3D Virtual Reconstruction and Visualization of Complex Architectures. Vol. XXXVIII-5/W1, Trento, Italy, (2009). Accessed March 07, 2014. [ Links ]

Park, Soon-Yong and Subbarao, Murali, "An accurate and fast point-to-plane registration technique", Pattern Recogn. Lett, (2003) 2967-2976, Accessed June 11, 2013. doi: 10.1016/S0167-8655(03)00157-0 [ Links ]

Pathak, Kaustubh; Birk, Andreas; Vaskevicius, Narunas and Poppinga, Jann. Fast registration based on noisy planes with unknown correspondences for 3-D mapping, IEEE Transactions Robotics, vol. 26, (2010), no. 3, 424-441, Accessed March 26, 2014. doi: 10.1109/TRO.2010.2042989 [ Links ]

Rabbani, Tahir; Dijkman, Sander; Heuvel, van den and Vosselman, George, "An Integrated Approach for Modelling and Global Registration of Point Clouds", ISPRS journal of Photogrammetry and Remote Sensing Vol. 61, (2007), 355-370. Accessed March 12, 2014 [ Links ]

Rusu, Radu Bogdan and Cousins, Steve, "3d is here: Point cloud library (PCL)". IEEE International Conference on Robotics and Automation (ICRA). Shanghai (2011) 1-4. Accessed March 12, 2014. doi: 10.1109/ICRA.2011.5980567 [ Links ]

Segal, Aleksandr V., Haehnel, Dirk and Thrun, Sebastian, "Generalized-ICP". Proceedings of Robotics: Science and Systems, Seattle, USA. (2009) 2660-2667. Accessed March 12, 2014. doi: 10.1109/ICRA.2011.5980322 [ Links ]

Shen, Y. Z., Chen, Y. and Zheng, D. H., "A quaternions-based geodetic datum transformation algorithm". Journal of Geodesy, (2006), 233-239. Accessed March 12, 2014. doi: 10.1007/s00190-006-0054-8 [ Links ]

Steinbruch, Alfredo; Winterle, Paulo. Geometria analítica. 2. ed. São Paulo: Pearson Makron Books, 2006. [ Links ]

Taguchi, Y.; Jian, Y.-D.; Ramalingam, S.; e Feng, C., "Point-plane SLAM for hand-held 3D sensors". IEEE International Conference on Robotics and Automation (ICRA), (2013), 5182 - 5189. Accessed March 27, 2014. doi: 10.1109/ICRA.2013.6631318 [ Links ]

Wang, L., Sohn, G., "Automatic co-registration of terrestrial laser scanning data and 2D floor plan". Proceedings of Canadian Geomatics Conference 2010 and ISPRS COM I Symposium, Calgary, Canada, (2010) 14-18. Accessed March 12, 2014. [ Links ]

Walker, Michael W., Shao, Lejun; Volz, Richard A., Estimating 3-D location parameters using dual quaternions. CVGIP: Image Understanding 54 (1991), 358-367. Accessed March 27, 2014. doi:10.1016/1049-9660(91)90036-O [ Links ]

Recebido: Julho de 2014; Aceito: Maio de 2015

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License