% env.w -- Envolvente Convex para MapBasic
%
% Copyright (C) 1997, Javier Goizueta <javier@@goizueta.info>
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.

% ===========================================================================
\documentclass[a4paper,oneside,english]{article}
\usepackage[english,read]{nwprog}
% ===========================================================================

\isodate

\newcommand{\ProgTitle}{Envolvente convexa}
\newcommand{\ProgAuth}{Javier Goizueta} 
\newcommand{\ProgDate}{Mayo 1997} 
\newcommand{\ProgVer}{1.0}
\newcommand{\ProgSource}{\ttfamily\bfseries tools.w}

\title{\ProgTitle}
\author{\ProgAuth} 
\date{\ProgDate}

% ===========================================================================


% Éste código fuente nuweb permite generar:
% Código:
%   env.def -- Fuente MapBasic, Declaraciones Interface
%   env.mb  -- Fuente MapBasic, Implementación
% Documentación
%   env.tex -- Fuente LaTeX, incluye EPS de env.mp
%                   (combinar env.dvi con envN.eps)
% el archivo env.mp contiene las ilustraciones de la documentación

\usepackage{latexsym}

%\usepackage[dvipdfm,colorlinks=true,linkcolor=blue]{hyperref}
%%\usepackage[html]{hyperref}

%\usepackage[dvips]{graphics}

\begin{document}


\maketitle

% opcional: - - - - - - - - - - - - 
%\tableofcontents
%\section{Introducción}
% - - - - - - - - - - - - - - - - -

Este módulo de MapBasic calcula la \emph{envolvente convexa} de una
nube de puntos (función \cd{ConvexHull}).
La envolvente convexa es un polígono convexo que contiene a todos
los puntos y cuyos vértices son puntos del conjunto; si colocamos
un clavo en cada punto y a su alrededor una goma elástica
que quede tensa, la forma que esta goma adopta es la de la
envolvente convexa.

% FIG1: nube de puntos y envolvente
\begin{figure}[!htbp]
\begin{center}
\includegraphics{env.1}
\caption{envolvente convexa}
\end{center}
\end{figure}

 Una aplicación útil de esta
función es el cálculo de ``áreas percentiles'', áreas que
contienen una proporción determinada de un conjunto de puntos
cercana a un punto dado. Esta función está implementada 
por \cd{Nearest}. 
También se define una forma de acceso interactiva a esta función
mediante \cd{NearDlg} que toma los datos de un cuadro de diálogo.

%====================================================================

\section{Programa}

 Estructura general del módulo: archivo \cd{env.def} que
declara el \textit{interface}, y archivo \cd{env.mb} con
la implementación.


@o env.def
@{@%
@<Declaraciones P\'ublicas@>
@}

@o env.mb
@{@%
@<Interfaces Usados@>
@<Definiciones@>
@<Declaraciones@>
@<Implementaci\'on@>
@}

 Como en la mayoría de los programas de MapBasic, utilizaremos
definiciones básicas de \cd{mapbasic.def}.
@d Interfaces Usados
@{@%
Include "mapbasic.def"
@}

 También incluimos el interface propio, para declarar las funciones
públicas a implementar.
@d Interfaces Usados
@{@%
Include "env.def"
@}



% ====================================================================

\subsection{Trigonometría}

 En primer lugar se precisan algunas funciones trigonométricas.
Tendremos que usar frecuentemente el valor de la constante $\pi$,
para ello definimos \cd{Pi}. También definimos el identificador
\cd{HalfPi} con el valor $\pi \over 2$ por comodidad y eficiencia:
@d Definiciones
@{@%
Define Pi 3.14159265358979323846
Define HalfPi 1.57079632679489661923
@| Pi HalfPi @}

% --------------------------------------------------------------------

 La función \cd{ATan2} es equivalente a las funciones \cd{atan2}
de Fortran o C; calcula el ángulo (argumento) de un vector
(equivalentemente, de un número complejo).
El valor devuelto es un ángulo en radianes con valor en
$[-\pi,\pi]$.

En primer lugar, la declaración, en el lugar apropiado.
@d Declaraciones
@{@%
Declare Function ATan2(ByVal y As Float, ByVal x As Float) As Float
@}

 Y ahora la definición; cuando el vector es $(0,0)$, el valor
de \cd{Atan2} es indeterminado. Está función no genera ningún tipo de
error; en tal caso de vuelve el valor de la variable
(global visible sólo dentro de este módulo) \cd{ATan2\_Indet}, que
por defecto es 0.

% FIG2: ángulo de un vector, etiquetar theta, x, y
\begin{figure}[!htbp]
\begin{center}
\includegraphics{env.2}
\caption{ángulo de un vector}
\end{center}
\end{figure}

@d Implementaci\'on
@{@%
Dim ATan2_Indet As Float
Function ATan2(ByVal y As Float, ByVal x As Float) As Float
  Dim xsgn, ysgn As Integer
  xsgn = Sgn(x)
  ysgn = Sgn(y)
  Dim a As Float
  a = ATan2_Indet
  If (xsgn<>0 Or ysgn<>0) Then
    x = Abs(x)
    y = Abs(y)
    If (y>x) Then
      a = HalfPi - ATn(x/y)
    Else
      a = ATn(y/x)
    End If
    If (xsgn<0) Then a = Pi-a End If
    If (ysgn<0) Then a = -a End If
  End If
  ATan2 = a
End Function
@| ATan2_Indet ATan2 @}

% --------------------------------------------------------------------

 La función \cd{AngleFromTo} calcula el valor del ángulo que va de
un vector \cd{(x1,y1)} a otro \cd{(x2,y2)} % $\(x_1,y_1\)$ $\(x_2,y_2\)$

% FIG3:  ángulo de un vector a otro
\begin{figure}[!htbp]
\begin{center}
\includegraphics{env.3}
\caption{ángulo entre dos vectores}
\end{center}
\end{figure}

@d Declaraciones
@{@%
Declare Function AngleFromTo(
          ByVal x1 As Float, ByVal y1 As Float,
          ByVal x2 As Float, ByVal y2 As Float
                            ) As Float
@}
@d Implementaci\'on
@{@%
Function AngleFromTo(
       ByVal x1 As Float, ByVal y1 As Float,
       ByVal x2 As Float, ByVal y2 As Float) As Float
  AngleFromTo = ATan2(y2*x1-x2*y1,x2*x1+y2*y1)
End Function
@| AngleFromTo @}


% --------------------------------------------------------------------

% ====================================================================

\subsection{Algoritmo}

 El algoritmo usado para calcular la envolvente convexa de
un conjunto de puntos $(x_n,y_n)$ es el siguiente:

\begin{itemize}

\item
En primer lugar se calcula un punto que será interior a la envolvente;
el centroide o media $(x_c,y_c)$ de los puntos.

\item
A continuación se toma el punto
$({x_i}_0,{y_i}_0)$ más alejado del centroide; este punto pertenece
a la envolvente, y partiremos de él para construirla.

\item
Una vez tenemos un punto $({x_i}_j,{y_i}_j)$ de la envolvente,
el siguiente punto $({x_i}_{j+1},{y_i}_{j+1})$ será aquél que maximice
o minimice el valor de $\alpha$. (maximizando se recorre la envolvente
en el sentido de las agujas del reloj ---el interior queda a la derecha
del perímetro; minimizando se invierte el sentido)

% FIG4: ángulo alpha
% z1 = (3,8)  z2 = (4,3)
% vector (0,0)--z1, etiquetado $p_i_j$ a la derecha de z1
% vector (0,0)--z2, etiquetado $p_i_{j+1}$ a la derecha de z2
% vector z1--z2, etiquetado $p_i_{j+1}-p_i_j$ a la derecha del punto medio
% línea punteada extensión de (0,0)--z1 i.e. de z1 a 1.33z1
% arco orientado de circunferencia de la extensión hasta (z1--z2)
% en sentido contrario agujas reloj, etiquetado a la izda $\alpha$
\begin{figure}[!htbp]
\begin{center}
\includegraphics{env.4}
\caption{parámetro $\alpha$}
\label{fig:alpha}
\end{center}
\end{figure}


\item
El algoritmo finaliza una vez que se vuelve a encontrar el punto inicial.

\end{itemize}


El conjunto de puntos es pasado a la función como una tabla.


El parámetro $\alpha$ a minimizar depende del centroide,
del último punto de la envolvente, $({x_i}_j,{y_i}_j)$, y
de cada punto a considerar. Para evitar pasar un gran
número de argumentos y simplificar las expresiones,
definiremos una función auxiliar que permita fijar el centroide
y el primer punto en variables globales (visibles sólo en el módulo)
@d Declaraciones
@{@%
Dim Alpha_xc, Alpha_yc As Float ' centro de la envolvente
Dim Alpha_x1, Alpha_y1 As Float ' primer vértice
Declare Function Alpha(ByVal x As Float, ByVal y As Float) As Float
@| Alpha_xc Alpha_yc Alpha_x1 Alpha_y1@}

@d Declaraciones P\'ublicas
@{@%
Declare Sub ConvexHull(ByVal point_tab As String, env As Object)
@}

@d Implementaci\'on
@{@%
Sub ConvexHull(ByVal point_tab As String, env As Object)
  Dim xc,yc As Float
  @<Calcular el centroide de \cd{point\_tab} en \cd{xc,yc}@>

  Create Region Into Variable env 1 0
    Center (xc,yc) ' por desgracia esto es ignorado

  Dim x0, y0, xn, yn As Float
  @<Calcular el primer v\'ertice de la envolvente en \cd{x0,y0}@>

  ' Fijamos el centroide para Alpha
  Alpha_xc = xc
  Alpha_yc = yc

  Dim n As SmallInt ' contador de vértices

  ' vértice inicial
  n = 1
  Print Str$(n)+" vértice"
  xn = x0
  yn = y0
  Alter Object env Node Add (xn,yn)

  Do
    ' Fijamos el punto inicial para Alpha
    Alpha_x1 = xn
    Alpha_y1 = yn

    ' vértice siguiente
    n = n+1
    Print Str$(n)+" vértices"
    @<Calcular el punto de \cd{Point\_Tab} que minimice \cd{Alpha} en \cd{xn,yn}@>
    Alter Object env Node Add (xn,yn)

  Loop Until(xn=x0 And yn=y0)

  ' Forzamos a que se recalcule el centroide, ya que si no queda
  ' como uno de los vértices introducidos
  env = ConvertToRegion(ConvertToPLine(env))
  ' también podríamos haber hecho:  env = OverlayNodes(env,env)

End Sub
@| ConvexHull @}

@d Calcular el centroide de \cd{point\_tab} en \cd{xc,yc}
@{@%
  Dim n_points As Integer
  n_points = TableInfo(point_tab,TAB_INFO_NROWS)
  Select Sum(CentroidX(obj)) "X", Sum(CentroidY(obj)) "Y"
    From point_tab Into Tmp_Env
  Fetch First From point_tab
  xc = Tmp_Env.X / n_points
  yc = Tmp_Env.Y / n_points
  Close Table Tmp_Env
@}

@d Calcular el primer v\'ertice de la envolvente en \cd{x0,y0}
@{@%
  Select
    Distance(xc,yc,CentroidX(obj),CentroidY(obj),"m"),
    CentroidX(obj) "X", CentroidY(obj) "Y"
    From Point_Tab
    Order By col1 Desc
    Into Tmp_Env
  Fetch First From Tmp_Env
  x0 = Tmp_Env.X
  y0 = Tmp_Env.Y
  Close Table Tmp_Env
@}

  Al minimizar $\alpha$ hay que tener en cuenta que entre los puntos
\cd{Point\_Tab} se encontrará el punto inicial \cd{Alpha\_x1,Alpha\_y1}.
Si la desigualdad de la función \cd{AngleAlpha} no hubiera sido fijada
correctamente de acuerdo con la operación a realizar,
habría que establecer
el valor \cd{Atan2\_Indet} como \cd{Pi} en caso de minimizar o como \cd{0}
si se va a maximizar $\alpha$.

@d Calcular el punto de \cd{Point\_Tab} que minimice \cd{Alpha} en \cd{xn,yn}
@{@%
   'ATan2_Indet = Pi ' valor para minimizar Alpha

   Select
     Alpha(CentroidX(obj),CentroidY(obj)),
     CentroidX(obj)"X",CentroidY(obj) "Y"
     From Point_Tab
     Order By Col1 Asc
     Into Tmp_Env
   Fetch First From Tmp_Env
   xn = Tmp_Env.X
   yn = Tmp_Env.Y
   Close Table Tmp_Env
@}


% ====================================================================

  Función que calcula parámetro a maximizar, fijados
  el punto central \cd{Alpha\_xc,Alpha\_yc}
  y el punto inicial \cd{Alpha\_x1,Alpha\_y1}.

@d Implementaci\'on
@{@%
Function Alpha(ByVal x As Float, ByVal y As Float) As Float
  Alpha = AngleAlpha(Alpha_x1-Alpha_xc,Alpha_y1-Alpha_yc,
                     x-Alpha_xc,y-Alpha_yc)
End Function
@| Alpha @}

% --------------------------------------------------------------------

% ref{fig:alpha} no funciona, \pageref{fig:alpha} sí. (?)
 Ángulo $\alpha$ a minimizar. (ver figura \ref{fig:alpha}, página \pageref{fig:alpha})
@d Declaraciones
@{@%
Declare Function AngleAlpha(
              ByVal x1 As Float, ByVal y1 As Float,
              ByVal x2 As Float, ByVal y2 As Float) As Float
@}

@d Implementaci\'on
@{@%
Function AngleAlpha(ByVal x1 As Float, ByVal y1 As Float,
                    ByVal x2 As Float, ByVal y2 As Float) As Float
  x2 = x2-x1
  y2 = y2-y1
  Dim a As Float
  a = AngleFromTo(x1,y1,x2,y2)
  If a<=0 Then ' cambiar por a<0 si se va a maximizar Alpha
    a = a+2*Pi
  End If
  AngleAlpha = a
End Function
@| AngleAlpha @}

% ====================================================================

 La función \cd{AuxColVal} resulta muy útil para referirse
a columnas de una tabla cuando el nombre de la tabla o de la
columna se halla en una variable de texto.

Ésta función coincide con la función \cd{ColVal} del módulo
\cd{colval} del programa RTP. Para independizar esté módulo de
aquél programa la implemento también aquí, cambiando el
nombre; cuando 
esté preparada mi biblioteca de utilidades MapBasic, será
incluida allí.

@d Declaraciones
@{@%
Declare Function AuxColVal(ByVal tab As String,
                           ByVal col As String) As Alias
@}

 La implementación original, \cd{AuxColVal=tab+"."+col} da
problemas en algunas expresiones.

@d Implementaci\'on
@{@%
Function AuxColVal(ByVal tab As String,
                   ByVal col As String) As Alias
  'AuxColVal = tab+"."+col ' esto da problemas en ciertos casos
  Dim a As Alias
  a = tab+"."+col
  AuxColVal = a
End Function
@| AuxColVal @}

% ====================================================================

\subsection{Puntos cercanos}

 La función \cd{Nearest} calcula el área que contiene a la
fracción más cercana a \cd{(xc,yc)} de los puntos de la tabla
\cd{tabname}.

Antes de usar esta función se debe establecer el sistema de
coordenadas apropiado para \cd{(xc,yc)} mediante una instrucción
del tipo \cd{Set CoordSys...}.

@d Declaraciones P\'ublicas
@{@%
Declare Sub Nearest
  (ByVal tabname As String,              ' tabla de puntos
   ByVal f As Float,                     ' fracción de puntos
   ByVal xc As Float, ByVal yc As Float, ' centro
   env As Object)                        ' resultado
@}
@d Implementaci\'on
@{@%
Sub Nearest(ByVal tabname As String,
            ByVal per As Float,
            ByVal xc As Float, ByVal yc As Float,
            env As Object)
   If per<1 Then
     Select Distance(xc,yc,CentroidX(obj),CentroidY(obj),"m")
        From tabname
        Order By Col1
        Into Tmp_Nearest0
     Dim nrows, threshold As Integer
     nrows = TableInfo(Tmp_Nearest0,TAB_INFO_NROWS)
     threshold = nrows*per
     Select * From Tmp_Nearest0
       Where RowID<=threshold Into Tmp_Nearest
   Else
     Select * From tabname Into Tmp_Nearest
   End If
   Call ConvexHull("Tmp_Nearest",env)
   Close Table "Tmp_Nearest"
   If per<1 Then
     Close Table "Tmp_Nearest0"
   End If
End Sub
@| Nearest @}


% --------------------------------------------------------------------
\subsection{Interfaz de Usuario}

 La función \cd{NearDlg} presenta un {em interface} interactivo
para \cd{Nearest}. La envolvente creada la coloca en la capa
descriptiva de la ventana de mapa seleccionada.

El uso de esta función sería así:
\begin{itemize}
\item
Si no se dispone de una tabla que contenga únicamente el punto central,
seleccionar este punto (crearlo si es necesario en una capa descriptiva)
\item
Ejecutar la utilidad que ejecute \c{NearDlg}
\item
En centro elegir la tabla que contiene el centro (e.g. la selección)
\item
En puntos elegir la tabla de los puntos a envolver
(si la tabla contiene elementos no puntuales se envolverán sus centroides)
\item
Introducir el porcentaje de puntos a envolver
\item
Elegir una ventana de mapa, en cuya capa descriptiva se colocará
   la envolvente generada.
\end{itemize}


@d Declaraciones P\'ublicas
@{@%
Declare Sub NearDlg
@}

@D Implementaci\'on
@{@%
Sub NearDlg

  Dim i As SmallInt

  ' Vector de tablas mapificables
  Dim maptabs(1) As String
  Dim n_maptabs As SmallInt

  @<Llena vector de tablas mapificables, \cd{n\_maptabs},\cd{maptabs} @>

  ' Vector de ventanas de mapa
  Dim mapwins(1) As Integer
  Dim maptitles As String
  Dim n_mapwins As SmallInt

  @<Llena el vector de ventanas de mapa, \cd{n\_mapwins,mapwins,maptitles}@>

  ' Valores seleccionables en el diálogo
  Dim center_i, points_i, map_i As SmallInt
  Dim percent As Float

  @<Establece valores iniciales de \cd{center\_i}, \cd{points\_i}, \cd{maps\_i}@>

  ' establece fracción inicial (50%)  
  percent = 50  

  @<Cuadro de di\'alogo@>

  If CommandInfo(CMD_INFO_DLG_OK) Then
    ' Se debe haber pulsado OK...
    If center_i>0 And points_i>0 And map_i>0 Then
      ' ...con selecciones válidas...

      ' nombre de la tabla que contiene el punto central
      Dim tabname As String
      tabname = maptabs(center_i)
      If tabname="Selección" Then
        tabname = "Selection"
      End If          

      If TableInfo(tabname,TAB_INFO_NROWS)=1 Then
        Fetch First From tabname
        If Not AuxColVal(tabname,"Obj") Then
          Note "El centro debe ser un único punto"
        Else
          @<Crear \'area en el mapa \cd{map\_i} con \cd{tabname}, \cd{percent}@>
        End If
      Else
        Note "El centro debe ser un único punto"
      End If
    End If
  End If

End Sub
@| NearDlg @}


% --------------------------------------------------------------------

@d Llena vector de tablas mapificables, \cd{n\_maptabs},\cd{maptabs}
@{@%
  Dim n_tabs As SmallInt

  n_tabs = NumTables()
  ReDim maptabs(n_tabs+1)
  n_maptabs = 0

  ' Si hay una selección activa y mapificable, se establece
  ' como primera tabla, con el nombre "Selección"
  If SelectionInfo(SEL_INFO_NROWS)>0 Then
     If TableInfo(SelectionInfo(SEL_INFO_TABLENAME),TAB_INFO_MAPPABLE) Then
        n_maptabs = n_maptabs+1
        maptabs(n_maptabs) = "Selección"
     End If
  End If

  For i = 1 to n_tabs
    If TableInfo(i,TAB_INFO_MAPPABLE) Then
      n_maptabs = n_maptabs+1
      maptabs(n_maptabs) = tableinfo(i,TAB_INFO_NAME)  
    End If
  Next

  ReDim maptabs(n_maptabs)
@}

% --------------------------------------------------------------------

@d Llena el vector de ventanas de mapa, \cd{n\_mapwins,mapwins,maptitles}
@{@%
  Dim n_wins As SmallInt
  Dim win_id As Integer

  n_wins = NumWindows()
  ReDim mapwins(n_wins)
  n_mapwins = 0

  For i = 1 to n_wins
    win_id = WindowID(i)
    If WindowInfo(win_id,WIN_INFO_TYPE)=WIN_MAPPER Then
      n_mapwins = n_mapwins+1
      mapwins(n_mapwins) = win_id
      maptitles = maptitles + WindowInfo(win_id,WIN_INFO_NAME) + ";"
    End If
  Next
  ReDim mapwins(n_mapwins)
@}

% --------------------------------------------------------------------

@d Establece valores iniciales de \cd{center\_i}, \cd{points\_i}, \cd{maps\_i}
@{@%
  ' centro: selección o primera tabla mapificable
  center_i = Minimum(1,n_maptabs)

  ' puntos: primera tabla mapificable, no la selección
  points_i = center_i
  If maptabs(points_i)="Selección" Then
    points_i = Minimum(points_i+1,n_maptabs)
  End If

  ' mapa: primera ventana de mapa
  map_i = Minimum(1,n_mapwins)
@}


% --------------------------------------------------------------------

@d Crear \'area en el mapa \cd{map\_i} con \cd{tabname}, \cd{percent}
@{@%
 ' punto central; ---objeto de la tabla tabname
 Dim cnt As Object
 cnt = AuxColVal(tabname,"Obj")

 ' coordenadas del punto central
 Dim cntx, cnty As Float
 cntx = CentroidX(cnt)
 cnty = CentroidY(cnt)

 ' fracción de los puntos; ---percent%
 Dim fract As Float
 fract = percent*0.01

 ' área buscada
 Dim env As Object
 Call Nearest(maptabs(points_i), fract, cntx, cnty, env)

 ' Insertamos el área en la capa descriptiva del mapa; ---ventana map_i
 Dim restab As String
 restab = LayerInfo(mapwins(map_i),0,LAYER_INFO_NAME)
 Insert Into restab (Obj) Values (env)
@}

% --------------------------------------------------------------------
El cuadro de diálogo emplea las variables \cd{maptabs}, \cd{maptitles}, 
valores \cd{center\_i}, \cd{points\_i}, \cd{percent} y \cd{map\_i}.

@D Cuadro de di\'alogo
@{@%
  Dialog
    Title "Área Envolvente"
    Width 178 Height 94

    Control StaticText
      Title "Centro"
      Position 5,11
    Control PopUpMenu
       Value center_i
       Into center_i
       Title From Variable maptabs
       Position 60,10 Width 115 'Height 12

    Control StaticText
      Title "Puntos"
      Position 5,26
    Control PopUpMenu
       Value points_i
       Into points_i
       Title From Variable maptabs
       Position 60,25 Width 115 'Height 12

    Control StaticText
      Title "Fracción de los puntos (%)"
      Position 5,41          
    Control EditText
       Value percent
       Into percent
       Position 95,40 Width 80 'Height 12

    Control StaticText
      Title "Resultado en"
      Position 5,56          
    Control PopUpMenu
       Value map_i
       Into map_i
       Title From Variable maptitles
       Position 60,55 Width 115 'Height 12

    Control OKButton Title "&Aceptar"
       Position 85,74
    Control CancelButton Title "&Cancelar"
       Position 130,74
@}

\section{Índices}

\subsection{Archivos}
@f

\subsection{Fragmentos}
@m

\subsection{Identificadores}
@u





\end{document}