Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Locate which triangle a point is in #1

Open
HugoGranstrom opened this issue Sep 13, 2021 · 24 comments
Open

Locate which triangle a point is in #1

HugoGranstrom opened this issue Sep 13, 2021 · 24 comments

Comments

@HugoGranstrom
Copy link

I have interests in using Delaunay Triangulation for interpolation, but I can't seem to find any functions suitable for that. locatePoint was the one I found, but it just returns the nearest edge (2 out of 3 points). Am I just blind or is this still not implemented? And in such case, how hard would it be to add? 😄
@bassmit
Copy link

bassmit commented Sep 13, 2021 via email

@bassmit
Copy link

bassmit commented Sep 13, 2021 via email

@HugoGranstrom
Copy link
Author

HugoGranstrom commented Sep 13, 2021 •
edited
Loading

Hi bas (bas ≠ Stefan?), I basically want to get the 3 points making up the triangle, and then use them to calculate an interpolation for a point in the triangle. But it's the triangulation I'm interested in here. It feels like this should be a utility function many users of cdt would want. It feels a bit unnecessary if a user who only want that functionality would have to understand quadedges. It's the exact meanings of the conventions and exactly how they take form in the direction that confuses me. Is the queried point p always on the left side of the directed edge?

So psuedocode would be something like this?

(edge, loc) = dt.locatePoint(p)
next_edge = edge.dNext
points = (edge.org, edge.dest, next_edge.dest)

Edit: reading the doc-comment above dNext made me confused. The dest is the same(?), so I should do next_edge.org to get the third point?

@bassmit
Copy link

bassmit commented Sep 13, 2021 via email

@bassmit
Copy link

bassmit commented Sep 13, 2021 via email

@HugoGranstrom
Copy link
Author

Sorry about the mix-up, it's here https://github.com/StefanSalewski/cdt/blob/master/src/cdt/edges.nim

Yeah I found those, but there are so many variants and it isn't really obvious which combinations of them that I'm looking for in 😅

@bassmit
Copy link

bassmit commented Sep 13, 2021 via email

@HugoGranstrom
Copy link
Author

You sure? I'm so confused by "with the same dest" from below. That suggest dest is one of the verticies of my original edge. I haven't been able to test it on my computer yet (so sorry if I'm just stupid asking the same question over and over 🙈)

# Returns the next Edge about the destination with the same desination.
func dNext*(e: Edge): Edge =

@bassmit
Copy link

bassmit commented Sep 13, 2021 via email

@HugoGranstrom
Copy link
Author

Ok, thanks for the answer. 😄 So the documentation comment is misleading? (It's the comment I'm so confused by now) It says we will get the same dest while you say we get a different dest.

@bassmit
Copy link

bassmit commented Sep 13, 2021 via email

@HugoGranstrom
Copy link
Author

Ahh, that one makes sense! Thanks a lot for taking your time digesting this with me ⭐ I could test this and throw together a PR for the point →3 vertices or 2 edges function if you or Stefan see it as fitting in cdt

@bassmit
Copy link

bassmit commented Sep 13, 2021 via email

@HugoGranstrom
Copy link
Author

🙈 then I'm even more grateful for your help 😄

@StefanSalewski
Copy link
Owner

Thanks for reporting, I will try to investigate soon. Have not touched the code in the last two years, so it will took me some time to become familiar again.

@HugoGranstrom
Copy link
Author

No hurry, I think I've got something to start with now at least 😄

@StefanSalewski
Copy link
Owner

Fine that Mr. bassmit could have helped you already. I will add a function to get the point you desired, I may have to consult the papers on which the module is based for that. Indeed all the code and most comments are based on the various papers. Interpolation was not really a goal of the implementation and was not mentioned in the papers. Main goal was to get a topological structure for PCB trace routing. The funny fact is that the most difficult part for me was to start a new triangulation with only one edge, as that is no real quad edge structure yet. I have seen some different ways to start with one edge, I think the one I finally selected works at least. Personally for me the biggest restriction is that the triangulation can not grow unlimited, but that we always have to start with four points building a bounding box for all other points. CGAL has not this restriction. And as you may have seen, I created about one year ago a cdt2 module with slightly modified API. The initial implementation and its API followed strictly the papers, but I got the feeling when using it for the first time for the topological rubber band router that the API is not well suited for OOP style of coding.

@HugoGranstrom
Copy link
Author

Many thanks :) Delauney Triangulation has many uses, and with this now I should at least be able to do linear interpolation on non-grided data. The plan then is to use even more than the 3 points in the triangle, but I suspect that's something that's not very generally applicable. So I'll probably just have to learn how to jump between edges then. I looked into rolling my own CDT-impl but was scared by it so kudos to you for making it work! For my usecase, we will likely have all the data from the beginning so the bounding box could be calculated from it. But yeah, a bit cumbersome if you have no clue on how big it will be.

Hadn't seen cdt2 actually. Do you recommend using it instead of this repo? Or perhaps it doesn't matter for my case? I will only create the DT, query it, and perhaps jump some edges.

@StefanSalewski
Copy link
Owner

rolling my own CDT-impl

One without constrained edges and without support for dynamically inserting and removing points and edges is very simple, see

https://github.com/Nycto/DelaunayNim

But for a constrained fully dynamically one it is not that easy. I tried to look at some other implementations on github in Java or C++, but the amount of code is generally very large and understanding the implementations is hard. For PCB routing one obviously needed constrained edges, and dynamically inserting and removing points and edges is really a benefit. The Kallmann paper describes all that well, and with a few other papers I was finally able to create the implementation from scratch, but it took me about 3 months full time. And it is not much tested yet, I did a fast test with https://github.com/StefanSalewski/RBR.

Next step would be to work on a plain schematics and PCB editor, see https://github.com/StefanSalewski/SDT before it would make sense to continue with the rubberband router.

The CDT2 module is fully identical in the code, only the API is changed a bit. Finally I may join both APIs, but that will happen not in the next 3 years. I think when you do not OOP style programming with subclassing, then the CDT module is fine for you. As we know Araq does not like OOP, and people like treeform likes to produce modules with do not support subclassing too, Araq thinks that is fine. We have a forum thread about that. But OOP is not really that bad, in some situations avoiding OOP can make live really hard.

@HugoGranstrom
Copy link
Author

That impl wasn't too bad actually. Still, a fair bit of theory to crunch through to get to that point though.

That's neat! Hadn't even heard of PCB routing before so TIL.

will happen not in the next 3 years

3 years of stable API in other words 😎
I'll go with ordinary CDT then. Not getting into the debate on this, I just use whichever tool is best suited for each individual problem.

@StefanSalewski
Copy link
Owner

StefanSalewski commented Sep 13, 2021 •
edited
Loading

I just had a short look at the CDT module. Unfortunately documentation is really missing -- I guess I was too tired after creating that module. The RTree module has at least some more documentation. Explaining CDT is not really easy, and chances that someone really wants to use that modules in very small, at least in the last 2 years no one asked for it.

I think someone who wants to use the CDT module may have to look at the code actually, the code is not that large and there are some comments. Now looking at it I really wonder why it took me so long to write it and to got it working :-)

Unfortunately the locatePoint() proc is not even exported, se we have to add an export marker first, for the LocPos enum too. Point locate is well defined, we get the nearest edge for which the triangular area e.org, e.dest, p is positive. org and dest is origin and destination of the edge. In edges module there are a lot tiny procs to traverse the triangulation.

# Edge Traversal Operators -- counter clockwise (CCW) order

As p has a positive orientation to the returned edge, it seems to make sense to use oNext() to get the next edge with same origin. dest() should then be the opposite point you are locking for. I did a short test without graphical output to proof that, and it seems to make some sense. To do a better verify one can sketch some points on a piece of paper, enter that points then into the code and then do the point locate and verify it on the paper. Or one can graphical output on the computer, I have provided a test program with cairo graphics output in the cdt package, but have to remember how it works.

Here is my short test code:

import cdt/[dt, vectors, edges, types]

const
  minX = 0
  minY = 0
  maxX = 100
  maxY = 100

proc main =
  var dt: DelaunayTriangulation = initDelaunayTriangulation(initVector3(minX, minY), initVector3(maxX, maxY))
  discard dt.insert(Vector(x: 5, y: 10)) # insert a single point
  let cid = dt.insert(Vector(x: 20, y: 30), Vector(x: 30, y: 90)) # insert a constrained line segment
  #discard dt.removeConstraint(cid) # remove an segment

  for e in dt.edges:
    echo e.org, " ", e.dest

  var (e, lp) = dt.locatePoint(Vector(x: 5, y: 10)) # we have an edge (5.0, 10.0) (20.0, 30.0)
  echo e.org, " ", e.dest, " ", lp # 5.0, 10.0) (20.0, 30.0) lpOrg

  (e, lp) = dt.locatePoint(Vector(x: 6, y: 10)) # point below above edge
  echo e.org, " ", e.dest, " ", lp # (5.0, 10.0) (100.0, 0.0) lpFace

  # edge e is well defined, see comment for xlocatePoint(): triArea(e.org.point, e.dest.point, p) >= 0
  # so I would strongly assume that e.onext.dest is the oposite point

  var x = e.onext.dest
  echo x # (20.0, 30.0)

  (e, lp) = dt.locatePoint(Vector(x: 17, y: 30)) # point above  edge
  echo e.org, " ", e.dest, " ", lp # (5.0, 10.0) (20.0, 30.0) lpFace

  x = e.onext.dest
  echo x # (0.0, 100.0)

main()

@StefanSalewski
Copy link
Owner

So the documentation comment is misleading?

I think and hope not really misleading, but I easily get confused.

I think that these are the same points:

e.onext.dest
e.dprev.org

May do a more serious test tomorrow.

@StefanSalewski
Copy link
Owner

OK, here is the graphical test program. I took the file test.nim from the cdt package, had to patch

$ diff ~/.nimble/pkgs/cdt-0.1.1/cdt/dt.nim ~/cdt/src/cdt/dt.nim 
200c200
<   LocPos* = enum
---
>   LocPos = enum
203c203
< proc locatePoint*(this: DelaunayTriangulation; p: Vector; startEdge: Edge = nil): (Edge, LocPos) =
---
> proc locatePoint(this: DelaunayTriangulation; p: Vector; startEdge: Edge = nil): (Edge, LocPos) =

$ diff ~/.nimble/pkgs/cdt-0.1.1/cdt/drawingarea.nim ~/cdt/src/cdt/drawingarea.nim 
109c109
<   return gdk.Event_Stop
---
>   return SignalEventStopPropagation

as gintro has changed a bit in the last years. That drawingarea is based on gintro/examples/gtk3/drawingarea.nim. We should compare both for latest fixes in gintro/examples/gtk3/drawingarea.nim -- but well it seems to work.

Then I added the pointLocate() function. I store the coordinates of a random point, the located edge and the opposite point in global variables and draw them. Point to locate in red, located edge end points in blue, opposite point in green.

When you run the program (after nimble install gintro) a window appears, you can click into the window to add more constrained edges. Your desired points should be displayed in colors. When you continue clicking, constrained edge are removed again.

I still don't know how your interpolation is intended to work, I guess you have 2D points with a Z value? Well using more neightbored points for the interpolation should be not that difficult, the edges module should provide the functions which you need to find all the other points.

For the next version of the cdt module I should export the pointLocate() function. But the real task would be to write documentation, and that is some work.

import cdt/[dt, vectors, edges, types]
import gintro/cairo
import cdt/drawingarea
import random, parseutils, OS, math, tables

const
  initVector = initVector3

proc main =
  if paramCount() == 1:
    var r: int
    discard parseInt(paramStr(1), r)
    echo r
    randomize(r)
  else:
    randomize()
    let r = rand(1000)
    echo r
    randomize(r)

  proc extents(): (float, float, float, float) =
    (-10.0, -10.0, 180.0, 120.0)

  let
    minX = 0.0
    minY = 0.0
    maxX = 160.0
    maxY = 100.0

  var dt: DelaunayTriangulation = initDelaunayTriangulation(initVector(minX, minY), initVector(maxX, maxY))
  
  var el: Edge
  var pl: Vector
  var pOpp: Vertex
  var lp: LocPos
  
  proc myrand: array[2, Vector] =
    result[0].x = rand(130.0) + 15
    result[0].y = rand(70.0) + 15
    result[1].x = result[0].x + rand(20.0) - 10
    result[1].y = result[0].y + rand(20.0) - 10

  proc update =
    var i {.global.}: int
    if i == 100:
      for j in 0 .. 25:
        discard dt.insert(myrand()[0])
        discard dt.insert(myrand()[0])
        discard dt.insert(myrand())
      return

    case i:
      of 0: discard dt.insert(Vector(x: 10, y: 10), Vector(x: 20, y: 10), Vector(x: 20, y: 90), Vector(x: 10, y: 90), Vector(x: 10, y: 10))
      of 1: discard dt.insert(Vector(x: 30, y: 10), Vector(x: 40, y: 10), Vector(x: 40, y: 90), Vector(x: 30, y: 90), Vector(x: 30, y: 10))
      of 2: discard dt.insert(Vector(x: 20, y: 30), Vector(x: 30, y: 90))
      of 3: discard dt.insert(Vector(x: 20, y: 10), Vector(x: 30, y: 70))
      of 4: discard dt.insert(Vector(x: 50, y: 10), Vector(x: 60, y: 10), Vector(x: 60, y: 90), Vector(x: 50, y: 90), Vector(x: 50, y: 10))
      of 5: discard dt.insert(Vector(x: 80, y: 10), Vector(x: 90, y: 10), Vector(x: 80, y: 90), Vector(x: 70, y: 90), Vector(x: 80, y: 10))
      of 6: discard dt.insert(Vector(x: 90, y: 10), Vector(x: 100, y: 10), Vector(x: 110, y: 90), Vector(x: 100, y: 90), Vector(x: 90, y: 10))
      of 7: discard dt.insert(Vector(x: 120, y: 10), Vector(x: 130, y: 10), Vector(x: 120, y: 90), Vector(x: 110, y: 90), Vector(x: 120, y: 10))
      of 8: discard dt.insert(Vector(x: 130, y: 10), Vector(x: 140, y: 10), Vector(x: 150, y: 90), Vector(x: 140, y: 90), Vector(x: 130, y: 10))
      else: discard
    inc(i)
    if i == 10: i = -i + 1

    if i < 0:
      if not dt.removeConstraint(ConstraintID(10 + i)):
        i = 100
    block: # test the locate point, storing results in global variables
      pl = myrand()[0]
      (el, lp) = dt.locatePoint(pl)
      pOpp = el.onext.dest
      
  proc draw(cr: cairo.Context) =
    cr.setSource(1, 1, 1) 
    cr.paint
    cr.setSource(0, 0, 0, 0.5)
    cr.setLineWidth(0.3)
    cr.setLineCap(LineCap.round)
    for e in unconstrainedEdges(dt):
      cr.moveTo(e.org.point[0], e.org.point[1])
      cr.lineTo(e.dest.point[0], e.dest.point[1])
    cr.stroke
    cr.setLineWidth(0.8)
    cr.setSource(1, 0, 0, 0.6)
    for e in constrainedEdges(dt):
      cr.moveTo(e.org.point[0], e.org.point[1])
      cr.lineTo(e.dest.point[0], e.dest.point[1])
    cr.stroke
    cr.setSource(0, 0, 0, 0.5)
    for v in values(dt.subdivision.vertices):
      cr.newSubPath
      cr.arc(v.point[0], v.point[1], 1, 0, math.Tau)
    cr.stroke
    block: # draw the located points
      if el == nil: return # initially undefined!
      cr.newPath
      cr.setSource(1, 0, 0)
      cr.arc(pl.x, pl.y, 3, 0, math.Tau)
      cr.stroke # red circ for point to locate
      cr.newPath
      cr.setSource(0, 0, 1)
      cr.arc(el.org.point[0], el.org.point[1], 3, 0, math.Tau)
      cr.stroke # blue circ for points of located edge
      cr.newPath
      #cr.setSource(0, 0, 1)
      cr.arc(el.dest.point[0], el.dest.point[1], 3, 0, math.Tau)
      cr.stroke # other point of located edge
      cr.newPath
      cr.setSource(0, 1, 0)
      cr.arc(pOpp.point[0], pOpp.point[1], 3, 0, math.Tau)
      cr.stroke # green circ for point opposite to located edge
    
  var data: PDA_Data
  data.draw = draw
  data.update = update
  data.extents = extents
  data.windowSize = (800, 600)
  newDisplay(data)

main()

@HugoGranstrom
Copy link
Author

As p has a positive orientation to the returned edge, it seems to make sense to use oNext() to get the next edge with same origin. dest() should then be the opposite point you are locking for. I did a short test without graphical output to proof that, and it seems to make some sense. To do a better verify one can sketch some points on a piece of paper, enter that points then into the code and then do the point locate and verify it on the paper. Or one can graphical output on the computer, I have provided a test program with cairo graphics output in the cdt package, but have to remember how it works.

Ok, that cleared up a lot of my confusion, and judging from the graphical demo (thank goodness I installed Linux on my old laptop last week, couldn't get GIntro installed in Windows) it seems to be correct :D

I think and hope not really misleading, but I easily get confused.

I think it was just me and has who were confused. I'm not used to the cdt-jargon yet, so the meaning of "next" was what confused me. But I guess "next" always means counter-clockwise rotation around one of the points? e.dNext would give us an edge that isn't part of the triangle of interest in other words.

When you run the program (after nimble install gintro) a window appears, you can click into the window to add more constrained edges. Your desired points should be displayed in colors. When you continue clicking, constrained edge are removed again.

Very nice actually, might reuse this example when I want to practice edge-jumping for my later usecases :)

I still don't know how your interpolation is intended to work, I guess you have 2D points with a Z value? Well using more neightbored points for the interpolation should be not that difficult, the edges module should provide the functions which you need to find all the other points.

Yes exactly, I'll save the 2d points in the DT and the z-values in a Table or something. I'll experiment my way forward from here on until I get it working. It can't be harder than writing a CDT library at least ;)

For the next version of the cdt module I should export the pointLocate() function. But the real task would be to write documentation, and that is some work.

Yeah, some concrete examples would definitely be nice to have. But at the same time, as you said, you shouldn't put too much time into a project that's not very used unless you really want to. So slow and steady is the way.

Many thanks for your help Stefan!😄🌟

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
3 participants
@bassmit @HugoGranstrom @StefanSalewski