Il Bar dell'Ingegneria

Algoritmi: Punto interno o esterno ad un poligono

« Older   Newer »
 
  Share  
.
  1.     +1   -1
     
    .
    Avatar

    Advanced Member

    Group
    Administrator
    Posts
    8,163
    Reputation
    +294

    Status
    Offline
    inside001


    RIleggendo questo vecchio ed appassionate topic ( https://bar-ingegneria.forumfree.it/?t=8970619 ) risalente a qualche anno fa (giugno del 2006), mi sono ricordato di aver letto in questi giorni una pagina di Paul Bourke dove è riportato l'algoritmo che risolve la questione.

    http://paulbourke.net/geometry/insidepoly/

    Prendendo lo spunto da questo algoritmo ho anche pensato di scrivere un codice minimo in VBA per autocad per stabilire se un punto è interno o esterno ad un poligono sfruttando una particolare funzione di alcuni oggetti grafici di Autocad

    Edited by afazio - 26/3/2013, 14:37
     
    Top
    .
  2.     +1   -1
     
    .
    Avatar

    Advanced Member

    Group
    Administrator
    Posts
    8,163
    Reputation
    +294

    Status
    Offline
    Il metodo a cui mi riferisco è il intersectswith method di un elemento grafico con altro elemento grafico (nel caso specifico una linea con una polilinea).

    La macro che intendo creare, dovrà consentire la selezione di una polilinea chiusa, quindi il clik su un punto qualsiasi e quindi ottenere l'informazione richiesta (se il punto è interno o esterno) attraverso una finestra di avviso.

    Vediamo la particolarità del metodo .intersectswith. Il metodo richiede come parametri i due elementi grafici di cui cercare le intersezioni e retituisce la lista delle coordinate dei punti di intersezione trovati.
     
    Top
    .
  3.     +1   -1
     
    .
    Avatar

    Advanced Member

    Group
    Administrator
    Posts
    8,163
    Reputation
    +294

    Status
    Offline
    Iniziamo.
    Apriamo un nuovo disegno in autocad.
    Disegniamo una qualsiasi polilinea chiusa sul foglio.
    Adesso entriamo in ambiente programmazione.
    Inseriamo un nuovo modulo
    Inseriamo un nuovo form

    Nel form inseriamo:
    - il pulsante di uscita
    - un pulsante di avvio della macro
    - una serie di testi che informano come si deve procedere

    Al pulsante di uscita (a cui io ho dato i nome "Esci") leghiamo il codice:
    CODICE
    Private Sub Esci_Click()
       Me.hide
    End Sub


    Andiamo nello spazio codice del modulo che abbiamo inserito e scriviamo la procedura che ci visualizza il form

    CODICE
    Sub esegui()
       UserForm1.Show
    End Sub


    Adesso non ci resta che dedicarci al codice che ci permette di sapere se un punto clikkato è interno o esterno ad un poligono che abbiamo clikkato

    Io ho organizzato il dialogo come nella immagine che segue:

    insidepoly01
     
    Top
    .
  4.     +1   -1
     
    .
    Avatar

    Advanced Member

    Group
    Administrator
    Posts
    8,163
    Reputation
    +294

    Status
    Offline
    Volendo mettere in codice l'algoritmo riportato, dato il poligono ed il punto di test, prima costruirsi una semiretta arbitraria con origine nel punto di test (possiamo pertanto scegliere un secondo punto a piacere discosto dal punto test di un deltax=1).
    In autocad è possibile disegnare una semiretta infinita (di costruzione) col comando _ray. La stessa semiretta puo essere disegnata da VBA per autocad col metodo:

    Set oggetto_ray= ThisDrawing.ModelSpace.AddRay(punto_inizio, punto_a_piacere)

    Sfruttando questo comando possiamo far disegnare temporaneamente la semiretta

    Adesso ci seve conoscere il numero delle intersezioni tra questa semiretta e la polilinea. Come anticipato, si ricorre al metodo:

    coordinate_punti_intersezione =Oggetto_ray.IntersectWith(Oggetto_Polilinea, acExtendNone)

    Il parametro acExtendNone serve per evitare di trovare intersezioni sugli ideali lati estesi sia della semiretta che del poligono.

    LA variabile coordinate_punti_intersezione dove saranno depositate una dietro l'altra le coordinate x,y,z di tutti i punti di intersezione deve essere di tipo Variant e conterrà un numero di elementi pari a 3*N dove N sono il numero delle intersezioni.

    Abbiamo tutto. Associamo al tato Avvia il seguente codice:

    CODICE
    Private Sub Avvia_Click()

    Dim Semiretta As AcadRay
    Dim Oggetto_generico As AcadObject
    Dim Polilinea As AcadPolyline

    Dim VarPunto As Variant
    Dim Punto_click As Variant

    Dim coOrdinate As Variant
    Dim numero_intersezioni As Long

    Dim punto1(0 To 2) As Double
    Dim punto2(0 To 2) As Double

    Dim errore As Boolean

    UserForm1.hide
    errore = False
    While Not (errore)
       On Error Resume Next
       ThisDrawing.Utility.GetEntity Oggetto_generico, VarPunto, "Seleziona una polilinea Chiusa"
       If Err <> 0 Then
           Err.Clear
           errore = True
       Else
           If Oggetto_generico.EntityName = "AcDbPolyline" Then
               MsgBox ("l'oggetto clickato è una polilinea")
               ' devo vedere se la polilinea è chiusa
               Set Polilinea = Oggetto_generico
               If Polilinea.Closed Then
                   MsgBox ("la polilinea è chiusa")
                   Punto_click = ThisDrawing.Utility.GetPoint(, "clik sul punto")
                   ' stabilisco i due punti per disegnare la semiretta
                   ' il primo punto è il punto di test e coincide col punto clik
                   punto1(0) = Punto_click(0)
                   punto1(1) = Punto_click(1)
                   punto1(2) = Punto_click(2)
                   ' il secondo punto è arbitrario, io ho scelto un punto sulla stessa
                   ' ordinata dal primo e distante di una unità
                   punto2(0) = punto1(0) + 1
                   punto2(1) = punto1(1)
                   punto2(2) = punto1(2)
                   ' disegno la semiretta che parte dal primo punto, passa dal secondo e
                   ' si estende all'infinito
                   Set Semiretta = ThisDrawing.ModelSpace.AddRay(punto1, punto2)
                   Update
                   
                   'mi determino le intersezioni tra la semiretta e la polilinea
                   coOrdinate = Semiretta.IntersectWith(Oggetto_generico, acExtendNone)
                   ' se vi sono intersezioni allora il vettore delle coordinate non è vuoto
                   
                   If VarType(coOrdinate) <> vbEmpty Then
                       numero_intersezioni = (UBound(coOrdinate) + 1) / 3
                       MsgBox ("trovate " & numero_intersezioni & " intersezioni ")
                   Else
                       numero_intersezioni = 0
                       MsgBox ("nessuna intersezione trovata")
                   End If
                   ' determino se il numero è pari o dispari attraverso
                   ' il resto della divisione intera
                   If numero_intersezioni Mod 2 = 0 Then
                       MsgBox ("il punto clickato è esterno alla polilinea")
                   Else
                       MsgBox ("il punto clickato è interno alla polilinea")
                   End If
                   'cancello la semiretta
                    Semiretta.Delete
               End If
          End If
       End If
       Update
    Wend
    UserForm1.Show
    End Sub


    Salviamo ed eseguiamo la macro vedendone il funzionamento
    Tutto funziona perfettamente quando i punti test sono nettamente fuori o nettamente dentro il poligono, mentre comincia a dare risultati strani quando il punto test è su un vertice.

    Comunque questo codice non differenzia i casi possibili come ero riuscito a fare col codice scritto sei anni fa in VBA per excel:
    - punto interno
    - punto esterno
    - punto su un vertice
    - punto su un lato.

    Se avessi necessità di una funzione che mi distingua lo stato del punto rispetto al poligono, dovrei pensare ad aggiustamenti

    Edited by afazio - 21/8/2012, 20:57
     
    Top
    .
  5.     +1   -1
     
    .
    Avatar

    Advanced Member

    Group
    Administrator
    Posts
    8,163
    Reputation
    +294

    Status
    Offline
    Dopo essermici sbattuto parecchio nel tentativo di districare tutti i vari casi adottando l'algoritmo del Ray Scan, scivolando sempre proprio quando pensavo di averlo risolto, su quel caso imprevisto, ho deciso di abbandonarlo e cercare la soluzione adottando l'algoritmo della somma degli angoli sottesi dai lati rispetto al punto di cui vogliamo sapere se è interno o esterno.

    L'algoritmo si basa sulla constatazione che sommando gli angoli sottesi dai lati rispetto al punto, otteniamo un angolo giro (positivo o negativo a seconda del clock del poligono) se il punto è interno, mentre otteniamo zero se il punto è esterno.

    basta dare un'occhiata alle due figure che seguono per rendersene immediatamente conto

    Punto interno al poligono
    inside001

    Punto esterno al poligono
    inside002



    Quindi la chiave dell'algoritmo consiste nel riuscire a determinarci l'angolo formato tra i due vettori P-Pi e P-Pk, essendo Pi e Pk gli estremi del genrico lato del poligono.

    Propongo il codice ed a seguire lo commento.

    CODICE
    Function IsInsidePoly(Polygon() As TipoPunto, Point As TipoPunto) As Boolean

    Dim i As Integer
    DIm k as Integer
    Dim AngoloTotale As Double
    Const Tolleranza As Double = 0.0000001
       
    AngoloTotale = 0
    For i = LBound(Polygon) To UBound(Polygon)
       k = k + 1: If i = UBound(Polygon) Then k = LBound(Polygon)
       AngoloTotale = AngoloTotale + AngoloTra2Vettori(Point, Polygon(i), Polygon(k))
    Next
    ' l'angolo totale sarà 2*pi o -2*pi se il punto è interno,  oppure zero se il punto è esterno
    IsInsidePoly = (Abs(AngoloTotale) > Tolleranza)

    End Function


    come vedete la cosa è abbastanza semplice e lascia tutto il compito alla funzione:

    AngoloTra2Vettori(Point, Polygon(i), Polygon(k))

    che deve calcolare l'angolo (con segno) formato tra i due vettori Point-Polygon(i) e Point-Polygon(k)




    L'angolo sotteso da un lato rispetto al punto d'esame non puo' essere mai ne maggiore di un angolo piatto ne minore (nel caso di clock al contrario) dell'angolo piatto cambiato di segno

    cioè vale la relazione:

    -pi <= alfa <= pi
    quindi la funzione AngoloTradueVettori deve essere in grado di restituirmi un angolo con segno compreso tra i due valori specificati.

    La funzione è la seguente e simula la funzione di foglio Atan2 (che però è assente in VBA)
    CODICE
    '------------------------------------------------------------------------------------------
    ' Calcola l'angolo formato tra  P1-P-P2. (calcola l'angolo in P)
    ' La funzione ritorna un valore dell'angolo ridotto compreso tra -pigreco e +pigreco.
    ' Input: i tre punti
    ' output: L'angolo con segno formato tra i due vettori
    ' Autore: Afazio - Ottobre 2012
    '------------------------------------------------------------------------------------------
    Function AngoloTra2Vettori(P1 As TipoPunto, P2 As TipoPunto, P As TipoPunto) As Double

    Dim ProdottoScalare As Double
    Dim ProdottoVettore As Double
       ' calcolo il prodotto scalare e vettoriale
       ProdottoScalare = ProdottoScalare3P(P1, P2, P)
       ProdottoVettore = ProdottoVettore3P(P1, P2, P)

       ' Calcolo l'angolo tra i due vettori
       AngoloTra2Vettori = ATan2(ProdottoVettore, ProdottoScalare)
    End Function


    Qui è da notare il ricorso ai prodotti tra due vettori in particolare al prodotto vettoriale ed a quello scalare.

    inside003

    Ricordando che si ha:

    modulo del prodotto vettoriale: ProdottoVettoriale = |v1|*|v2|*sen(alfa
    Modulo del prodotto scalare : ProdottoScalare= |v1|*|v2|*cos(alfa)

    Dividendo membro a membro la prima equazione per la seconda

    abbiamo:

    tan(alfa) = ProdottoVettoriale/ProdottoScalare

    Ma i moduli dei due prodotti li possiamo calcolare come somme o differenze di prodotti tra le componenti dei due vettori:

    CODICE
    '---------------------------------------------------------------------
    ' la funzione che segue restituisce il modulo del prodotto scalare
    ' tra due vettori definiti attraverso tre punti P1P2 e P1P
    ' input: tre punti secondo struttura TipoPunto
    '---------------------------------------------------------------------
    Public Function ProdottoScalare3P(P1 As TipoPunto, P2 As TipoPunto, P As TipoPunto) As Double
       ProdottoScalare3P = (P2.X - P1.X) * (P.X - P1.X) + (P2.Y - P1.Y) * (P.Y - P1.Y)
    End Function
    '---------------------------------------------------------------------


    CODICE
    '---------------------------------------------------------------------
    ' la funzione che segue restituisce il modulo del prodotto vettoriale
    ' tra due vettori definiti attraverso tre punti P1P2 e P1P
    ' input: tre punti secondo struttura TipoPunto
    '---------------------------------------------------------------------
    Public Function ProdottoVettore3P(P1 As TipoPunto, P2 As TipoPunto, P As TipoPunto) As Variant
       ProdottoVettore3P = (P2.X - P1.X) * (P.Y - P1.Y) - (P2.Y - P1.Y) * (P.X - P1.X)
    End Function
    '---------------------------------------------------------------------


    Basta mettere insieme il tutto ed il codice è fatto.

    Enjoy

    Edited by afazio - 11/10/2012, 15:52
     
    Top
    .
  6.     +1   -1
     
    .
    Avatar

    Advanced Member

    Group
    Administrator
    Posts
    8,163
    Reputation
    +294

    Status
    Offline
    Infine ecco la Simulazione della funzione di foglio Atan2

    CODICE
    '------------------------------------------------------------------------------------------------
    ' la funzione ATan2() restituisce l'angolo avente tangente pari a Numeratore/Denominatore
    ' Il valore ottenuto è compreso tra  pigreco e -pigreco.
    ' Praticamente simula la funzione di foglio Atan2()
    '---------------------------------------------------------------------------------------------------
    Function ATan2(Numeratore As Double, Denominatore As Double) As Double
    Dim Angolo As Double
    Dim Pi As Double
    Pi = 4 * Atn(1)
    Const Tolleranza As Double = 0.0000001
       
       '  controllo che il denominatore sia diverso da zero
       If Abs(Denominatore) < Tolleranza Then
           Angolo = Pi / 2
       Else
           Angolo = Abs(Atn(Numeratore / Denominatore))
       End If
       ' controllo se sono nel  quadrante 2 o 3.
       If Denominatore < 0 Then
           ' l'angolo è > pi/2 oppure < -pi/2.
           Angolo = Pi - Angolo
       End If
       ' controllo se sono nel  quadrante 3 o 4.
       If Numeratore < 0 Then
           Angolo = -Angolo
       End If
       ATan2 = Angolo
       
    End Function
     
    Top
    .
  7.     +1   -1
     
    .
    Avatar

    Advanced Member

    Group
    Administrator
    Posts
    8,163
    Reputation
    +294

    Status
    Offline
    Qui trovate il solito file di test con aggiunte le ultime funzioni.

    https://www.box.com/s/p252utwnkykoyykje52j
     
    Top
    .
  8.     +1   -1
     
    .
    Avatar

    Advanced Member

    Group
    Administrator
    Posts
    8,163
    Reputation
    +294

    Status
    Offline
    Il caso di un poligono definito in senso antiorario col punto interno
    angolie

    Quando il punto "torna indietro", l'angolo calcolato è negativo (angolo in rosso) e si sottrae alla somma. Superata l'ansa. continuando a seguire l'ordine dei vertici, l'angolo torna a sommarsi
     
    Top
    .
  9.     +1   -1
     
    .
    Avatar

    Advanced Member

    Group
    Member
    Posts
    2,939
    Reputation
    +187

    Status
    Offline
    Grazie dello schema afazio.

    Ho anche provato, con la stessa figura a posizionare il punto all'esterno del poligono, ma nel 'becco'. In quel caso con gli opportuni segni la sommatoria degli angoli è nulla c.v.d.
     
    Top
    .
  10.     +1   -1
     
    .
    Avatar

    Advanced Member

    Group
    Administrator
    Posts
    8,163
    Reputation
    +294

    Status
    Offline
    Lo stesso poligono col punto esterno

    angoli01

    La somma algebrica di tutti gli angoli negativi (quelli in bianco) e di tutti gli angoli positivi (quelli in rosso) da come risultato zero. Provare a percorrere con la mente il poligono sommando l'angolo in bianco e sottraendo quello in rosso.

    CITAZIONE (zax2010 @ 12/10/2012, 10:59) 
    Grazie dello schema afazio.

    Ho anche provato, con la stessa figura a posizionare il punto all'esterno del poligono, ma nel 'becco'. In quel caso con gli opportuni segni la sommatoria degli angoli è nulla c.v.d.

    Se poi vai a guardare il codice necessario per l'implemetazione di questo algoritmo noti che si riduce davvero a poche righe. Con l'algoritmo del Ray scan il processo è piu lungo poiche oltre a dover determinare le intsersezioni con tutti i lati dei poligoni hai i problemi nei casi in cui: il ray incoccia un vertice, l'intersezione sara contata una volta appartenente ad un lato ed una volta appartenente al successivo e devi pertanto inserire qualche controllo che ne escluda uno, quando il ray è coincidente con un lato hai ancora due intersezioni, una nel vertice del lato precedente ed una nel vertice del lato successivo.
    Ovviamente tutto si può implementare anche questo algoritmo ma a costo di qualche controllo in piu'

    Edited by afazio - 12/10/2012, 11:32
     
    Top
    .
  11.     +1   -1
     
    .
    Avatar

    Advanced Member

    Group
    Member
    Posts
    2,939
    Reputation
    +187

    Status
    Offline
    Ho fatto di più afazio.
    (essendo come San Tommaso)

    Ho fatto il disegnino in autocad, mi sono misurato tutti gli angoli e li ho sommati algebricamente. 360° per i punti interni, 0° per quelli esterni.
     
    Top
    .
  12.     +1   -1
     
    .
    Avatar

    Advanced Member

    Group
    Administrator
    Posts
    8,163
    Reputation
    +294

    Status
    Offline
    CITAZIONE (zax2010 @ 12/10/2012, 11:11) 
    Ho fatto di più afazio.
    (essendo come San Tommaso)

    Ho fatto il disegnino in autocad, mi sono misurato tutti gli angoli e li ho sommati algebricamente. 360° per i punti interni, 0° per quelli esterni.

    Esatto, ma quel 360 diventa -360 quando cambi il clok del poligono.

    Basta quindi una buona funzione che ti calcola l'angolo con segno sapendo che comunque l'angolo non puo' essere magiore in valore assoluto a 180 ed hai risolto il tutto.
     
    Top
    .
  13.     +1   -1
     
    .
    Avatar

    Advanced Member

    Group
    Administrator
    Posts
    8,163
    Reputation
    +294

    Status
    Offline
    Una funzione in pascal che determina se un punto è interno o esterno ad un poligono:

    http://wiki.lazarus.freepascal.org/Geometry_in_Pascal

    CODICE
    //  The function will return True if the point x,y is inside the polygon, or
    //  False if it is not.
    //
    //  Original C code: http://www.visibone.com/inpoly/inpoly.c.txt
    //
    //  Translation from C by Felipe Monteiro de Carvalho
    //
    //  License: Public Domain


    function IsPointInPolygon(AX, AY: Integer; APolygon: array of TPoint): Boolean;
    var
     xnew, ynew: Cardinal;
     xold,yold: Cardinal;
     x1,y1: Cardinal;
     x2,y2: Cardinal;
     i, npoints: Integer;
     inside: Integer = 0;
    begin
     Result := False;
     npoints := Length(APolygon);
     if (npoints < 3) then Exit;
     xold := APolygon[npoints-1].X;
     yold := APolygon[npoints-1].Y;
     for i := 0 to npoints - 1 do
     begin
       xnew := APolygon[i].X;
       ynew := APolygon[i].Y;
       if (xnew > xold) then
       begin
         x1:=xold;
         x2:=xnew;
         y1:=yold;
         y2:=ynew;
       end
       else
       begin
         x1:=xnew;
         x2:=xold;
         y1:=ynew;
         y2:=yold;
       end;
       if (((xnew < AX) = (AX <= xold))         // edge "open" at left end
         and ((AY-y1)*(x2-x1) < (y2-y1)*(AX-x1))) then
       begin
         inside := not inside;
       end;
       xold:=xnew;
       yold:=ynew;
     end;
     Result := inside <> 0;
    end;
     
    Top
    .
  14.     +1   -1
     
    .
    Avatar

    Member

    Group
    Member
    Posts
    766
    Reputation
    +23

    Status
    Offline
    Io ho recuperato questa funzione in vba che determina se un punto è interno od esterno ad un poligono oltre a segnalare se appartiene ad un lato oppure coincide con un vertice:

    CODICE
    Private Sub INSIDE(XAPEX() As Double, YAPEX() As Double, napex As Integer, xPT As Double, yPT As Double, loc As Integer)
    '
    '      inside          definisce la posizione di un punto rispetto
    '                      a un poligono (interno o esterno)
    '
    '      input           xapex = ascisse vertici del poligono
    '                      yapex = ordinate vertici poligono
    '                      napex = nr. vertici poligono
    '                      xpt & ypt = coord. punto
    '
    '      output          loc
    '
    '      loc = 0         n. vertici < 3
    '      loc = 1         il punto e' esterno al poligono
    '      loc = 2         il punto coincide con un vertice
    '      loc = 3         il punto sta sulla frontiera del poligono
    '      loc = 4         il punto e' interno al poligono
    '
    Dim i As Integer, ANGLE As Double
    Dim theta As Double
    Dim Alfa As Double
    Dim xnew As Double
    Dim ynew As Double
    Dim pi As Double
    pi = 4# * Atn(1#) 'atan(1#)

    XAPEX(napex) = XAPEX(0)   '  forza la chisura del poligono
    YAPEX(napex) = YAPEX(0)   '

    If (napex < 3) Then
       loc = 0          ' poligono inesistente
       Exit Sub
    End If

    For i = 0 To napex
      If (xPT = XAPEX(i) And yPT = YAPEX(i)) Then
           loc = 2
           Exit Sub
      End If
    Next i


    ANGLE = 0#
    For i = 0 To napex - 1
           theta = atan2(YAPEX(i) - yPT, XAPEX(i) - xPT)
           xnew = (XAPEX(i + 1) - xPT) * Cos(theta) + (YAPEX(i + 1) - yPT) * Sin(theta)
           ynew = (YAPEX(i + 1) - yPT) * Cos(theta) - (XAPEX(i + 1) - xPT) * Sin(theta)
           Alfa = atan2(ynew, xnew)
           If Alfa = pi Then
                   loc = 3       ' punto sul lato del poligono
                   Exit Sub
           End If
           ANGLE = ANGLE + Alfa
    Next i
    If (Abs(ANGLE) < 0.0001) Then
           loc = 1              ' Punto esterno
           Exit Sub
    ElseIf (Abs(ANGLE) > pi) Then
           loc = 4              ' punto interno
           Exit Sub
    End If
    End Sub



    Private Function atan2(Y As Double, X As Double) As Double
    '
    '       in Visual basic  la atan2 NON esiste !
    '       uguale alla funzione fortran atan2
    '
    Dim pi
    Dim hpi
    pi = 4 * Atn(1#)    ' pigreco        3.14.....    (180)

    hpi = pi / 2#       ' pigrecomezzi   1.57......   ( 90)
    If X = 0# Then
       If Y > 0# Then
           atan2 = hpi
          ElseIf Y < 0# Then
           atan2 = -hpi
          Else
           atan2 = 0#
       End If
    ElseIf Y = 0# Then
       If X < 0# Then
         atan2 = pi
        Else
         atan2 = 0#
       End If
    Else
       If X < 0# Then
          If Y > 0# Then
            atan2 = Atn(Y / X) + pi
           Else
            atan2 = Atn(Y / X) - pi
          End If
        Else
         atan2 = Atn(Y / X)
       End If
    End If
    End Function


    Altro codice equivalente sviluppato interamente da me su un concetto base di computer graphics (una semiretta che parte dal punto dato se interseca il poligono per un numero pari di volte vuol dire che il punto dato è esterno altrimenti è interno al poligono stesso):
    CODICE
    '///////////////////////////////////////////////////////////////////////////
    'Determina se un punto di coordinate P(Xp,Yp)
    'è interno al poligono indicato tenendo conto anche delle cavità
    '///////////////////////////////////////////////////////////////////////////
    Public Function Punto_Interno_poligono(xP As Double, yP As Double, X() As Double, Y() As Double, nv As Integer) As Boolean

    'xP,yP --> coordinate punto P(Xp,Yp)
    'X()   --> vettore delle ascisse dei vertici del poligono
    'Y()   --> vettore delle ordinate dei vertici del poligono
    'nv    --> numero dei vertici del poligono

       Dim j As Integer
       Dim k As Integer
       Dim xt() As Double
       Dim ntlx As Integer
       Dim yt() As Double
       Dim ntly As Integer
       Dim Xc As Double
       Dim Yc As Double
       Dim denom As Double
       Dim kdx As Integer
       Dim ksx As Integer
       Dim kup As Integer
       Dim kdn As Integer
       

           ntlx = 0
           'calcolo punti di intersezione con il poligono dato
           For j = 1 To nv
           
               If j = nv Then k = 1 Else k = j + 1
           
               'considero solo i lati intersecanti la retta orizzontale passante per yp
               If (yP < Y(j) And yP >= Y(k)) Or (yP >= Y(j) And yP < Y(k)) Then
               
                   'trovo le coordinate del punto di intersezione
                   denom = (Y(k) - Y(j))
                   'le linee parallele non si intersecano mai pertanto non la calcolo
                   'oltretutto genera una divisione per zero non possibile
                   If denom <> 0# Then

                       Xc = ((X(k) - X(j)) * yP + X(j) * Y(k) - X(k) * Y(j)) / denom
                   
                       'aggiungo spazio per un punto
                       ntlx = ntlx + 1
                       ReDim Preserve xt(ntlx)
                       
                       'inserisco il punto
                       xt(ntlx) = Xc
               
                   End If
               
               End If
           
           Next j
                                                                           
           ntly = 0
           'calcolo punti di intersezione con il poligono dato
           For j = 1 To nv
           
               If j = nv Then k = 1 Else k = j + 1
                   
               'considero solo i lati intersecanti la retta verticale passante per xp
               If (xP < X(j) And xP >= X(k)) Or (xP >= X(j) And xP < X(k)) Then
               
                   'trovo le coordinate del punto di intersezione
                   denom = (X(k) - X(j))
                   'le linee parallele non si intersecano mai pertanto non la calcolo
                   'oltretutto genera una divisione per zero non possibile
                   If denom <> 0# Then

                       Yc = ((Y(k) - Y(j)) * xP + Y(j) * X(k) - Y(k) * X(j)) / denom
                   
                       'aggiungo spazio per un punto
                       ntly = ntly + 1
                       ReDim Preserve yt(ntly)
                       
                       'inserisco il punto
                       yt(ntly) = Yc
                   
                   End If
               
               End If
           
           Next j


           'conto i punti di intersezione con il perimetro a destra o a sinistra del punto P
           'se il numero di intersezioni è pari oppure zero significa che il punto P
           'è esterno al contorno.
           kdx = 0
           ksx = 0
           For j = 1 To ntlx
               If xt(j) > xP Then kdx = kdx + 1
               If xt(j) < xP Then ksx = ksx + 1
           Next j

           'conto i punti di intersezione con il perimetro sopra e sotto del punto P
           'se il numero di intersezioni è pari oppure zero significa che il punto P
           'è esterno al contorno.
           kup = 0
           kdn = 0
           For j = 1 To ntly
               If yt(j) > yP Then kup = kup + 1
               If yt(j) < yP Then kdn = kdn + 1
           Next j
                   
           If kdx Mod 2 = 0 Or ksx Mod 2 = 0 Or kup Mod 2 = 0 Or kdn Mod 2 = 0 Then
               Punto_Interno_poligono = False
           Else
               'gestione di caso particolare molto raro
               If punto_su_vertice(xP, yP, X(), Y(), nv) = True Then
                   Punto_Interno_poligono = False
               Else
                   Punto_Interno_poligono = True
               End If
           End If
                   
    End Function

    'Verifichiamo se un punto è anche vertice del poligono
    Private Function punto_su_vertice(xP As Double, yP As Double, Xv() As Double, Yv() As Double, nv As Integer) As Boolean
       
       Dim i As Integer
       
       For i = 1 To nv
           If xP = Xv(i) And yP = Yv(i) Then
               punto_su_vertice = True
               Exit Function
           End If
       Next i
       
       'default
       punto_su_vertice = False

    End Function


    Edited by texitaliano64 - 25/4/2013, 22:17
     
    Top
    .
  15.     +1   -1
     
    .
    Avatar

    Advanced Member

    Group
    Administrator
    Posts
    8,163
    Reputation
    +294

    Status
    Offline
    credo che l'ultimo codice entri in errore in questo caso:

    erroreun

    Infatti, considerando la retta orizzontale passante per il punto abbiamo due intersezioni a destra ed uno a sinistra.
    Per le intersezioni a destra dovremmo avere come risultato "punto esterno" mentre per quelle a sinistra dovremmo avere come risultato "punto interno". Ma è e facile disegnare un poligono che abbia due intersezioni a destra e due a sinistra due verso l'alto e due verso il basso ed avere come risultato "punto esterno" mentre il punto è interno.

    Che succede per esempio nel caso che illustro a seguire?

    errore1p
     
    Top
    .
21 replies since 21/8/2012, 11:07   4906 views
  Share  
.