function inventorycube(img_top, img_left, img_right, img_dst)
  calculate world to projection matrix A
  draw_3d_quad(img_right, img_dst, color_right, A,
      (+0.5,+0.5,-0.5), (+0.5,-0.5,-0.5), (+0.5,-0.5,+0.5), (+0.5,+0.5,+0.5))
  draw_3d_quad(img_top, img_dst, color_top, A,
      (-0.5,+0.5,+0.5), (-0.5,+0.5,-0.5), (+0.5,+0.5,-0.5), (+0.5,+0.5,+0.5))
  draw_3d_quad(img_left, img_dst, color_left, A,
      (-0.5,+0.5,-0.5), (-0.5,-0.5,-0.5), (+0.5,-0.5,-0.5), (+0.5,+0.5,-0.5))
end

function draw_3d_quad(img_src, img_dst, color, TP, w00, w01, w11, w10)
  -- map world to projection space
  p00 = first 2 coordinates of (w00*A)
  p01 = first 2 coordinates of (w01*A)
  p11 = first 2 coordinates of (w11*A)
  p10 = first 2 coordinates of (w10*A)
  -- generate texture matrix
  T = 3x3 rotation/scale/translation matrix that maps
      p00 to (0, 0),
      p01 to (0, img_src:height()-1),
      p11 to (img_dst:width()-1, img_src:height()-1))
  -- generate and draw polygon
  -- note: clipping can be skipped if we know the quad lies entirely within img_dst
  pol = make_convex_polygon([p00, p01, p11, p10])
  draw_convex_polygon(img_src, img_dst, color, pol, T)
end


--[[

all indices are 0 based, index -1 denotes the last entry of a vector/list

convex polygon data structure:
- v2f[] vl;
- v2f[] vr;
properties:
- vl[0], vr[0] are topmost vertices of the polygon (smallest y coordinate);
  vl[-1], vr[-1] are bottommost vertices (biggest y coordinate);
  the following (in)equalities hold:
    vl[0].y == vr[0].y
    vl[0].x <= vr[0].x
    vl[-1].y == vr[-1].y
    vl[-1].x <= vr[-1].x
- vl describes the left boundary of the polygon (counterclockwise);
  vertices in vl bound the polygon from the left
- vr describes the right boundary of the polygon (clockwise);
  vertices in vr bound the polygon from the right
- there are no horizontal line segments descibed by consecutive vertices in vl or vr;
  the polygon has a horizontal top boundary if and only if vl[0].x < vr[0].x;
  the polygon has a horizontal bottom boundary if and only if vl[-1].x < vr[-1].x

]]--


function make_convex_polygon(l)
  -- takes a list l containing the vertices of a convex polygon
  -- l must contain the vertices in order, but does not need to start at the topmost vertex

  i = index of element of l such that its y coordinate is minimal
  topmost = l[i]
  rest = l[i+1 .. -1] + l[0 .. i-1]
  -- in case vertices were not given in CCW order, reverse rest
  if angle(rest[0], topmost, rest[-1]) < 0 then
    rest = rest:reverse()
  end
  j = first index of element of rest such that rest[k+1].y < rest[k].y
  pol = {
    vl = [topmost] + rest[0 .. j],
    vr = [topmost] + rest[j .. -1]:reverse()]
  }

  -- remove horizontal line segments from the top of the polygon
  while #pol.vl >= 2 and pol.vl[0].y == pol.vl[1].y do
    pol.vl = pol.vl[1 .. -1]
  end
  while #pol.vr >= 2 and pol.vr[0].y == pol.vr[1].y do
    pol.vr = pol.vr[1 .. -1]
  end

  -- remove horizontal line segments from the bottom of the polygon
  while #pol.vl >= 2 and pol.vl[-2].y == pol.vl[-1].y do
    pol.vl = pol.vl[0 .. -2]
  end
  while #pol.vr >= 2 and pol.vr[-2].y == pol.vr[-1].y do
    pol.vr = pol.vr[0 .. -2]
  end

  return pol
end

function draw_convex_polygon(img_src, img_dst, color, pol, T)
  -- draw convex polygon pol onto img_dst using source pixels from img_src
  -- multiplied by constant color (used to simulate lighting)
  -- T maps coordinates in img_dst to coordinates in img_src

  clip_xmin = 0
  clip_xmax = img_dst:width() - 1
  clip_ymin = 0
  clip_ymax = img_dst:height() - 1

  miny = floor(pol.vl[0].y)
  maxy = ceil(pol.vl[-1].y)
  pl = pol.vl[0]
  pr = pol.vr[0]
  jl = 1
  jr = 1

  while jl < #pol.vl and jr < #pol.vl do
    -- find the next trapezoid
    if pol.vl[jl].y <= pol.vl[jr].y then
      pl2 = pol.vl[jl]
      jl2 = jl + 1
      pr2 = intersection of line segment pol.vr[jr-1]..pol.vr[jr] with line y = pl2.y
      jr2 = jr
    else
      pr2 = pol.vr[jr]
      jr2 = jr + 1
      pl2 = intersection of line segment pol.vl[jl-1]..pol.vl[jl] with line y = pr2.y
      jl2 = jl
    end

    -- note that pl.y == pr.y and pl2.y == pr2.y
    -- draw trapezoid bounded by the corners pl, pl2, pr2, pr
    ymin = max(floor(pl.y), clip_ymin)
    ymax = min(ceil(pl2.y), clip_ymax)
    for y = ymin to ymax do
      xmin = max(pl.x + (pl2.x - pl.x) * (y - pl.y) / (pl2.y - pl.y), clip_xmin)
      xmax = min(pr.x + (pr2.x - pr.x) * (y - pr.y) / (pr2.y - pr.y), clip_xmax)
      for x = xmin to xmax do
        img_dst[x, y] = img_src[[x, y]*T] * color
      end
    end

    -- update clipping rect to avoid drawing to a scanline twice
    clip_ymin = ymax + 1

    -- update loop variables
    pl = pl2
    pr = pr2
    jl = jl2
    jr = jr2
  end
end