Essays/Hidden Line Removal
Hidden Line Removal
Hidden Line Removal is used in wireframe graphics, to delete lines that are hidden by nearer polygons.
Here is a J verb to perform hidden line removal. The x operand is a table of triangles in 3D (shape nx3x3), and the y operand is an edge consisting of 2 points in 3D (shape 2x3). The result is a table containing the portions of the edge y that remain visible. There may be 0, 1, or more such portions.
This code assumes that the (x,y) values are the screen positions of the lines and triangles. The z values are for occlusion calculation only and will be either the orthographic z-value or the perspective pseudo-z; in either case, the triangle is planar in xyz.
Frontfacing triangles should have CCW winding in (x,y). Higher z values occlude lower z values.
usedtocull =: 1 : 'u#]'
ifany =: ^:(*@#@])
NB. y is a line, start and end points 2x3
NB. x is a nx3x3 array of triangles, frontfaces in CCW order
NB. Result is a table of visible 2-D segments for the line, with Z component removed
NB. (2 3 3 $0 0 1 1 0 1 0 1 1 1 0 1 1 1 1 0 1 1) visiblesegments 2 3 $ 0 0 0 2 2 3
visiblesegments =: 4 : 0"3 2
line2d =. 2 {."1 line =. y
NB. Discard triangles that do not intersect the bounding square for the line
NB. Calculate bounding square x,y, min,:max
NB. Exclusion if min>max or (-max)>(-min)
NB. If any exclusion, discard the tri
qtri =. line ([: -. ((+./@,"2)@:(>"2&:(1 _1&*"2)) |.)~&:((<./ ,: >./)"2@:(2&{."1))) usedtocull x
NB. Calculate volume to each endpoint (+ volume means the face is in front
NB. of the point, having higher z, i. e. line is occludable) nx2
epvol =. qtri -/ . *@:(-"1)"2 1/ line
NB. Keep faces that are in front of at least one endpoint (or have 0 volume, which
NB. is most likely the face that contributed the line)
ftri =. (+./"1 epvol > 0) # qtri
epvol =. (+./"1 epvol > 0) # epvol
NB. Create an interval of applicability for each face. This is [0,crosspoint)
NB. or (crosspoint,0] for interpenetrations, or [0,1] for faces in front of
NB. the edge
NB. 0,1 if both +
NB. 0,cross if +,-
NB. cross,1 if -,+
applintvl =. 0 1 +"1 (0&<. % -/"1) epvol
NB. For each triangle, calculate 2-D area from each point to each edge. nx3x2
NB. The sign will be + if the point is on the inside of the edge
areas =. 0 2 1 |: ftri (-/ . *)@:(0 1 2&(}.@|."0 _)"2)@:(-"1"2 1"2)&:(2&{."1) line
NB. Convert the areas to intervals in (0,1) giving the part of the line that
NB. is inside the edge of the triangle nx3x2
edgeintvls =. 0 1 +"1 (0&<. % -/"1) areas
NB. Take the intersection of the intervals for each triangle, to give the
NB. interval of occlusion for the triangle
NB. Intersect the occlusion interval with the applicability interval to give
NB. the total occluded interval for each triangle nx2
occludeintvl =. (>./&.:(1 _1&*"1))"2 edgeintvls ,"2 1 applintvl
NB. Union the occluded intervals for all triangles to get the total occluded area.
NB. We go through pairs in ascending order of start point. Result at each point
NB. is, if overlap, lower start,larger end; if no overlap, keep the second intvl
NB. Then discard intervals that share a start point
occludeintvl =. \:~ </"1 usedtocull occludeintvl NB. Discard empty or invalid intervals
occludeintvl =. ({./.~ {."1) ({.@[ , >./&{:)^:({:@[ >: {.@])~/\. ifany occludeintvl
NB. Complement the occluded areas to produce final result: list of visible segments, in 2D
({. line) +"1 (-~/line) */~ (0.001 < -~/"1) usedtocull _2 ]\ 0 , (,occludeintvl) , 1
)
This verb can be used to create a hidden-surface-removed plot of a parametric surface, as in the following example.
load 'plot'
endtoend =: 1 : ';@:(<@u)'
createlines =: 3 : 0
NB. Create the polygons for the object, in uv space
uvgrid =. (,. {."_1) (, {.) (2p1 * (%~ i.) 60) ,"0/ (2p1 * (%~ i.) 20)
NB. Apply the function to get mesh in world space
worldgrid =. ((cos@[ * 2 + cos@]) , (sin@[ * 2 + sin@]) , (3 * sin@[ * cos@]))/"1 uvgrid
NB. Calculate viewing matrix (world-to-view)
NB. Apply azimuth first, then elevation. nxmx3
azmtx =. ((cos,sin,0:) , (-@sin,cos,0:) ,: (0,0,1:)) 2p1*214%360
elmtx =. ((cos,0,-@sin) , (0,1,0:) ,: (sin,0,cos)) 2p1*220%360
viewgrid =. (elmtx mp azmtx)&mp&.|:"2 worldgrid
NB. Create triangles (CCW winding). nxmx2x3x3
trimesh =. 2 2 3 ((0 1 3,:3 2 0) { ,/);._3 viewgrid
NB. Cull backfaces. nx3x3
fronttris =. (#~ 0 <: -/ . *@:(}. -"1 {.)"2@:(0 1&{"1)) ,/^:3 trimesh
NB. Create table of lines. nx2x3
trilines =. ~. /:~"2 ,/ (2 ]\ (,{.))"2 fronttris
NB. Produce the visible portion of each line
vislines =. fronttris visiblesegments endtoend trilines
NB. Truncate to 2D
2 {."1 vislines
)
'type line;color 0 0 0' plot ~. /:~"1 ,/^:(0 >. 2 -~ #@$) j./"1 createlines''
About 80% of the time in this example is taken calculating qtri in visiblesegments, which suggests that a spatial index on the triangles would provide a performance improvement.
Contributed by Henry Rich