Calculo del área cubierta por una hoja doblada. Visualización en geogebra: http://www.geogebratube.org/student/m52908?mobile=true

hoja_doblada.py

"""Calculo del area cubierta por una hoja doblada."""
from __future__ import division

import math
import numpy as np


def solve(w, h, x1, x2):
    """Resuelve el calculo de area cubierta de una hoja doblada.

            |---- w ----|
                    x1
            +-------/---+   -
            |      /    |   |
            |     /     |   h
            |    /      |   |
            +---/-------+   -
              x2

    Al doblarse el lado izquierdo del eje x1-x2, se obtiene el lado simetrico
    (o reflexion) respecto el eje x1-x2. Suficiente calcular el area de
    interseccion de la reflexion del lado izquierdo y el lado derecho del eje
    x1-x2.
 
    La respuesta buscada es el area total menos el area de interseccion.

    Parametros:
        w : ancho de la hoja
        h : alto de la hoja
        x1 : distancia superior del eje doblez
        x2 : distancia inferior del eje de doblez
    """
    # calculo directo
    if x1 == x2:
        # el area cubierta es la mayor region cortada por x1-x2
        return max(w - x1, w - x2) * h

    # asumimos x1 > x2 dado que el area de interseccion es la misma
    if x1 < x2:
        x1, x2 = x2, x1

    return w * h - inter_area(w, h, x1, x2)


def inter_area(w, h, x1, x2):
    """Calcula el area de interseccion de las regiones divididas por la recta
    x1-x2 al doblar la hoja por el misma recta como eje.

           c       x1     d
            +------q-----+  -
            |     /      |  |
            |    /       |  | h
            |   /        |  |
            +--o---------+  -
           a  x2          b
            |------------|
                  w
    """
    # calculamos los vectores tomando como origen (x2, 0)
    w2 = w - x2
    o = np.array([0, 0])
    a = np.array([-x2, 0])
    b = np.array([w2, 0])
    c = np.array([-x2, h])
    d = np.array([w2, h])
    q = np.array([x1 - x2, h])

    # calculamos las reflexiones de las esquinas
    c2 = reflection(c, q)
    a2 = reflection(a, q)

    # El segmento q-c2 puede tener interseccion solo con un lado de la hoja o
    # ninguno (punto interior).
    # El segmento o-a2 no tiene interseccion con la hoja dado las suposiciones
    # iniciales.
    # El segmento a2-c2 puede tener solo una interseccion con un lado o dos
    # intersecciones con ambos lados o bien ninguna interseccion.

    if c2[0] <= w2 and 0 <= c2[1] <= h:
        # Punto dentro la hoja
        r = c2
        if (a2 == o).all():
            # El area de interseccion es o-q-r
            return tri_area(o, q, r)
        else:
            # El segmento c2-a2 intersecta o-b
            s = intersect(c2, a2, o, b)
            # El area de interseccion esta definido por o-q-r-s
            return tri_area(o, q, r) + tri_area(o, r, s)
    else:
        # Punto fuera de la hoja
        # El segmento q-c2 puede o no intersectar con b-d
        r = intersect(q, c2, b, d)
        if r is not None:
            # Si intersecta, el segmento c2-a2 puede o no intersectar b-d
            s = intersect(c2, a2, b, d)
            if s is not None:
                # El segmento c2-a2 intersecta o-b
                t = intersect(c2, a2, o, b)
                # El area de interseccion es o-q-r-s-t
                return tri_area(o, q, r) + tri_area(o, r, s) + tri_area(o, s, t)
            else:
                # La region de interseccion esta contenida dentro o-q-c2-a2
                # El area de interseccion es o-q-r-b
                return tri_area(o, q, r) + tri_area(o, r, b)
        else:
            # Si no intersecta, q-c2 intersecta si o si a o-b
            r = intersect(q, c2, o, b)
            # La region de interseccion esta contenida en o-q-c2-a2
            # El area de interseccion es o-q-r
            return tri_area(o, q, r)


def reflection(v, w):
    """Calcula la reflexion de v por la recta de direccion w.

                w
        v *    /
              /
             /    * v'
            /
           o

    >>> v = np.array([1, 2]); w = np.array([3, 4])
    >>> reflection(reflection(v, w), w) == v
    array([ True,  True], dtype=bool)

    """
    u = np.dot(v, w) / np.dot(w, w) * w
    return v + 2 * (u - v)


def tri_area(x, y, z):
    """Calcula el area de un triangulo.

    >>> x = np.array([0, 0]); y = np.array([1, 0]); z = np.array([1, 1])
    >>> round(tri_area(x, y, z), 2)
    0.5
    >>> x = np.array([0, 0]); y = np.array([1, 1]); z = np.array([2, 0])
    >>> round(tri_area(x, y, z), 2)
    1.0

    """
    # Heron's formula
    a = np.linalg.norm(y - x)
    b = np.linalg.norm(z - x)
    c = np.linalg.norm(z - y)
    s = (a + b + c) / 2
    return math.sqrt(s * (s - a) * (s - b) * (s - c))


def intersect(v1, v2, w1, w2):
    """Calcula el punto de interseccion de los segmentos definidas por v1-v2 y
    w1-w2.

    >>> v1 = np.array([0, 0]); v2 = np.array([2, 2])
    >>> w1 = np.array([0, 2]); w2 = np.array([2, 0])
    >>> intersect(v1, v2, w1, w2)
    array([ 1.,  1.])

    """
    d1 = v2 - v1
    d2 = w2 - w1
    A = np.array([d1, d2]).T
    b = w1 - v1
    k, _ = np.linalg.solve(A, b)
    # nos interesa punto dentro el segmento
    if 0 <= k <= 1:
        return v1 + k * d1


if __name__ == '__main__':
    print solve(7, 100, 1, 1)
    print solve(13, 5, 10, 3)
    print solve(13, 5, 3, 10)
    print solve(100, 7, 90, 1)
    print solve(10, 10, 0, 10)