{+doc
+NOMBRE: IntPoint
+CREACION: 31.1.98
+AUTORES:rch
+REGISTRO:
+TIPO: Unit PASCAL
+PROPOSITO:
   Algoritmo que determina si un punto es interior o no a un
   polígono definido como un vector de complexs.

+PROYECTO:

+REVISION:
+AUTOR:
+DESCRIPCION:
-doc}

{ La funci¢n InternalPoint retorna TRUE si el punto Q es interior
al Poligono que forman los puntos del vector X pasado como par metro }
unit IntPoint;

{$MODE Delphi}

interface

uses
  Math, xMatDefs, AlgebraC, MatCpx;

// recibe un puntero al primer Punto de un Poligono y la cantidad NPuntos del polígono
// que deben estar en direcciones contiguas de memoria.
// La función determina si el punto Q es interior al Poligono o no.
// El método utilizado para determinar si el punto es interior consiste en
// contar las veces que los segmentos del Polígono cortan la semirecta que
//  x = Q.r e y >= Q.i
// Si la cantidad de cortes es PAR el punto es exterior. Si la cantidad es impar
// el punto es interior.
function InternalPoint(Q: NComplex; pPrimerPunto: PNComplex; NPuntos: integer): boolean;
  overload;

// lo mismo que la anterior, pero sobre un Poligono pasado como un array dinámico de complejos
function InternalPoint(Q: NComplex; const poligono: TVarArrayOfComplex): boolean;
  overload;

// productoVectorial( u, v )
// res= ro_u * ro_v * sin( alfa_u - alfa_v ) = img( u * cc( v ) )
function productoVectorial(const u, v: NComplex): NReal;

// angulo_uv( u, v )
// res = fase( u * cc( v ) )
function angulo_uv(const u, v: NComplex): NReal;

// dado un Punto Q y un Poligono, retorna los dos vértices Va, Vb, del Poligono que
// forman el mayor ángulo con Q. Si Q es exterior al polígono, el ángulo formado
// en sentido anti-horario por las semirectas de origen Q y que pasan por Va y Vb
// contienen al polígono.
// Se asume que el polígono es convexo y que NPuntos >= 2
procedure AnguloAlConvexo(var kA, kB: integer;
  // vertices derecho e izquierdo del ángulo al Poligono ( 0..NPUntos-1)
  var alfa_P0_A, alfa_P0_B: NReal; const Q: NComplex; pPrimerPunto: PNComplex;
  NPuntos: integer // el poligono Convexo
  ); overload;

procedure AnguloAlConvexo(var kA, kB: integer;
  // vertices derecho e izquierdo del ángulo al Poligono
  var alfa_P0_A, alfa_P0_B: NReal; const Q: NComplex;
  const poligono: TVarArrayOfComplex); overload;


// AmpliarConvexoAlPunto
// Amplia un convexo dado para incluir el punto dado
// El resultado es un nuevo convexo que incluye al punto Q.
procedure AmpliarConvexoAlPunto(var Convexo: TVarArrayOfComplex; const Q: NComplex);

// ContornoConvexo
// Dada una Nube de Complejos, retorna el Poligono que describe el contorno
// convexo más ajustado a la nube.
function ContornoConvexo(const Nube: TVarArrayOfComplex): TVarArrayOfComplex;

implementation

// res= ro_u * ro_v * sin( alfa_u - alfa_v ) = img( u * cc( v ) )
function productoVectorial(const u, v: NComplex): NReal;
begin
  Result := u.r * v.i - u.i * v.r;
end;

// angulo_uv( u, v )
// res = fase( u * cc( v ) )
function angulo_uv(const u, v: NComplex): NReal;
begin
  Result := fase(pc(u, cc(v)^)^);
end;


procedure AnguloAlConvexo(var kA, kB: integer;
  // vertices derecho e izquierdo del ángulo al Poligono ( 0..NPUntos-1)
  var alfa_P0_A, alfa_P0_B: NReal; const Q: NComplex; pPrimerPunto: PNComplex;
  NPuntos: integer // el poligono convexo
  );
var
  minAlfa, maxAlfa, alfa: NReal;
  p: PNComplex;
  QP0, QPk, QPk_ant: NComplex;
  buscando: boolean;
  k: integer;
  estado: integer; // 0 no determiné el sentido, 1 sentido determinado, 2 reversa, 3 fin
  sentido: integer; // 1 alfa crece con k, -1 alfa decrece con k

  procedure buscar_alfa_creciente_inck;
  begin
    buscando := True;
    while buscando and (k < (NPuntos - 1)) do
    begin
      Inc(p);
      Inc(k);
      QPk_ant := QPk;
      QPk := rc(p^, Q)^;
      alfa := angulo_uv(QPk, QP0);
      if alfa > maxAlfa then
        maxAlfa := alfa
      else
      if alfa = maxAlfa then
      begin
        if mod2(QPk) < mod2(QPk_ant) then
          buscando := False;
      end
      else
        buscando := False;
    end;
  end;

  procedure buscar_alfa_decreciente_inck;
  begin
    buscando := True;
    while buscando and (k < (NPuntos - 1)) do
    begin
      Inc(p);
      Inc(k);
      QPk_ant := QPk;
      QPk := rc(p^, Q)^;
      alfa := angulo_uv(QPk, QP0);
      if alfa < minAlfa then
        minAlfa := alfa
      else
      if minAlfa = alfa then
      begin
        if mod2(QPk) < mod2(QPk_ant) then
          buscando := False;
      end
      else
        buscando := False;
    end;
  end;

  procedure buscar_alfa_creciente_deck;
  begin
    buscando := True;
    while buscando and (k > 0) do
    begin
      Dec(p);
      Dec(k);
      QPk_ant := QPk;
      QPk := rc(p^, Q)^;
      alfa := angulo_uv(QPk, QP0);
      if alfa > maxAlfa then
        maxAlfa := alfa
      else
      if maxAlfa = alfa then
      begin
        if mod2(QPk) < mod2(QPk_ant) then
          buscando := False;
      end
      else
        buscando := False;
    end;
  end;

  procedure buscar_alfa_decreciente_deck;
  begin
    buscando := True;
    while buscando and (k > 0) do
    begin
      Dec(p);
      Dec(k);
      QPk_ant := QPk;
      QPk := rc(p^, Q)^;
      alfa := angulo_uv(QPk, QP0);
      if alfa < minAlfa then
        minAlfa := alfa
      else
      if minAlfa = alfa then
      begin
        if mod2(QPk) < mod2(QPk_ant) then
          buscando := False;
      end
      else
        buscando := False;
    end;
  end;

begin
  minAlfa := 0;
  maxAlfa := 0;
  QP0 := rc(pPrimerPunto^, Q)^;

  // búsqueda k_Creciente
  buscando := True;
  k := 0;
  p := pPrimerPunto;
  QPk := QP0;

  // busco el primer vértice que me permita determinar el sentido
  while buscando and (k < (NPuntos - 1)) do
  begin
    Inc(p);
    Inc(k);
    QPk_ant := QPk;
    QPk := rc(p^, Q)^;
    alfa := angulo_uv(QPk, QP0);
    if (alfa > 0) and (alfa < pi) then
      sentido := 1
    else if alfa < 0 then
      sentido := -1;
  end;

  // si no pude determinar el sentido quiere decir que el polígono
  // es una recta que pasa por Q. En este sentido, debiera ser un
  // segmento. Devolvemos como puntos extremos el Primero y El ultimo;
  if buscando then
  begin
    kA := 0;
    alfa_P0_A := 0;
    kB := NPuntos - 1;
    alfa_P0_B := alfa;
    exit;
  end;

  // ahora continuamos incrementando k hasta que se invierta el sentido
  // o lleguemos al final
  if sentido > 0 then
  begin
    buscar_alfa_creciente_inck;
    kB := k;
    p := pPrimerPunto;
    Inc(p, NPuntos); // apunto el último pues voy en reversa.
    k := 0;
    buscar_alfa_decreciente_deck;
    kA := k;
  end
  else
  begin
    buscar_alfa_decreciente_inck;
    kB := k;
    p := pPrimerPunto;
    Inc(p, NPuntos); // apunto el último pues voy en reversa.
    k := 0;
    buscar_alfa_creciente_deck;
    kA := k;
  end;
end;

procedure AnguloAlConvexo(var kA, kB: integer;
  // vertices derecho e izquierdo del ángulo al Poligono
  var alfa_P0_A, alfa_P0_B: NReal; const Q: NComplex;
  const poligono: TVarArrayOfComplex);
begin
  AnguloAlConvexo(kA, kB, alfa_P0_A, alfa_P0_B, Q, @poligono[0], length(poligono));
end;

function InternalPoint(Q: NComplex; pPrimerPunto: PNComplex; NPuntos: integer): boolean;
var
  par: boolean;
  k: integer;
  amb1, amb2: integer;
  p1, p2: PNComplex;

  function intersec(var A, B: NComplex): integer;
  var
    QA, QB: NComplex;
    t1, t2: NReal;
  begin
    QA := rc(A, Q)^;
    QB := rc(B, Q)^;
    t1 := QA.r;
    t2 := QB.r;
    if productoVectorial(QA, QB) > 0 then
      if (t1 < 0) and (t2 > 0) then
        intersec := 1
      else if (t1 = 0) or (t2 = 0) then
        if (t1 < 0) or (t2 > 0) then
          intersec := 2
        else
          intersec := 0
      else
        intersec := 0
    else
    if (t1 > 0) and (t2 < 0) then
      intersec := 1
    else if (t1 = 0) or (t2 = 0) then
      if (t1 > 0) or (t2 < 0) then
        intersec := -2
      else
        intersec := 0
    else
      intersec := 0;
  end;

  procedure proc1;
  begin
    case intersec(p1^, p2^) of
      1: par := par xor True;
      2: Inc(amb1);
      -2: Inc(amb2);
      else {nada}
    end;
  end;

begin
  par := False;
  amb1 := 0;
  amb2 := 0;


  p1 := pPrimerPunto;
  p2 := p1;
  Inc(p2);
  for k := 1 to NPuntos - 1 do
  begin
    proc1;
    p1 := p2;
    Inc(p2);
  end;
  p2 := pPrimerPunto;
  proc1;
  amb1 := amb1 - amb2;
  internalPoint := par xor (amb1 <> 0);
end;


function InternalPoint(Q: NComplex; const poligono: TVarArrayOfComplex): boolean;
begin
  Result := InternalPoint(Q, @poligono[0], length(poligono));
end;


procedure AmpliarConvexoAlPunto(var Convexo: TVarArrayOfComplex; const Q: NComplex);
var
  res: TVarArrayOfComplex;
  kA, kB, k: integer;
  kTest: integer;
  flgTestTramo, eliminarTramo: integer;
  alfa_P0_Q_A, alfa_P0_Q_B: NReal;
  uAB, uAQ, uAP0: NComplex;
  //  alfa_Q, alfa_P0: NReal;
  d2AQ, d2AB, d2QB: NReal;

  k1, k2: integer;
  uAP: NComplex;
  alfa_BAP, alfa_BAQ: NReal;

begin
  if length(Convexo) < 2 then
  begin
    if length(convexo) = 0 then
    begin
      setlength(convexo, 1);
      convexo[0] := Q;
    end
    else
    begin
      uAB := rc(convexo[0], Q)^;
      if mod2(uAB) > 0 then
      begin
        setlength(Convexo, Length(Convexo) + 1);
        res[high(res)] := Q;
      end;
    end;
    exit;
  end;

  // Si no salí antes es que el convexo tiene por lo menos 3 puntos.

  if (InternalPoint(Q, Convexo)) then
    exit;

  // Si el punto no es interior, puede ser exterior o estar en el contorno.
  AnguloAlConvexo(kA, kB, alfa_P0_Q_A, alfa_P0_Q_B, Q, Convexo);
  uAB := rc(Convexo[kB], Convexo[kA])^;
  uAQ := rc(Q, Convexo[kA])^;

  // detección del caso especial en que Q pertece a la recta AB
  // para ello, alfa_P0_Q_A y alfa_P0_Q_B son CERO o PI
  // si ambos tiene el mismo valor entonces Q es externo al segmento (AB)
  // y deberá sustituir uno de ambos. Si son diferentes Q está entre A y B
  // y entoces está incluído en el segmento y no hay que hacer nada.
  if ((alfa_P0_Q_A = 0) or (alfa_P0_Q_A = pi)) and
    ((alfa_P0_Q_B = 0) or (alfa_P0_Q_B = pi)) then
  begin
    if (alfa_P0_Q_A = alfa_P0_Q_B) then
    begin
      d2AQ := mod2(uAQ);
      d2QB := mod2(rc(Convexo[kB], Q)^);
      if d2AQ < d2QB then
        convexo[kA] := Q
      else
        convexo[kB] := Q;
    end;
    exit;
  end;

  // Si llegamos hasta aquí, Q es externo al convexo (incluyendo su contorno)
  // y además no está sobre la recta AB.


  assert(kA <> kB, 'AmpliarConvexoAlPunto kA<>kB FALLO!');

  if kA < kB then
  begin
    k1 := kA;
    k2 := kB;
  end
  else
  begin
    k1 := kB;
    k2 := kA;
  end;


  if (length(convexo) = 2) then
  begin
    // simplemente agregamos el punto Q.
    setlength(convexo, 3);
    convexo[2] := Q;
    exit;
  end;


  // El Convexo está dividido en 2 tramos por el segmento P_k1 P_k2
  // Un tramo es el de los vértices k1, k1+1, k1+2 , ... k2 (Tramo 1)
  // El otro tramo es el de los vértives k2, k2+1, ..NPuntos-1, 0, 1, ... k1 (Tramo 0)
  // O sea el Tramo 0, "cierra en anillo"
  // justo sea k1 = 0 o k2 = 0.
  // De los dos tramos, hay uno que está "mas cercano a Q" y que debe ser
  // eliminado del contorno del Convexo pues es sutituído por los lados
  // (Q, P_k1) y (Q, P_k2 )

  if (k2 - k1) > 1 then // hay por lo menos un punto en el tramo superior
  begin
    // testeo sobre el tramo superior
    kTest := k1 + 1;
    flgTestTramo := 1;
  end
  else
  begin
    //testeo sobre el tramo inferior
    flgTestTramo := 0;
    if (k1 > 0) then
      kTest := 0
    else
    if k2 < high(convexo) then
      kTest := k2 + 1;
  end;

  uAP := rc(Convexo[kTest], Convexo[kA])^;
  alfa_BAP := angulo_uv(uAB, uAP);
  alfa_BAQ := angulo_uv(uAB, uAQ);

  if (sign(alfa_BAP) = sign(alfa_BAQ)) then
    eliminarTramo := flgTestTramo
  // si coinciden los ángulos el tramo a eliminar es el bajo test
  else
    eliminarTramo := 1 - flgTestTramo;
  // si no coinciden es el tramo que no está bajo test

  if eliminarTramo = 1 then
  begin
    // 0..k1 , Q , k2...
    setlength(res, (k1 + 1) + 1 + (high(convexo) - k2 + 1));
    for k := 0 to k1 do
      res[k] := convexo[k];
    res[k1 + 1] := Q;
    for k := k2 to high(convexo) do
      res[k - k2 + k1 + 2] := convexo[k];
    convexo := res;
  end
  else
  begin
    // Q, k1 ... k2
    setlength(res, 1 + (k2 - k1 + 1));
    res[0] := Q;
    for k := k1 to k2 do
      res[k - k1 + 1] := convexo[k];
    convexo := res;
  end;
end;

function ContornoConvexo(const Nube: TVarArrayOfComplex): TVarArrayOfComplex;
var
  res: TVarArrayOfComplex;
  NPuntos: integer;
  k: integer;

begin
  setlength(res, 0);
  for k := 0 to high(nube) do
    AmpliarConvexoAlPunto(res, Nube[k]);
  Result := res;
end;

end.
