Método para determinar se um segmento de linha é uma borda externa de uma triangulação de Delauney?

5

Eu criei uma triangulação de Delauney de um conjunto de pontos. Agora eu quero iterar sobre a triangulação e remover quaisquer segmentos / bordas de linha no exterior para os quais os seguintes itens são verdadeiros:

  1. A borda externa é a hipotenusa do triângulo
  2. Chame a hipotenusa de "base" e determine a "altura" do triângulo. Então, se base / altura > x (onde x é algum valor de critério, digamos 4), exclua essa aresta / triângulo (dependendo da estrutura de dados / definições de objeto).

Essencialmente, estou tentando remover bordas da triangulação para criar uma forma mais característica, excluindo bordas longas, finas e do lado de fora da triangulação. Eu poderia adicionar um cheque mais tarde para o tamanho do triângulo em relação à média do resto também.

Minha pergunta é: existe uma maneira fácil de verificar se uma determinada borda está no exterior da triangulação?

Eu imaginei que poderia incorporar ao próprio processo de triangulação um design que faria com que cada segmento de linha fosse seu próprio objeto, então cada segmento de linha "conhecesse" seus polígonos vizinhos e, se tivesse um slot vizinho vazio, seria externo . No entanto, isso parecia ineficiente, e o código atualmente expõe um monte de objetos triangulares com os índices dos vértices armazenados. Calcular os segmentos de linha não deve ser muito ruim, mas eu só quero fazer o processo de exclusão em bordas externas em oposição a todas elas.

Não tenho certeza se essa pergunta pertence a matemática ou programação, peço desculpas se ela não deveria estar aqui ou se eu forneci pouca ou muita informação.

Estou trabalhando em Python, mas minha preocupação aqui é sobre a abordagem / algoritmo teórico em oposição ao código, então não acho necessário incluir código.

EDIT para anotar: Eu irei querer iterar sobre o exterior várias vezes, no caso em que, após excluir um triângulo externo, a primeira rodada revela um novo exterior que se ajusta aos parâmetros de exclusão. Portanto, eu gostaria de evitar recalcular completamente as bordas externas de cada vez - portanto, a pergunta sobre um rápido e fácil, se existir.

    
por acm_myk 16.01.2015 / 18:20
fonte

3 respostas

5

Uma triangulação do seu conjunto de pontos fornece a cada ponto uma lista de pontos "adjacentes", os vizinhos aos quais um ponto está conectado.

Depois de calcular a triangulação, você pode determinar as bordas externas aplicando um algoritmo semelhante ao clássico algoritmo de embrulho de presente para cascos convexos. Você começa com, por exemplo, o ponto "mais à direita" (que é parte do casco convexo e, portanto, na borda externa). Em seguida, vá de um ponto a outro usando apenas as bordas externas da triangulação. A única diferença para o algoritmo de "casco convexo" é, quando tendo pontos p_ (n-1) e p_ (n), você escolhe p_ (n + 1) como o ponto entre todos os pontos adjacentes minimizando o ângulo entre p_ (n) - > p_ (n-1) e p_ (n) - > p_ (n + 1).

Espero que isso ajude.

    
por 16.01.2015 / 19:14
fonte
3

É uma sorte poder começar com uma lista de objetos triangulares, porque isso facilita:

  • flatMap a lista de triângulos em uma lista de arestas. / li>
  • Suas bordas externas aparecem apenas uma vez na lista.

Você pode ver como isso funciona olhando para uma foto. Cada borda interna é compartilhada por exatamente dois triângulos. Cada borda externa pertence exatamente a um triângulo.

Quando você remove uma aresta, as outras duas arestas desse triângulo sempre se tornam arestas externas, se já não estiverem. Basta adicioná-los ao conjunto de arestas externas antes de remover o triângulo da lista. Isso pode ser acelerado mantendo um mapa das bordas de volta aos triângulos.

    
por 16.01.2015 / 20:24
fonte
-2

Para futuros leitores, aqui está uma função python (usando numpy) que cria uma lista dos índices de pontos no limite.

Estou começando com uma rede geral de pontos conectados (não necessariamente triangulados), então meus dados estão na forma de (xy, NL, KL) com NL sendo a lista de vizinhos e KL sendo a lista de conectividade (veja abaixo) .

def extract_boundary_from_NL(xy,NL,KL):
    '''
    Extract the boundary of a 2D network (xy,NL,KL).

    Parameters
    ----------
    xy : #pts x 2 float array
        point set in 2D
    NL : #pts x max # nearest neighbors int array
        Neighbor list. The ith row contains the indices of xy that are the bonded pts to the ith pt.
        Nonexistent bonds are replaced by zero.
    KL : #pts x max # nearest neighbors int array
        Connectivity list. The jth column of the ith row ==1 if pt i is bonded to pt NL[i,j].
        The jth column of the ith row ==0 if pt i is not bonded to point NL[i,j].

    Returns
    ----------
    boundary : #points on boundary x 1 int array
        indices of points living on boundary of the network

    '''
    # Initialize the list of boundary indices to be larger than necessary
    bb = np.zeros(len(xy), dtype=int)

    # Start with the rightmost point, which is guaranteed to be 
    # at the convex hull and thus also at the outer edge.
    # Then take the first step to be along the minimum angle bond
    rightIND = np.where(xy[:,0]== max(xy[:,0]))[0]
    # If there are more than one rightmost point, choose one
    if rightIND.size >1:
        rightIND = rightIND[0]
    # Grab the true neighbors of this starting point
    neighbors = NL[rightIND,np.argwhere(KL[rightIND]).ravel()]
    # Compute the angles of the neighbor bonds 
    angles = np.mod(np.arctan2(xy[neighbors,1]-xy[rightIND,1],xy[neighbors,0]-xy[rightIND,0]).ravel(), 2*np.pi)
    # Take the second particle to be the one with the lowest bond angle (will be >= pi/2)
    nextIND = neighbors[angles==min(angles)][0]
    bb[0] = rightIND

    dmyi = 1
    # as long as we haven't completed the full outer edge/boundary, add nextIND
    while nextIND != rightIND:
        bb[dmyi] = nextIND
        n_tmp = NL[nextIND,np.argwhere(KL[nextIND]).ravel()]
        # Exclude previous boundary particle from the neighbors array
        # since its angle with itself is zero!
        neighbors = np.delete(n_tmp, np.where(n_tmp == bb[dmyi-1])[0])
        angles = np.mod( np.arctan2(xy[neighbors,1]-xy[nextIND,1],xy[neighbors,0]-xy[nextIND,0]).ravel() \
                - np.arctan2( xy[bb[dmyi-1],1]-xy[nextIND,1], xy[bb[dmyi-1],0]-xy[nextIND,0] ).ravel(), 2*np.pi)
        nextIND = neighbors[angles==min(angles)][0]
        dmyi += 1

    # Truncate the list of boundary indices
    boundary = bb[0:dmyi]

    return boundary

Para colocar uma triangulação no formulário (NL, KL) para usar o acima, você poderia fazer:

BL = TRI2BL(TRI)
NL,KL = BL2NLandKL(BL,NN=10)



def BL2NLandKL(BL, NP='auto', NN=6):
    '''Convert bond list (#bonds x 2) to neighbor list (NL) and connectivity list (KL) for lattice of bonded points.
    Returns KL as ones where there is a bond and zero where there is not.
    (Even if you just want NL from BL, you have to compute KL anyway.)
    Note that this makes no attempt to preserve any previous version of NL, which in the philosophy of these simulations should remain constant during a simulation.
    If NL is known, use BL2KL instead, which creates KL according to the existing NL. 

    Parameters
    ----------
    BL : array of dimension #bonds x 2
        Each row is a bond and contains indices of connected points
    NP : int
        number of points (defines the length of NL and KL)
    NN : int
        maximum number of neighbors (defines the width of NL and KL)

    Returns
    ----------
    NL : array of dimension #pts x max(#neighbors)
        The ith row contains indices for the neighbors for the ith point.
    KL :  array of dimension #pts x (max number of neighbors)
        Spring constant list, where 1 corresponds to a true connection while 0 signifies that there is not a connection.
    '''
    if NP=='auto':
        if BL.size>0:
            NL = np.zeros((max(BL.ravel())+1,NN),dtype=np.intc)
            KL = np.zeros((max(BL.ravel())+1,NN),dtype=np.intc)
        else:
            raise RuntimeError('ERROR: there is no BL to use to define NL and KL, so cannot run BL2NLandKL()')
    else:
        NL = np.zeros((NP,NN),dtype=np.intc)
        KL = np.zeros((NP,NN),dtype=np.intc)

    if BL.size > 0:
        for row in BL:
            col = np.where(KL[row[0],:]==0)[0][0]
            NL[row[0],col] = row[1]
            KL[row[0],col] = 1
            col = np.where(KL[row[1],:]==0)[0][0]
            NL[row[1],col] = row[0]
            KL[row[1],col] = 1

    return NL, KL

def TRI2BL(TRI):
    '''
    Convert triangulation index array (Ntris x 3) to Bond List (Nbonds x 2) array.

    Parameters
    ----------
    TRI : Ntris x 3 int array
        Triangulation of a point set. Each row gives indices of vertices of single triangle.

    Returns
    ----------
    BL : Nbonds x 2 int array 
        Bond list

    '''
    # each edge is shared by 2 triangles unless at the boundary.
    # each row contains 3 edges.
    # An upper bound on the number bonds is 3*len(TRI)
    BL = np.zeros((3*len(TRI),2),dtype = int)

    dmyi = 0
    for row in TRI:
        BL[dmyi] = [row[0], row[1]]
        BL[dmyi+1] = [row[1], row[2]]
        BL[dmyi+2] = [row[0], row[2]]
        dmyi += 3

    # Sort each row to be ascending
    BL_sort = np.sort(BL, axis=1)
    BLtrim = unique_rows(BL_sort)

    return BLtrim
    
por 08.02.2016 / 14:47
fonte