How lit is this room? pt. 1

Question

Related to this question.

A room is defined to be a (not necessarily convex) non-intersecting polygon, expressed as an ordered list of 2-dimensional coordinates. A sufficiently bright lightbulb is placed at a specific point inside the room, and emits light in every direction. Your task is to find the total illuminated area of the room. You may take in input in any reasonable format. The points on the polygon/room as well as the coordinates of the light source are rational numbers. They can be taken in clockwise or counterclockwise, either format is fine. The test case in the problem is given counterclockwise.

The following picture shows two example rooms, where the purplish dot represents the light source, and the shaded region represents the illuminated area.

Test case:

(1/2, 18)
(1,3)
(5,1/2)
(7,5)
(12,7)
(16,3)
(15,11)
(8,19)
(3,7)
Light source located at (5,8)
Answer: 815523/6710 ≈ 121.538

Here is a graphical depiction of the solution to that test case. The two points that define the solution that are not on the original polygon are (55/61, 363/61) and (856/55, 357/55).

This formula may be helpful in calculating the area.

Since this is code-golf, the shortest code in bytes wins.

Answer 1

Python 3, 388 398 408 409 415 417 493 bytes

---

To make it more accurate, increase n

from random import*
u=uniform
c=lambda A,B,C:(C[1]-A[1])*(B[0]-A[0])>(B[1]-A[1])*(C[0]-A[0])
I=lambda A,B,C,D:c(A,C,D)!=c(B,C,D)and c(A,B,C)!=c(A,B,D)
def a(l,v,n=9**6,s=0):
 g=lambda i:(min(x[i]for x in v),max(x[i]for x in v))
 for _ in'x'*n:
  h=((u(*g(0)),u(*g(1))),l);s+=any([I(*f,*h)for f in list(zip(v,v[1:]+[v[0]]))])^1
 return(abs(g(0)[0]-g(0)[1])*abs(g(1)[0]-g(1)[1]))*float(s/n)

Basic Monte-Carlo approach. Steps listed below.

1. Find x and y ranges that the shape occupies.
2. Create a list of edges created by the vertices
3. Iterate a large number of times (the more the better)
4. Create a random point (j,k) inside the x,y range.
5. Check if any of the edges intercept with the line-segment created by the light and the random point. If any of the edges intercept, increment the variable s
6. Divide s by the total numbers, then multiply by the total range area.

Ungolfed version:

import random

def ccw(A,B,C):
    return (C[1]-A[1])*(B[0]-A[0]) > (B[1]-A[1])*(C[0]-A[0])

def intersect(A,B,C,D):
    return ccw(A,C,D) != ccw(B,C,D) and ccw(A,B,C) != ccw(A,B,D)

def lit_area(light, vertices):
    # points: list of points
    # i     : x => i=0
    #       : y => i=1
    get_range = lambda i: (min(x[i] for x in vertices), max(x[i] for x in vertices))
    xr = abs(get_range(0)[0] - get_range(0)[1])
    yr = abs(get_range(1)[0] - get_range(1)[1])

    edges = list(zip(vertices, vertices[1:] + [vertices[0]]))

    num_sims = 1000000

    num_successes = 0
    for _ in range(num_sims):
        guess_x = random.uniform(*get_range(0))
        guess_y = random.uniform(*get_range(1))

        light_guess_line = ((guess_x, guess_y), light)

        if not any([intersect(*e, *light_guess_line) for e in edges]):
            num_successes += 1
    return float(num_successes / num_sims) * (xr * yr)


if __name__ == "__main__":
    points = [
    (1/2, 18),
    (1,3),
    (5,1/2),
    (7,5),
    (12,7),
    (16,3),
    (15,11),
    (8,19),
    (3,7)
    ]
    light_source = (5,8)
    print("Area lit by light: %f"% lit_area(light_source, points))

Try it online!

Credit for line intersection algorithm

Also, credit to all the helpful commenters on how to golf this even further.

Answer 2

Haskell, 559 _{^{618 632}} bytes

r(a:b)=b++[a]
s=zip<*>r
(?)a=sum.zipWith(*)a
o(a,b)=r a?b-a?r b
(a,b)!(c,d)=(c-a,d-b)
(a,b)#(c,d)=a*d-b*c
x i a@(e,f)b j c d|let k@(g,h)=a!b;l=c!d;m=c!a;n=l#k;o=m#l/n;p=m#k/n;q|i>0=o<0||o>1|let=o<=0||o>=1;r|n==0||q||p<0||p*j>1=[]|let=[(e+o*g,f+o*h)]=r
(a&b)(c:e@(d:_))|let(f,g)=span(/=d)b;h=zip f$r$f++[d]=concat[[k,l]|(i,j)<-h,[[k],[l]]<-[x 1 i j 0 a<$>[c,d]],and[x 0 m n 1 a o==[]|o<-[k,l],(m,n)<-h,(m,n)/=(i,j)]]++(a&g)e
(_&_)_=[]
z a b=sum[o$unzip[c,a,d]|e@(f:_)<-[[c|c<-b,and[all(==c)$x 1 d e 1 a c|(d,e)<-s b]]],(c,d)<-s$a&until((f==).head)r b$e++[f]]/2

Exact solution (barring bugs). Haskell has built-in exact rational arithmetic. Try it online!

Note that this gives 815523/6710, not 814643/6710, for the example room, and the first wall intersection is calculated as (55/61, 363/61). I'm fairly sure this is correct because the Monte Carlo entry (slowly) converges to the same result.

Legend:

z light roomPoints
    -- Main function, returns lit area.
    -- Compute list of visible corners in the room, then calls (&).
(&) light roomPoints' visibleCorners
    -- Compute visibility polygon. visibleCorners is the subset of points
    -- that are visible from the light. The first point of roomPoints'
    -- must coincide with the first visibleCorner.
x pEndpoints p1 p2 qSegment q1 q2
    -- Intersect line segments (p1, p2) and (q1, q2).
    -- If pEndpoints, exclude endpoints p1, p2.
    -- If not qSegment, allow intersection to extend past q2 (i.e. raycast).
r   -- Rotate list by one, used to construct closed loops etc.
s   -- Construct closed loop
(!) -- Vector between two points
(?) -- Dot product
(#) -- Cross product
o   -- Polygon area

Bonus: Gloss GUI for testing. Click next to points to move them.

import qualified Graphics.Gloss as G
import qualified Graphics.Gloss.Interface.IO.Interact as GI

solnPoly a b|let c@(d:_)=[c|c<-b,and[all(==c)$x 1 d e 1 a c|(d,e)<-s b]]=a&until((d==).head)r b$c++[d]
solnArea = z

main =
  let fromRatP (x, y) = (fromRational x, fromRational y)
      displayScale = 10
      scalePoints = G.scale (fromInteger displayScale) (fromInteger displayScale)
      displayMode = G.InWindow "" (512, 512) (0, 0)
      drawBasePoly pointSz ps =
          mconcat $ G.lineLoop ps :
                    [G.translate x y (G.circleSolid pointSz) | (x, y) <- ps]
      drawVisPolyOf light ps =
          G.color G.blue $ drawBasePoly 0.2 $ map fromRatP $ solnPoly light ps
      drawLight (x, y) =
          G.translate x y $
          G.color G.yellow (G.circleSolid 0.5) <> G.circle 0.5
      draw (light, ps) =
          mconcat [
              scalePoints $ drawLight (fromRatP light),
              scalePoints $ drawBasePoly 0.4 (map fromRatP ps),
              scalePoints $ drawVisPolyOf light ps,
              G.translate (-200) (-50) $ G.scale 0.2 0.2 $
                G.color G.blue $ G.text $ "Lit area: " ++ show (solnArea light ps)
          ]
      event (GI.EventKey (GI.MouseButton GI.LeftButton) GI.Down _ (curx_, cury_)) (light, ps) =
          let dist (x,y) (x',y') = (x'-x)^2 + (y'-y)^2
              curx = curx_ / fromInteger displayScale
              cury = cury_ / fromInteger displayScale
              cursorR = (fromInteger$round curx, fromInteger$round cury)
              maxDist = 3
              snapAmount = 1
              (d, i) = minimum [(dist p cursorR, i) | (p, i) <- zip (light : ps) [0..]]
              snapTo n a = fromInteger$n*round(a/fromInteger n)
              snapCursor = (snapTo snapAmount curx, snapTo snapAmount cury)
              light' | i == 0 && d < maxDist^2 = snapCursor
                     | otherwise = light
              ps' | i > 0 && d < maxDist^2 = take (i-1) ps ++ [snapCursor] ++ drop i ps
                  | otherwise = ps
          in (light', ps')
      event _ state = state
      state0 =
        ((2, 2), [(0, 0), (10, 0), (10, 5), (20, 0), (20, 20), (15, 5),
                  (10, 10), (6, 10), (10, 12), (0, 12), (4, 10), (0, 10)])
  in G.play displayMode G.white 60
            state0
            draw
            event
            (\_ -> id)

Answer 3

APL+WIN

This is an ungolfed version of this interesting challenge which I offer to demonstrate my logic. My ancient version of APL+WIN is not well suited to golfing nested control structures. More modern APLs could well do better - challenge?

If readers validate the logic I will have a go at golfing this solution. If the logic is wrong I will simply delete.

r←b Room v

⍝Separate x and y coordinates of vertices               
x←v[;1] ⋄ y←v[;2]

⍝Intercept and slope of each line segment and ray through each vertex
s←(y,¨1⌽y)⌹¨(1E¯9+1,[1.1]¨x,¨1⌽1E¯9+x)
l←(y,¨b[2])⌹¨(1E¯9+1,[1.1]¨x,¨b[1]+1E¯9)                          

⍝Coordinates of vertices
x←x,¨1⌽x ⋄ y←y,¨1⌽y                                                  

⍝Initialise intersection matrix
r←((⍴s),0)⍴0

⍝Evaluate intersection matrix 
:for i :in ⍳⍴l 
    t←0⍴0
    :for j :in ⍳⍴s
        t←t,⍎1⍕∊((↑∊l[i])-↑∊s[j])÷((1↓∊s[j])-1↓∊l[i]) 
    :endfor
    z←r←r,t
:endfor 

⍝Identify x coordinates of valid intersections in the direction of the ray
:for i :in ⍳⍴l 
    t←(r[i;i])
    :for j :in ⍳⍴s
        :if t<b[1] 
            r[j;i]←r[j;i]×(r[j;i]<t)^r[j;i]>⌊/∊x[i]
        :else
            r[j;i]←r[j;i]×(r[j;i]>t)^r[j;i]<⌈/∊x[i]
        :endif
    :endfor
 :endfor

⍝Identify the edges intersected
e←(+/r≠0)/⍳⍴l 

⍝Intersection x coordinates
cx←(+/r)[e]

⍝Intersection y coordinates
cy←⍎1⍕+/¨(s[e])× 1,¨(+/r)[e]

⍝Replace original coordinates that are in shadow
x[e]←(1↓¨x[e]),¨cx
y[e]←(1↓¨y[e]),¨cy

⍝Calculate lit area
r←+/.5×(|-/¨x)×|-/¨y                                               

Answer 4

R, 296 255 bytes

function(s,l,u=cbind(s,s[,1]),d=t(diff(t(u))),q=l-u,r=apply(s,1,range),g=c(diff(r)))mean(replicate(1e6,!any((q[i<-1:ncol(s)*2]*(p=runif(2)*g+r[1,]-u)[j<-i-1]>p[i]*q[j])!=(q[i+2]*p[i+1]>q[i+1]*p[i+2])&(p[i]*d[j]>p[j]*d[i])!=(q[i]*d[j]>q[j]*d[i]))))*prod(g)

Try it online!

This is a further reduced version of the Python answer. The core Monte Carlo method is the same, but I rearranged some of the functions to make them shorter. In my first iteration, I'd been over-aggressive in rearrangement, and I then realised I could optimise both length and speed by returning to a version of the intersection algorithm closer to the python.

Here's an ungolfed version that also plots the results:

find_lit_ungolf <- function(shape, light, plot = TRUE) {
  lines <- cbind(shape ,shape[,1])
  diffs <- t(diff(t(lines)))
  light_minus_lines <- light - lines
  shape_range <- apply(s,1,range)
  shape_range_diff <- c(diff(shape_range))
  successes <- t(replicate(
    n = 1e5,
    {
      random_point <- runif(2) * shape_range_diff + shape_range[1, ]
      random_minus_lines <- random_point - lines
      q <- light_minus_lines
      p <- random_minus_lines
      d <- diffs
      i <- 1:ncol(s)*2
      success <-
        !any((q[i]*p[i-1]>p[i]*q[i-1])!=(q[i+2]*p[i+1]>q[i+1]*p[i+2])&(p[i]*d[i-1]>p[i-1]*d[i])!=(q[i]*d[i-1]>q[i-1]*d[i]))
      c(random_point, success)
    }))
  colnames(successes) <- c("x", "y", "success")
  if (plot) {
    shape <- t(shape)
    colnames(shape) <- c("x", "y")
    print(ggplot(as_tibble(successes), aes(x, y)) +
      geom_point(aes(colour = factor(success)), alpha = 0.3) +
      geom_polygon(data = as_tibble(shape), alpha = 0.2) +
      annotate("point", light[1], light[2], col = "yellow"))
  }
  mean(successes[, 3]) * prod(shape_range_diff)
}
find_lit_ungolf(s, l)