docs: def projection_matrix(point, normal, direction=None,
                        perspective=None, pseudo=False):
      """Return matrix to project onto plane defined by point and normal.
  
      Using either perspective point, projection direction, or none of both.
  
      If pseudo is True, perspective projections will preserve relative depth
      such that Perspective = dot(Orthogonal, PseudoPerspective).
  
      >>> P = projection_matrix((0, 0, 0), (1, 0, 0))
      >>> numpy.allclose(P[1:, 1:], numpy.identity(4)[1:, 1:])
      True
      >>> point = numpy.random.random(3) - 0.5
      >>> normal = numpy.random.random(3) - 0.5
      >>> direct = numpy.random.random(3) - 0.5
      >>> persp = numpy.random.random(3) - 0.5
      >>> P0 = projection_matrix(point, normal)
      >>> P1 = projection_matrix(point, normal, direction=direct)
      >>> P2 = projection_matrix(point, normal, perspective=persp)
      >>> P3 = projection_matrix(point, normal, perspective=persp, pseudo=True)
      >>> is_same_transform(P2, numpy.dot(P0, P3))
      True
      >>> P = projection_matrix((3, 0, 0), (1, 1, 0), (1, 0, 0))
      >>> v0 = (numpy.random.rand(4, 5) - 0.5) * 20.0
      >>> v0[3] = 1.0
      >>> v1 = numpy.dot(P, v0)
      >>> numpy.allclose(v1[1], v0[1])
      True
      >>> numpy.allclose(v1[0], 3.0-v1[1])
      True
  
      """
      M = numpy.identity(4)
      point = numpy.array(point[:3], dtype=numpy.float64, copy=False)
      normal = unit_vector(normal[:3])
      if perspective is not None:
          # perspective projection
          perspective = numpy.array(perspective[:3], dtype=numpy.float64,
                                    copy=False)
          M[0, 0] = M[1, 1] = M[2, 2] = numpy.dot(perspective-point, normal)
          M[:3, :3] -= numpy.outer(perspective, normal)
          if pseudo:
              # preserve relative depth
              M[:3, :3] -= numpy.outer(normal, normal)
              M[:3, 3] = numpy.dot(point, normal) * (perspective+normal)
          else:
              M[:3, 3] = numpy.dot(point, normal) * perspective
          M[3, :3] = -normal
          M[3, 3] = numpy.dot(perspective, normal)
      elif direction is not None:
          # parallel projection
          direction = numpy.array(direction[:3], dtype=numpy.float64, copy=False)
          scale = numpy.dot(direction, normal)
          M[:3, :3] -= numpy.outer(direction, normal) / scale
          M[:3, 3] = direction * (numpy.dot(point, normal) / scale)
      else:
          # orthogonal projection
          M[:3, :3] -= numpy.outer(normal, normal)
          M[:3, 3] = numpy.dot(point, normal) * normal
      return M