Adaptación e implementación en Matlab del método de Ghosh y Mount para el cálculo del grafo de visibilidad

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo de Fin de Máster Máster en Ingeniería Aeronáutica Adaptación e implementación en Matlab del método de Ghosh y Mount para el cálculo del grafo de visibilidad Autor: Guillermo Vallejo Soto Tutor: Antonio Franco Espín Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2024 Trabajo de Fin de Máster Máster en Ingeniería Aeronáutica Adaptación e implementación en Matlab del método de Ghosh y Mount para el cálculo del grafo de visibilidad Autor: Guillermo Vallejo Soto Tutor: Antonio Franco Espín Profesor Titular Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2024 Trabajo de Fin de Máster: Adaptación e implementación en Matlab del método de Ghosh y Mount para el cálculo del grafo de visibilidad Autor: Guillermo Vallejo Soto Tutor: Antonio Franco Espín El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos L a realización de este trabajo ha supuesto todo un reto a nivel personal. Parafraseando al ilustre filósofo Ortega y Gasset: Yo soy yo y mi circunstancia, y en esta ocasión salvarla ha implicado la entrada al mercado laboral, emigrar al extranjero y aprender un nuevo idioma desde sus bases. Todo esto no hubiera sido posible sin el apoyo de mi pareja, mi familia, mis amigos y la confianza mostrada por mi tutor, Antonio Franco. Muchas gracias, sin vosotros, esto no hubiera sido posible. Dankbarkeit ist das Gedächtnis des Herzens. Guillermo Vallejo Soto Máster en Ingeniería Aeronáutica Hannover, 2024 I Resumen E l grafo de visibilidad de un dominio es un grafo no direccionado cuyos vertices son los vértices del dominio y cuyas aristas son las líneas de visión directa entre los vértices no interrumpidas por los obstáculos del dominio, es por ello que su cálculo es crítico a la hora de resolver problemas de trazado óptimo de trayectorias. A lo largo del siguiente documento se hará una presentación acerca de las motivaciones que pudieran llevar a querer resolver este problema y una concisa revisión acerca del estado del arte. Se desarrollará una implementación en Matlab™realizada en base al algoritmo propuesto por Ghosh & Mount y se establecerán una serie de metodologías para ponerlo a prueba y comprobar su eficacia y efectividad. III XÍndice Almacenamiento y acceso a las familias 26 2.3 Bucle en ii 30 2.3.1 Nodos ii =1yii =230 2.3.2 Ejecución normal del bucle: Triangulación 30 2.3.3 Ejecución normal del bucle: SPLIT 40 Puntualización sobre la definición de visibilidad 43 Cálculo de u44 Calculo de u′50 Cálculo de qyr52 Cálculo de t55 Vector tk57 Últimas comprobaciones 60 2.3.4 Procesado de las familias 62 Procedimiento: Longitud de vectorrecorre es inferior a 2 63 Adición de los nodos asociados a las familias al diagrama de visibilidad extendido 63 3 Estudio de resultados 67 3.1 Comprobación del caso de estudio: Método de Lee en tormenta realista 67 3.2 Ejecución en casos propuestos: Recintos convexos 69 3.3 Resultados de la ejecución 71 3.3.1 Complejidad algorítimica real 73 3.3.2 Estimación matemática de E75 3.3.3 Desglose de tiempos utilizando profiler 76 4 Conclusiones y futuras líneas de trabajo 79 Apéndice A Códigos de Matlab 81 A.1 Relaciones entre los códigos 81 A.2 Programas parent 81 A.2.1 Ejecutaprograma 81 A.2.2 Estudio_parametrico 82 A.2.3 Compatibiliza 83 A.3 Programas Children 85 A.3.1 calcula_recinto 85 A.3.2 Main_v13_C 87 calcula_fronteras 87 test_de_paternidad_MAIN 89 cortapuntos 92 Visible_DEF 94 inter_visible_1 97 calculaCW_CCW:sucesores 100 checkIfInsideTriangle 101 A.3.3 add_familia2 102 futurasgeneraciones 112 generalistanodos 115 ordenalistasec 116 Índice de Figuras 117 Índice de Tablas 119 Índice XI Índice de Códigos 121 Bibliografía 123 Notación VMatriz de coordenadas WPs Vector de puntos de paso NNúmero de puntos de paso diferentes PDominio P ekConjunto de aristas perteneciente al dominio SRecta horizontal que pasa por un punto TConjunto de aristas intersecados por S x,yCoordenadas en espacio cartesiano xv Base izquierda de la familia de embudos (x,y) dividida por ii vy Base derecha de la familia de embudos (x,y) dividida por ii EVG Diagrama de Visibilidad extendido i j Segmento de recta desde el punto ial punto j  i j Vector director de iaj vVectorv ∥v∥Norma del vector v v·wProducto escalar de los vectores vyw v×wProducto vectorial de los vectores vyw qzVector unitario en la dirección perpendicular al plano z i,j,kÍndices comunes NaN Not a Number ∈Que se encuentra en |A|Determinante de la matriz cuadrada A det(A)Determinante de la matriz (cuadrada) A log Logaritmo de base natural eNúmero e sen Función seno tg Función tangente ⩽Menor o igual ⩾Mayor o igual ⋍Aproximadamente igual TM Trade Mark XIII 1 Motivación, estado del arte y objetivos del trabajo E xisten varias y diversas razones por las cuales resulta que en la práctica, la computación del diagrama de visibilidad es una necesidad ineludible. Dado que el contexto del redactor es aeroespacial, se partirá de esta base. 1.1 Importancia del trazado de ruta en el mundo aeronáutico Las operaciones aeronáuticas están fuertemente reguladas con el propósito de mantener los niveles de seguridad operacional y física a la vanguardia del mundo ingenieril, estableciendo estándares de calidad que superan las de otros sectores estratégicos. La preparación y ejecución de las rutas de los vuelos comerciales es una de estas actividades críticas, dado que un buen trazado proporciona beneficios a varios niveles: • Ecológicos: Con los objetivos de la Agenda 2030 en mente, enfatizando en los objetivos número 13 Acción por el clima y el número 9 Industria, Innovación y Crecimiento económico, la optimización de trayectorias de ruta es una alternativa factible para reducir el consumo de combustible combinada con otras estrategias. Reduciendo las emisiones de gases de efecto invernadero y haciendo la industria aeronáutica más sostenible. • Sociales: Existen zonas o bien cerradas al tráfico aéreo por razones de seguridad nacional, o bien por zonas de limitación de ruido dependientes del horario. Una ruta óptima deberá añadir estos factores a la ecuación para minimizar el impacto de las actividades aeronáuticas. • Económicos: Las rutas óptimas no sólo reducen el consumo de combustible, si no que, además, tienen un doble efecto de ahorro para las compañías que ofrecen estos servicios al gran público. El coste operacional de las aerolíneas es fuertemente dependiente de los gastos derivados del combustible empleado, [1]. Es por ello que, pequeñas reducciones en estos gastos pueden tener un gran efecto en la estrategia financiera. Por otro lado, reducir las emisiones de CO2 dispone en la Unión Europea a pagar cánones menores, disminuyendo por partida doble los costos operativos. Una buena planificación de ruta aplicada a todos los agentes de la industria puede evitar retrasos y los costes de los mismos para las aerolíneas. • Seguridad: Un trazado de rutas óptimo tiene impacto en la seguridad operacional, evitando accidentes que pudieran producirse debido a eventos climatológicos de alto impacto en las operaciones aeronáuticas, como las borrascas o nubes que pudieran reducir la visibilidad para los pilotos. Además conforma una herramienta de apoyo excelente para control de tráfico 1 2Capítulo 1. Motivación, estado del arte y objetivos del trabajo aéreo, ATM, pudiendo aumentar la capacidad de los sectores sin perjudicar los estándares de seguridad mencionados anteriormente. Esto es sólo el enfoque referente a la aviación comercial, pero a pequeña escala, por ejemplo, la industria de los vehículos aéreos no tripulados, el diseño y optimización de rutas acorde a la misión propuesta, permite una operación controlada y libre de riesgos de colisión. Estas aplicaciones serían de alto interés en la entrega urgente y autónoma de mercancía, como en el caso del hospital de la Paz en Madrid, que participa de forma activa en el desarrollo de un proyecto que permita estandarizar el uso de UAV con objeto de optimizar la entrega de material sanitario de alto interés, [2]. 1.2 Otras aplicaciones de interés Estas ventajas no se limitan exclusivamente al entorno aeroespacial, por lo que se hará a lo largo de esta sección un pequeño desglose acerca del impacto del cálculo de rutas, y por ende del grafo de visibilidad, en otras industrias. 1.2.1 Robótica A día de hoy, los robots autónomos no son ciencia ficción, y aunque distan de la visión popular que se tenía de los mismos, provocada en gran parte por la industria del entretenimiento, son muchas las empresas que los usan para optimizar sus procesos. Entre estos procesos se encuentran, por ejemplo, la gestión de almacenes y paquetería, donde el trazado de rutas óptimas en espacios reducidos y con presencia de otros agentes dinámicos, es imprescindible. Amazon y Tesla se presentan como referentes en la industria, tanto en sacarle ventaja al uso de estas tecnologías como en accidentes laborales en trabajadores de la empresa a causa de estándares de seguridad deficientes, [3, 4]. 1.2.2 Geoinformática y sistemas de información geográfica Un trazado de rutas óptimo puede llegar a tener impacto a la hora de representar rutas y mapear zonas accesibles o intransitables. Esto puede ser altamente beneficioso a la hora de evaluar zonas seguras en zonas propensas a sufrir deslizamientos, avalanchas o erupciones volcánicas, y trazado rutas óptimas en situaciones de emergencia y rescate. 1.2.3 Ingeniería en redes y telecomunicaciones Los grafos de visibilidad pueden ser empleados para determinar las posiciones óptimas para ubicar antenas de radio, telefonía y televisión, maximizando la cobertura mediante la elusión de obstáculos naturales y artificiales. Garantizando la cobertura de estos servicios en áreas pobladas, o disminuyendo los costes operacionales, facilitando que estos mismos lleguen a zonas más apartadas, como por ejemplo, la España vaciada. 1.2.4 Arquitectura y urbanismo De forma similar a como se hacía en la aplicación propuesta de gestión de almacenes, resulta de interés usar un enfoque similar para poder trazar rutas de accesibilidad y evacuación eficientes. Pudiéndose llegar a optimizar estas rutas de cara a la accesibilidad. Planteando, entre otras alternativas, rutas posibles para personas con movilidad reducida y mejorando la infraestructura urbana. Por estos y otros tantos motivos, el trazado y optimización de rutas, en especial el cálculo de grafos de visibilidad, es un problema cuyo perfeccionamiento es de alto interés estratégico de cara a los retos del mañana. 1.3 Estado del arte 3 1.3 Estado del arte En esta sección se repasarán diferentes métodos de cálculo de grafos de visibilidad empleados en la actualidad, junto a la complejidad asociada a los mismos. 1.3.1 Método Ingenuo A la hora de resolver un problema, suele ser interesante comprobar cual es la resolución posible que podría alcanzarse mediante el empleo de la fuerza bruta, para poder tener un caso de referencia y comprobar como de bueno es el algoritmo que se propone frente a él. Esta serie de métodos son comúnmente denominados como Ingenuos oNaive. Para el caso de la determinación del grafo de visibilidad el método Ingenuo propone recorrer cada punto de paso, denominados como Waypoints,WPs, comprobando si la visibilidad con el resto de puntos de paso,WPs se ve interrumpida por alguna arista obstáculo eP del dominio P, de forma adicional, es necesario comprobar si esta linea de visibilidad no cruza dentro de un obstáculo. Este problema puede simplificarse ligeramente si se consideran las propiedades simétricas del problema, si i j es visible, ji también lo será, y viceversa. Esto reduce el número de comprobaciones necesarias y acelera el proceso de cálculo. Se adjunta el pseudocódigo del programa que ejecuta lo anteriormente expresado, este método fue implementado satisfactoriamente en Matlab en el Trabajo de Fin de Grado de Narciso Valverde [5]. Pseudocódigo 2.1 Cálculo del grafo de Visibilidad mediante algoritmo Ingenuo Grafo de visibilidad = función_algoritmo_ingenuo(Matriz de coordenadas V) 1: para cada punto de paso i∈WPs,hacer 2: para cada punto de paso j∈(i+1,nWPs), hacer 3: para cada arista obstáculo k∈eP,hacer 4: si no hay intersección entre la arista i j yeky además i j no pasa dentro de un obstáculo contenido en P,entonces 5: Añadir al Grafo de visibilidad i j 6: devolver Grafo de visibilidad La comprobación de si la potencial línea de visibilidad pasa por dentro del polígono al que pertenecen puede resolverse si se tiene información de antemano acerca de los puntos del polígono ordenados. Esto forma parte de la premisa, por lo que se puede utilizar estos datos proporcionados junto con los productos vectoriales para minimizar el número de operaciones. El fundamento se basa en el siguiente precepto: todo punto en un polígono plano convencional está conectado a otros dos puntos, uno que lo precede y otro que lo sucede. Este arco que forman puede ser cóncavo o convexo, pero cualquier línea de visión que parta de ese punto saldrá hacia el exterior del polígono o hacia el interior del polígono y eso puede comprobarse usándo productos vectoriales. Se muestra un ejemplo visual en la figura 1.1. Determinar si el arco formado es cóncavo o convexo para un punto y su predecesor/sucesor puede hacerse realizando el producto vectorial entre ambos, teniendo en mente que están proporcionados los puntos del polígono respecto a su centroide y se disponen en el sentido contrario a las agujas del reloj, manteniendo la nomenclatura de la figura 1.1: a)v1×v2·qz>0,(1.1) b)v1×v2·qz<0,(1.2) 4Capítulo 1. Motivación, estado del arte y objetivos del trabajo a) Recinto cóncavo b) Recinto convexo i+1 i i-1 v1 v2 ij i+1 i i-1 v1 v2 ij Figura 1.1 Uso de productos vectoriales para determinar si una línea cruza por dentro de un polígono. donde el caso a) representa concavidad localizada y b) convexidad localizada. Una vez conocedores de esto, se puede distinguir a partir del uso de inecuaciones si pasa por dentro o no del polígono dependiendo del caso: 1. Caso concavidad localizada,  i j fuera del polígono:  i j ×v1·qz<0, i j ×v2·qz>0.(1.3) 2. Caso convexidad localizada,  i j dentro del polígono:  i j ×v1·qz>0, i j ×v2·qz<0.(1.4) Esta propiedad ha sido implementada en la algoritmia del método de Ghosh & Mount, [6], dado que permite estandarizar cálculos en recintos cóncavos con bajo coste computacional. Es uno de los códigos empleados en la sección 2.2.3. La complejidad algorítmica de este método propuesto en notación Big(O)es: O(N)≈N3,(1.5) debido a que es necesario anidar bucles tres veces recorriendo el conjunto de puntos completo, siendo Nel número de diferentes puntos de paso. Partiendo de esta base, cualquier algoritmo que mejore estas condiciones manteniendo la validez de los resultados representará una ventaja desde un punto de vista computacional. 1.3.2 Método de Lee El método de Lee es otra metodología que permite obtener el grafo de visibilidad empleando una estrategia que solo comprueba las visibilidades con la arista más cercana al punto que se está añadiendo al grafo. Toda la explicación desarrollada a continuación está basada en el trabajo de investigación realizado por Narciso Valverde [5], y sus resultados han sido utilizados como referencia a la hora de trazar comparaciones con el método de Ghosh & Mount. Para realizar el cálculo del grafo se recorren todos los puntos de paso y para cada uno de ellos se procesa el semiplano superior a la horizontal que pasa por el punto a añadir al grafo. Es posible 1.3 Estado del arte 5 realizar esto debido a que el problema es simétrico, es decir, los puntos que están en contacto visual de norte a sur, son los mismos que están en contacto visual de sur a norte. Por lo que procesando los semiplanos superiores para todos los puntos de paso queda cubierto completamente el dominio. S T ={e1,e3} e1 e2 e3 e4 e5 e6 e7 e8 e9 e10 e11 e12 e13 A B C D E F G H I J K L M a) Arranque de la iteración S e1 e2 e3 e4 e5 e6 e7 e8 e9 e10 e11 e12 e13 A B C D E F G H I J K L M b) Rotación CCW de S hasta vértice C 00 T ={e1,e3} i ¿intersección(iC,e1)? No, iC se añade al grafo 1 T ={} S e1 e2 e3 e4 e5 e6 e7 e8 e9 e10 e11 e12 e13 A B C D E F G H I J K L M c) Rotación CCW de S hasta vértice D S e1 e2 e3 e4 e5 e6 e7 e8 e9 e10 e11 e12 e13 A B C D E F G H I J K L M d) Rotación CCW de S hasta vértice E 2 T ={e4,e5} i ¿intersección(iE,e4)? No, iE se añade al grafo 3 T ={e7,e5} 1 T ={} 2 T ={e4,e5} Figura 1.2 Ejecución método de Lee para un punto de paso. Todo este proceso está reflejado en la figura 1.2, y se irán mencionando los diferentes puntos y etapas referenciando a este soporto visual. El siguiente procedimiento se ejecuta para cada punto de paso. Primero se traza una línea horizontal a través del punto de estudio i , se comprueban los cortes de esta recta S con el resto de polígonos del dominio, cada una de las aristas con las que corta se almacenan en orden de corte en el vector T , paso a) de la figura. Una vez inicializado este vector, se visitan en el sentido contrario a las agujas del reloj los puntos de paso ubicados en el semiplano superior. Cada vez que se visita un punto de paso se comprueba si existe un corte entre la línea de visión y el punto que se visita. En caso de que no existan cortes o no haya aristas almacenadas en T el segmento se añade al grafo de visibilidad. Después de esto, se actualiza el vector T , siempre que se visite un punto puede darse cualquiera de los tres siguientes casos: 1. El punto que se visita pertenece a un nuevo polígono, por lo que hay que añadir dos nuevas aristas a T. Esto puede observarse en la figura en el caso c. 12 Capítulo 2. Método de Ghosh & Mount ii +1sucede a ii. 72 %% Cálculo de las conexiones entre nodos de un polígono 73 % Esta matriz es importante ya que es la que regulará el cálculo de u vector, 74 % y es lo que evitará que se creen y consideren embudos que provocarían 75 % falsos positivos en el grafo de visibilidad. 76 77 78 %Consec_aristas registrará solo las aristas consecutivas 79 Consec_aristas=zeros(N,N); 80 for ii=2:length(V(:,1)) 81 if ~isnan(V(ii-1,1)) && ~isnan(V(ii,1)) 82 A=find(WPs(:,1)==V(ii-1,1)); 83 B=find(WPs(:,1)==V(ii,1)); 84 Consec_aristas(A,B)=1; 85 Consec_aristas(B,A)=1; 86 end 87 end 88 EVG.Consec_aristasrec=zeros(N,2); 89 for ii=1:N 90 ind=find(Consec_aristas(ii,:)==1); 91 EVG.Consec_aristasrec(ii,:)=ind; 92 end 93 94 Consec_aristas2=zeros(N,N); 95 for ii=2:length(V(:,1)) 96 if ~isnan(V(ii-1,1)) && ~isnan(V(ii,1)) 97 Consec_aristas2(WPs(:,1)==V(ii-1,1),WPs(:,1)==V(ii,1))=1; 98 end 99 100 end 101 102 EVG.Consec_aristas2=Consec_aristas2; 103 EVG.Consec_aristas=Consec_aristas; Código 2.3 Registro de consecución de los nodos en los polígonos. Nótese que de existir valores repetidos en la coordenada x , pudiera darse que para un punto de paso dado existieran más de dos aristas confluyendo hacia un mismo vértice. 2.2.2 Matrices de polígonos En esta sección se guardará información de rápido acceso acerca de cuantos polígonos hay y que vértices los componen. La información que se extraerá es la siguiente: • poligono2nodo: Es una matriz de dimensiones npol xN , siendo npol el número de polígonos, que pudieran causar interferencias en la visión, es decir, sin considerar el polígono exterior. El polígono que conforma el recinto exterior se considera el polígono 0 . Si se quisiera consultar que vértices tiene el polígono i , se accedería a la fila i de la matriz, y se recortarían los valores distintos de 0, que se usan de relleno para mantener consistencia en el tamaño de la matriz. • nodo2poligono: Es una matriz N x1 que proporciona la información inversa a poligono2nodo. Es decir, si se tuviera un punto de paso a que se quisiera saber a que polígono pertenece sólo sería necesario acceder a la información presente en la fila a. • limitespoligono: Es una matriz npol x2 , que registra cual son los puntos de paso de menor y mayor coordenada x para un polígono dado. Esto es interesante de cara a la identificación de aristas de intersección durante la triangulación. poligono2nodo Primero se definen dos variables bandera: f lag , que se mantiene bajada hasta que se encuentra el primer valor de NaN , que si está bien definida la matriz, será una vez pasadas las coordenadas 2.2 Arranque del programa 13 del polígono 0 ; por otro lado f lag2 sirve para indicarle al programa que cambie de fila una vez identifica que se ha repetido el punto inicial del polígono. Una vez definidas las banderas, se recorren todos los valores pertenecientes a la matriz completa de coordenadas V . Cuando f lag2 se levanta, se guarda la información en poligono2nodo junto con una fila de ceros para cubrir todas las columnas. 105 %% Matriz polígonos 106 %Esta matriz será tal que npol x N, contendrá información 107 %acerca de que aristas contiene cada poligono (Ojo, el recinto exterior no cuenta como polígono.) 108 flag=0; 109 npol=1; 110 ii=1; 111 jj=1; 112 113 while ii<=length(V(:,1)) 114 if flag==1 115 flag2=0; 116 while flag2==0 117 118 %Solo hay un valor de x en Waypoints que coincida 119 P=find(WPs(:,1)==V(ii,1),1); 120 121 if jj~=1 122 %Si no es el primer vértice del poligono, miramos si está repetido 123 if EVG.poligono2nodo(npol,1)==P 124 flag2=1;%Asi se sale del bucle que corre para cada secuencia de vértices que identifica el polígono 125 jj=1; 126 else 127 EVG.poligono2nodo(npol,jj)=P; 128 jj=jj+1; 129 end 130 else 131 EVG.poligono2nodo(npol,1)=P; 132 jj=2; 133 end 134 ii=ii+1; 135 end 136 EVG.poligono2nodo(npol,1:N)=[EVG.poligono2nodo(npol,:),zeros(1,N-length(EVG.poligono2nodo (npol,:)))]; 137 npol=npol+1; 138 end 139 140 if flag==0 && isnan(V(ii,1)) 141 flag=1; 142 end 143 ii=ii+1;%Vale tanto para avanzar en el caso del recinto externo como para salir del NaN en el que dejaría salir del bucle de identificación de vértices. 144 end Código 2.4 Identificación poligono2nodo. nodo2poligono Usando la información de poligono2nodo y recorriendo sus filas se puede registrar en un vector columna la relación de pertenencia de un nodo a un polígono. Es decir realizando la consulta nodo2poligono(i) , se obtendría el número de polígono (sin contar el recinto externo) al que pertenecería. Esto es especialmente útil para identificar a partir de un nodo concreto, cuales son los otros nodos que pertenecen al mismo polígono. Esto se emplea para calcular posibles interferencias de visibilidad con las aristas. Se recuerda al lector que todos los nodos están identificados por orden de xcreciente, tal y como están dispuestos en WPs. 146 EVG.nodo2poligono=zeros(N,1); 147 for ii=1:length(EVG.poligono2nodo(:,1)) 148 for jj=1:length(EVG.poligono2nodo(1,:)) 14 Capítulo 2. Método de Ghosh & Mount 149 if EVG.poligono2nodo(ii,jj)==0 150 break 151 end 152 EVG.nodo2poligono(EVG.poligono2nodo(ii,jj))=ii; 153 end 154 end Código 2.5 Identificación nodo2poligono. limitespoligono Sacando partido de las propiedades de nodo2poligono se recorren todos sus valores, no es necesario realizar este registro para el polígono 0 , en caso contrario, hay que comprobar para la fila asociada al polígono si hay un valor en la primera columna ya asociado, si no es así, se guarda, si ya había un valor en la columna izquierda, se guarda en la columna de la derecha. Esto es así porque al estar ordenados en orden creciente de coordenada x en el vector WPs , el primer valor que aparezca siempre será el más occidental, mientras que el último el más oriental. Para evitar meter varias comprobaciones redundantes, se sobrescriben datos en la columna de la derecha, la de mayor valor de coordenada x , hasta alcanzar el último valor, que provocará que a partir de ahí, no vuelvan a sobreescribirse valores. 156 EVG.limitespoligono=zeros(length(EVG.poligono2nodo(:,1)),2); 157 158 for ii=1:length(EVG.nodo2poligono) 159 160 if EVG.nodo2poligono(ii)~=0 161 if EVG.limitespoligono(EVG.nodo2poligono(ii),1)==0 162 EVG.limitespoligono(EVG.nodo2poligono(ii),1)=ii; 163 else 164 EVG.limitespoligono(EVG.nodo2poligono(ii),2)=ii; 165 end 166 end 167 end 168 % EVG.limitespoligono 169 %Esta variable limitespoligono guarda el primer punto y el último, en orden de polígono, de cada polígono Código 2.6 Identificación limitespoligono. 2.2.3 Caminos prohibidos Para optimizar los tiempos de procesamiento, se excluirán los caminos internos que atraviesan los polígonos, tanto en su totalidad como de manera parcial (en el caso de polígonos no convexos), ya que por definición son invisibles al cruzar áreas consideradas prohibidas. En esta fase del preprocesado, se extraerán los valores de los nodos vinculados a cada polígono y se enviarán a la subrutina calculafronteras, el cual se explicará detalladamente en el Apéndice A. Esbozando el funcionamiento del programa, se traza para cada par de vértices una linea auxiliar y se realizan las siguientes comprobaciones: 1. DoesCrossPolygon: Que comprueba si existen cortes de la línea auxiliar con alguna arista que no sean las que conforman el destino y el final de la línea auxiliar. Si se produce un corte bajo estas condiciones, ese camino ya no es factible. 2. isOutsideArc: Comprueba, a través de un calculo de productos vectoriales, si ese camino, pasa por dentro o no del polígono. Para que un camino sea posible es necesario que no cruce ninguna arista que pertenezca al polígono de estudio y que este camino siempre sea externo y nunca cruce por dentro del mismo. 2.2 Arranque del programa 15 calcula_f ronteras da como resultado un matriz de nodospolpolk x nodospolpolk , los valores iguales a 1 son los posibles caminos que cumplen las dos condiciones citadas anteriormente. Identificando estos puntos (teniendo en cuenta que naristas guarda los identificadores numéricos de los waypoints), es posible entonces establecer una matriz que proporcione la información contraria, es decir, que caminos son prohibidos. 172 %% Caminos prohibidos 173 Caminos_prohibidos=zeros(N,N); 174 175 for jj=1:length(EVG.poligono2nodo(:,1)) 176 %Extraemos poligono 177 npol=jj; 178 naristasextend=EVG.poligono2nodo(npol,:); 179 indcero=find(naristasextend==0,1); 180 naristas=naristasextend(1:indcero-1)’; 181 182 polygonVertices=WPs(naristas,:); 183 possiblePaths=calcula_fronteras(polygonVertices); 184 185 %Ahora se adapta esto al formato de caminos prohibidos 186 %possiblepaths(a,b) es 1 cuando el camino es posible, por lo que el camino no está prohibido 187 188 % possiblepaths tiene el formato de naristas por naristas 189 190 for kk=1:length(naristas) 191 for ll=kk+1:length(naristas) 192 if possiblePaths(kk,ll)==0 %O sea si el camino no es posible, está prohibido 193 a=naristas(kk); 194 b=naristas(ll); 195 Caminos_prohibidos(a,b)=1; 196 Caminos_prohibidos(b,a)=1; 197 end 198 end 199 end 200 end 201 202 % % 1 no permitido, 0 permitido 203 % %Las aristas consecutivas son caminos permitidos. Código 2.7 Identificación de Caminos_prohibidos. 2.2.4 Cálculo de las matrices de ordenación La siguiente sección entronca directamente con uno de los puntos claves del documento original de referencia de Ghosh & Mount, [6], que se basa en la idea de que la complejidad propuesta del algoritmo es alcanzable siempre y cuando los datos y operaciones que se narrarán a continuación sean calculables con complejidad O(1) . Todo esto se fundamenta en el Extended Visibility Graph, en la implementación realizada se ha creado una estructura de datos a la que se hace referencia en todas las funciones y subfunciones con las siglas en inglés EVG. Por ejemplo, nodo2poligono se ha almacenado dentro de la estructura de datos EVG , dado que contiene información de interés para varias funciones. El diagrama de visibilidad extendido: EVG Estrictamente hablando en el método originalmente propuesto se plantean dos posibilidades, computar la triangulación del dominio y luego el Grafo Extendido de Visibilidad ( EVG ), o computarlo todo de forma simultanea. En la implementación realizada se ha tomado la opción de procesarlo todo a la par para evitar almacenar información acerca de la triangulación y luego información extra acerca del propio Diagrama Extendido. Conceptualmente, una triangulación que respete las obstáculos que se plantean dentro del dominio, podría funcionar como un protodiagrama de visibilidad, proto− en el sentido de que existirían 16 Capítulo 2. Método de Ghosh & Mount muchos falsos negativos, pero funcionaría bajo ciertas situaciones concretas. La triangulación propuesta en el artículo original, [6], es una versión generalizada del método de Mehlhorn, [9], preparada para que triangule recintos con huecos en su interior. Figura 2.3 Ejemplo práctico de triangulación respetando obstáculos. Aquí ya podrían plantearse ciertas dudas respecto a la complejidad algorítmica propuesta en el artículo original, [6], donde se plantea que el método concreto es capaz de ejecutarse con una complejidad de O(E+nlogn) . Donde E es el número de lineas de visión en el grafo de visibilidad, que en el peor de los casos dado que el susodicho grafo puede ser trabajado como una matriz de dimensiones N x N , siendo N el número de puntos de paso, será una fracción ξ , que dependerá de la geometría de trabajo. Esto podría tender de forma aproximada como: O(ξN2)⋍O(N2).(2.2) Sumado a esto, hay que considerar que la complejidad algorítmica asociada a O(nlogn) está asociada a la computación de la triangulación de un entorno libre de obstáculos, es mencionado en el artículo que la triangulación que se consideró fue fundamentalmente similar a la de Mehlhorn, pero con una sencilla generalización, de la que no se aportan más detalles. Triangular un dominio sin agujeros es más sencillo que hacerlo con ellos, dado que se añade la necesidad de revisar si se está respetando adecuadamente la geometría de los susodichos agujeros. Este y otros aspectos serán analizados de forma crítica en la sección de revisión del método e implementación. El diagrama de visibilidad extendido plantea que las siguientes operaciones puedan ser realizadas en tiempo constante y con complejidad O(1) para un vértice v y una linea de visión (u,v) incidente en vdados: • Localización de las dos aristas pertenecientes al dominio triangulado incidentes en v . Que en la implementación se ha sacado más partido usando en su lugar Consecaristas , que no pertenecen al dominio triangulado sino que conforman las aristas que preceden y suceden a v . Se puede considerar que esta información está registrada en vectoru , véase la sección 2.3.2 para una explicación más detallada. •Los sucesores horarios y antihorarios de u,v. •Las extensiones horarias y antihorarias de u,v. •El reverso de (u,v). Que es trivial, dado que solo implicaría permutar las posiciones. Estas son las operaciones que se requiere que se puedan hacer sobre el grafo de visibilidad extendido, pero existen más datos que se almacenarán en el mismo, las familias de embudos y otras variables de interés precomputadas, que serán añadidas en la implementación en la sección 2.2.5. 2.2 Arranque del programa 17 Sucesores horarios y antihorarios Según la sección 7. Data Structure del artículo original, [6], se plantea que la implementación de la siguientes operaciones debe plantearse empleando una doble lista de adyacencia, tal que las entradas estén organizadas en el orden propuesto. Para simular este efecto, se han empleado estructuras vectoriales ordenadas para un punto dado, registrando así los nodos en un vector fila. Agrupando todos estos vectores fila en orden de punto de paso se obtiene la matriz Matrizorden. El proceso de ordenación se lleva a cabo al analizar cada punto de paso. Por esta razón, la matriz Matrizorden tendrá dimensiones N x N −1 , ya que un nodo no puede ordenarse angularmente en relación consigo mismo. Se utiliza un nodo como referencia, que corresponde al ángulo 0, ya sea el nodo 1 o 2 en el caso particular del nodo 1 , y se calculan los ángulos de ese nodo respecto a la referencia. Este ángulo se determina aplicando las propiedades fundamentales del producto escalar: a· b=|a|  barccos(θ),(2.3) para evitar las ambigüedades que comúnmente presenta la función inversa del coseno, es a adicionalmente comprobada la tercera componente del producto escalar. Todos estos valores angulares son ordenados usando el comando sortrows , dando como resultado el vector de nodos ordenados en sentido horario, que recorriendo de forma inversa el vector se conseguirá la versión en sentido antihorario. El uso de esta matriz para calcular sucesores será desarrollado en mayor detalle en el Apéndice A, sección A.3.2. Con el objetivo de no usar comandos find para buscar posiciones específicas de nodos dentro del vector fila, se usa la matriz Mat_ind_pos_enMatrizorden que contiene la información contraria, es decir en que posición del vector fila se encuentra el nodo de interés. Esto es interesante para calcular el sucesor horario o antihorario de un nodo en concreto, dado que si se quisiera saber cual es el sucesor dentro de la secuencia, sería necesario localizarlo y consultar el siguiente valor o el anterior en la secuencia. Es por ello que el disponer de esta información de antemano convierte el cálculo de O(N)en O(1). 207 %% Cálculo de las matrices de ordenación 208 % Esto es crucial para el desarrollo del algoritmo, se esecifica como dato 209 % del problema en la sección de estructura de datos. 210 %La matriz será NxN, donde el input es por filas, es decir, para una fila k 211 %dada, cada uno de los valores de esa fila son los otros nodos ordenados alrededor de 212 %el en sentido de las agujas del reloj visto desde el nodo k empezando por el primer nodo de waypoint (ya sorted). 213 Matrizorden=zeros(N,N-1); %N-1 por que el propio punto no se puede ordenar respecto de si mismo 214 filaaux=zeros(N,1); 215 216 for ii=1:N 217 %Formado elvector auxiliar a recorrer 218 if ii==1 219 filaaux=2:N; 220 elseif ii==N 221 filaaux=1:(N-1); 222 else 223 filaaux=[1:(ii-1),(ii+1):N]; 224 end 225 226 vbase=WPs(filaaux(1),:)-WPs(ii,:); 227 modvbase=norm(vbase);%Para calcular el ángulo 228 229 angulovec=zeros(2,N-1); 230 angulovec(2,:)=filaaux;%el primer nodo será el que corresponda al punto de base sobre el que medimos 231 232 for jj=2:(N-1)%Se ordenan los N-2 nodos restantes (Es decir, todos menos si mismo y el que usamos de referencia) 233 vpunto=WPs(filaaux(jj),:)-WPs(ii,:); 234 modvpunto=norm(vpunto); 235 anguloesc=acos(dot(vpunto,vbase)/(modvpunto*modvbase));%En radianes! 18 Capítulo 2. Método de Ghosh & Mount 236 237 %Discriminamos el sector para evitar ambigedades 238 qz=vpunto(2)*vbase(1)-vpunto(1)*vbase(2); % usamos el producto escalar para ver si queda a la izq o a la derecha del semiplano y asi discriminar el ángulo 239 %El producto vectorial es el de vbase sobre vpunto 240 241 if qz<0 %No puede ser igual a 0, implicaría que hay tres puntos en linea. 242 anguloesc=2*pi-anguloesc; 243 elseif qz==0 244 warning(’Tres puntos son colineares.’) 245 end 246 angulovec(1,jj)=anguloesc; 247 [180/pi*anguloesc, qz, filaaux(jj),vbase,vpunto]; 248 end 249 250 angulovec=sortrows(angulovec’,1)’;%Version ordenada de la matriz, la segunda fila son los vé rtices ordenados en sentido creciente de los ángulos 251 Matrizorden(ii,:)=angulovec(2,:); 252 end 253 254 255 %Matrizorden: Ordena los nodos en sentido horario partiendo desde el primer 256 %o el segundo nodo. 257 % Pongamos un ejemplo práctico, queremos saber cual es el nodo incidente en 258 % i que se encuentra a continuación de j en sentido horario: Eso se haría tal que Matrizorden(i,j +1). 259 % De forma análoga, el nodo incidente en i que se encuentra a continuación 260 % de j en sentido antihorario seria Matrizorden(i,j-1) 261 262 %% Matriz de correlación posición en Matrizorden con nodo 263 % Para evitar llamar varias veces al comando find dentro del algoritmo de 264 % cálculo de los CWs y los CCWs, precomputamos las posiciones de estas 265 % mismas, dado que esto es una cuestión geométrica. 266 267 Mat_ind_pos_enMatrizorden=nan(N,N); 268 269 for ii=1:length(Matrizorden(:,1)) 270 271 for jj=1:length(Matrizorden(1,:)) 272 valoraux=Matrizorden(ii,jj);%Da el valor del nodo 273 Mat_ind_pos_enMatrizorden(ii,valoraux)=jj;%Sería la operación inversa a matriz orden, en vez de decirnos la secuencia que tienen los nodos, dice la posición de la secuencia. 274 end 275 end 276 277 %Ejemplo de uso de Mat_ind_pos_enMatrizorden, queremos saber, que posición en Matrizorden tiene el 278 %nodo j para la secuencia de nodos alrededor de i, lo que debería hacerse 279 %es: posicion= Mat_ind_pos_enMatrizorden(i,j) Código 2.8 Cálculo de los sucesores CW/CCW. Nótese que ese postprocesado una vez más supera las expectativas de complejidad presentadas. Dado que se hacen dos bucles recorriendo el rango completo de N y al final de la segunda capa de bucle, se realiza una ordenación con el comando sortrows cuya complejidad algorítmica es, en el mejor de los casos, nlogn , es decir, la que correspondería a un Merging sort. A continuación se plantea la complejidad algorítmica de este cuello de botella presente en el preprocesado: O(N(N+NlogN)) ⋍O(N2logN).(2.4) Se adjunta en la figura 2.4, un ejemplo visual de cómo funciona la operación sucesores. Padres posibles En esta sección del código se revisa para cada combinación de punto de paso y polígono que nodos son factibles como padres, este concepto de paternidad cobrará sentido cuando se introduzcan los nodos en la sección 2.2.7. La función test_de_paternidadMain , cuya programación será desarrollada 2.2 Arranque del programa 19 v u CW(v,u) CCW(v,u) Figura 2.4 Sucesores CW y CCW. en el Apéndice A, sección A.3.2, revisa tres condiciones geométricas que validan si un nodo nodopropuesto perteneciente a un polígono tiene potencial para ser padre para un ancestro dado x : 1. Todos los nodos que pertenecen al polígono al que pertenezca x son potenciales padres. Es por ello que, 2. todos los nodos pueden ser padres para si mismos, dado que ellos son parte de la cadena de la que son ancestros. 3. Si nodopropuesto y x pertenecen a diferentes polígonos, solo podrán ser padres los que se encuentren en la cara oculta del polígono de nodopropuesto o bien la tangente al polígono ubicada más en el sentido de las agujas del reloj respecto a x. Todas estas comprobaciones se realizan en el código test_de_paternidadMain , en versiones tempranas del código esto se calculaba según se necesitaba disponer de él pero considerando que es una propiedad geométrica que no depende directamente del cálculo del grafo puede precomputarse. 281 %% Matriz padres posibles 282 EVG.MATposiblepadre=NaN(N,N);%Primera para nodox, segunda para el nodo a comprobar 283 EVG.WPs=WPs; 284 285 for jj=1:N 286 %Actualizar un nodo cambia todos los del polígono 287 288 for kk=[0 1:length(EVG.poligono2nodo(:,1))] 289 if kk==0 290 nodopropuestopadre=1; 291 else 292 nodopropuestopadre=EVG.poligono2nodo(kk,1); 293 end 294 295 [EVG]=test_de_paternidad_MAIN(jj,nodopropuestopadre,EVG); 296 end 297 end 298 299 EVG.MATposiblepadre(1,1)=1; 300 EVG.MATposiblepadre(2,2)=1; 301 EVG.MATposiblepadre(N-1,N-1)=1; 302 EVG.MATposiblepadre(N,N)=1; 303 EVG.MATposiblepadre(1,2)=1; Código 2.9 Identificación de padres posibles. 20 Capítulo 2. Método de Ghosh & Mount Extensiones horarias y antihorarias Esta operación no ha sido implementada directamente en el código, en la sección 2.2.7, se desarrolla una serie de criterios similares para discernir cuando es factible insertar nodos en familias, pero por resaltar las diferencias con el método original, se realizará la presentación de la operación. Para calcular la extensión horaria CX(u,v)es necesario seguir los siguientes pasos: 1. Se traza un arco de radio |uv|de 180 centrado en ven sentido horario. 2. Si la totalidad de ese arco pasa por el interior del dominio P , que es el conformado por la substracción de los agujeros poligonales al recinto exterior; la extensión será el siguiente segmento visible en sentido horario. Siempre que se mantenga todo el recorrido del arco dentro de P . En caso de que ese arco no estuviera dentro en su totalidad del dominio P , se considera la extensión no definida. La extensión antihoraria se extrae siguiendo los mismos pasos y restricciones salvo que trazando el arco en sentido antihorario. Se adjunta la imagen 2.5, con diferentes casos y aplicaciones para mayor claridad. Resaltadas en color verde se encuentran las extensiones factibles, y en roja un ejemplo de una extensión no definida. v1 u1 CCX (u,v) CX (u,v) v2 u2 Figura 2.5 Extensiones CX y CCX. 2.2.5 Formando la estructura de datos Este es todo el preprocesado que será realizado con la información proporcionada por la geometría. A continuación todas las variables computadas y otras que serán detalladas a continuación se añaden a la estructura de datos EV G. 304 %% Montando la estructura de datos 305 % Hay muchas variables precomputadas, para evitar entradas a funciones 306 % demasiado largas, las organizaremos todas en el Enhanced Visibility Graph 307 308 % EVG.WPs,Caminos_prohibidos,Matrizorden, Mat_ind_pos_enMatrizorden 309 310 Visibilidad=zeros(N,N); 311 312 EVG.Visibilidad=Visibilidad; 313 EVG.Caminos_prohibidos=Caminos_prohibidos; 314 EVG.Matrizorden=Matrizorden; 315 EVG.Mat_ind_pos_enMatrizorden=Mat_ind_pos_enMatrizorden; 316 EVG.TOT_famxv_glob={}; 2.2 Arranque del programa 21 317 EVG.TOT_famvy_glob={}; 318 EVG.nodosTOT_famvy_glob={}; 319 EVG.nodosTOT_famxv_glob={}; Código 2.10 Estableciendo el diagrama de visibilidad extendido EV G. 2.2.6 Antes de iniciar el bucle: Presentando la triangulación Una vez toda la información está preprocesada, es necesario presentar dos conceptos cruciales acerca del flujo de trabajo que tendrá el algoritmo implementado. La algoritmia entera se basa en los siguientes dos preceptos: 1. Expansión del dominio procesado mediante triangulación. 2. Determinación de visibilidad mediante recorrido de familias de embudos. Los detalles de la algoritmia de triangulación serán presentados junto con el código, pero siguiendo lo citado en la sección 2.2.4, se responderá al porqué plantear en un primer lugar un algoritmo de triangulación. Recordando la figura 2.3, existe cierta sinergia entre la triangulación de un dominio y la visibilidad entre diferentes puntos planteados en el mismo, por lo que es posible sacarle partido para estructurar el funcionamiento del código. El dominio P puede descomponerse en diferentes secciones disjuntas que conformen un todo: Pk=∪k i=1Tk,Pk=T1∪T2∪T3. .. ∪Tk−1∪TK.(2.5) Estas secciones Tk son las secciones trianguladas fruto de conectar el punto que se está añadiendo al grafo de visibilidad, que para ser congruentes con la implementación se denominará ii , con el dominio Pk−1 a través de los puntos de su envolvente convexa. Se adjunta la figura 2.6 para explicar el método de progresión aditiva. Se van visitando los diferentes puntos de paso en orden de coordenada xcreciente y se recoge esta información. T3 T4 T5 1 2 3 4 5 6 7 8 9 10 11 12 13 14 T6 T6 T7 T9 T9 T8 T11 T10 T13 T13 T14 T14 T14 Figura 2.6 Dominio triangulado compuesto. En la ilustración 2.6, se puede observar que algunas extensiones de dominio Tk pueden ser un espacio nulo, como es el caso de 1 o 2 , dado que no puede trazarse un área al añadirse un punto a un espacio nulo, y de manera acorde formar un área con dos puntos. Para el caso del nodo 12 , 28 Capítulo 2. Método de Ghosh & Mount xy ii Caso a) Validez de padres Padre válido Padre válido xy ii Caso b) Longitud de cadenas mínima Cadena válida Cadena no válida Conexión válida Conexión no válida xy ii Caso b) Longitud de cadenas mínima Figura 2.13 Casos ejemplificados de los criterios aplicados. x ii y familiaxv x y x' y' Figura 2.14 Reordenación de la familia partiendo de una base dada. 2.2 Arranque del programa 29 1function familiaembudo=calculafamilia(nodox,nodoy,EVG,ii) 2%El nodo más alto fue el último en computarse y es el que identifica si la 3%familia que se busca es xv o vy. Por ejemplo (15,14) es necesariamente una 4%familia computada en la iteración 15, de tipo vy 5if nodox>nodoy 6%Consulta a vy 7nodos_TOT=EVG.nodosTOT_famvy_glob{nodox}; 8else 9%Consulta a xv 10 nodos_TOT=EVG.nodosTOT_famxv_glob{nodoy}; 11 end 12 13 nodos_TOT(nodos_TOT(:,1)==nodox)=[]; 14 nodos_TOT(nodos_TOT(:,1)==nodoy)=[]; 15 16 if ~isempty(nodos_TOT) 17 posnodoy=EVG.Mat_ind_pos_enMatrizorden(nodox,nodoy); % Se saca la posición de nodoy en la secuencia de nodoorden 18 19 nodoorden=EVG.Matrizorden(nodox,:); % Se extrae el vector nodoorden que contiene los nodos ordenados en sentido CCW 20 nodoorden=[nodoorden(posnodoy+1:end),nodoorden(1:posnodoy)]; % Se ordena en CCW dejando primero a nodoy 21 nodoorden=nodoorden(end:-1:1); %Ahora y es primero en orden CW 22 [~, indicefiltrado]= ismember(nodoorden, nodos_TOT(:,1)); % Se filtran los valores posibles a partir de los nodos que pertenencian a subfamxv 23 nodos_TOT=nodoorden(indicefiltrado>0)’; % En este caso, se tienen todos los nodos ordenados respecto a ii 24 25 %nodos_TOT contiene ahora todos los nodos ordenados en un vector columna 26 27 %Se va a buscar descartar aquellos nodos que no sean factibles porque se 28 %encuentren por detrás de la familia 29 30 p2=EVG.WPs(nodox,:)’; 31 p3=EVG.WPs(nodoy,:)’; 32 v2=p3-p2; 33 34 35 for jj=1:length(nodos_TOT) 36 nodocheck=nodos_TOT(jj); 37 p1=EVG.WPs(nodocheck,:)’; 38 39 v1=p1-p2; 40 41 prodvec=v1(1)*v2(2)-v1(2)*v2(1); 42 43 if prodvec<0 %Esto significa que se encuentra en el semiplano válido de los nodos 44 if jj~=1 45 nodos_TOT=nodos_TOT(jj:end);%si el primer nodo es válido entonces, no hay necesidad de recortar los anteriores 46 end 47 break 48 elseif jj==length(nodos_TOT)&&prodvec>0 49 nodos_TOT=[];%ningun nodo es válido 50 end 51 end 52 %Se dispone en forma de familia 53 nodos_TOT=[nodos_TOT,nodos_TOT]; 54 end 55 56 %Se añaden los nodos x e y 57 nodos_TOT=[nodox,nodox;nodos_TOT;nodoy,nodox]; 58 59 %Se llama al procesado 60 [familiaembudo]=add_familia2(EVG,nodos_TOT); 61 end Código 2.11 Función calculafamilia. 30 Capítulo 2. Método de Ghosh & Mount 2.3 Bucle en ii Para realizar la triangulación del dominio de forma aditiva, véase la sección 2.2.6, es una necesidad recorrer todos los puntos de paso, añadiéndolos adecuadamente al grafo de visibilidad extendido. desde aquí, todo el código está metido en un gran bucle f or , que recorre todos los puntos de paso y se identifica el nodo que se procesa activamente con ii. 2.3.1 Nodos ii =1yii =2 El primer y el segundo punto de paso son los únicos que no pueden procesarse de manera normal, dado que no tienen área y se encuentran mirando enteramente fuera del dominio. Por lo que es necesario extender FNLextend y introducir los datos asociados a las familias a mano en EV G. Desde la perspectiva del punto de paso 1 , como las familias solo pueden incluir puntos de paso iguales o inferiores a ii dado que el proceso es aditivo, la familia está compuesta de un único embudo. Para 2 el único embudo posible, considerando que todos los puntos de paso estarán ubicados en valores superiores de coordenada x , es la familia 1,2 , que a su vez solo puede contenerse a si misma. 322 %% Loop general 323 % Por comodidad se procesan los primeros nodos por separado 324 325 ii=1; 326 EVG.TOT_famxv_glob{ii}=[1,1]; 327 EVG.TOT_famvy_glob{ii}=[1,1]; 328 EVG.nodosTOT_famvy_glob{1}=1; 329 EVG.nodosTOT_famxv_glob{1}=1; 330 331 ii=2; 332 EVG.TOT_famxv_glob{ii}=[1,1;2,1]; 333 EVG.TOT_famvy_glob{ii}=[2,2;1,2]; 334 EVG.nodosTOT_famvy_glob{2}=[1;2]; 335 EVG.nodosTOT_famxv_glob{2}=[1;2]; 336 337 338 EVG.Visibilidad(1,2)=1; 339 EVG.Visibilidad(2,1)=1; 340 341 %Para comprobar la visibilidad respecto a la envolvente 342 npol=length(EVG.poligono2nodo(:,1)); 343 vectorpol=[];%Inicialización de la variable Código 2.12 Procesado manual de los primeros puntos de paso. Por otro lado se inicializa la variable vectorpol , que permitirá ahorrar tiempo de procesado a la hora de determinar los puntos a visitar. 2.3.2 Ejecución normal del bucle: Triangulación A partir de este momento comienza el proceso de triangulación junto con el de determinación de visibilidad. Para ello se irá desarrollando conforme se muestra en el código, las decisiones que llevaron a la implementacíón final aquí desarrollado. Que se ha interpretado de manera más libre debido a que no se especifica en el documento de referencia [6], como se generaliza el método de Mehlhorn [9], para recintos con agujeros poligonales. Cabe destacar que es necesario terminar de preparar el terreno para el algoritmo de triangulación definiendo variables auxiliares, para ello se toma el periodo de ejecución del punto de paso 3 , del que se deducen sus propiedades siempre y cuando se respeten las directrices señaladas para el recinto externo (sección 2.2). 2.3 Bucle en ii 31 5%% Bucle en ejecución 6t1=toc; 7tic 8for ii=3:N 9%% %% Generación vector u 10 11 if ii==3 %No se puede sacar la envolvente de 2 puntos 12 13 %vectoru contiene la envolvente de puntos en la primera columna, vectoruextend en la segunda 14 %contiene los signos en formato + o - 1 15 vectoru=[1;2]; 16 vectoruextend=[1,1,1;2,1,1]; 17 % vectoruextend(ii,3) guarda la información de si es visible o no. 1 si, 0 18 % no 19 %vectorrecorre es el subfragmento a SPLITear 20 vectorrecorre=[1;2]; 21 22 %Valores postactualización 23 24 pol=EVG.nodo2poligono(3); 25 vectorpol=[vectorpol;pol];%añadimos a vectorpol el polígono 1, que corresponde por definición al que comienza en 3 26 vectoru=[1;3;2]; %Se añade 3 a la secuencia Código 2.13 Inicio de la triangulación. La variable vectoru conforma la envolvente de puntos que ya han sido procesado, nótese que no siempre va a contener todos los puntos procesados, solo el frente de aquellos que tengan los mayores valores de coordenada x . Por otro lado vectorpol , considera los polígonos que pudieran hacer pantalla a los puntos de paso que se procesen en futuras iteraciones. Por ejemplo, si el punto 4 está en el mismo polígono, sería necesario revisar que por cuestiones geométricas ninguna arista de ese polígono cortara la visibilidad con alguno de los puntos de la envolvente. Igual sería si 4 perteneciera a otro polígono, dado que sería el primer punto que se estuviera procesando sería imposible que el polígono al que pertenezca entorpeciera la linea de visión de la envolvente. Para ilustrar esta casuística se añade la figura 2.15. 1 2 3 3 44 Envolvente no convexa (vectoru) Contacto con envolvente válido Contacto con envolvente no válido 3 y 4 en el mismo polígono 1 2 3 y 4 en diferentes polígonos Figura 2.15 Ejemplo de procesamiento de vectorpol, bloqueo de visión respecto a la envolvente. 32 Capítulo 2. Método de Ghosh & Mount 368 else 369 370 %Atención: LOS VALORES DE U SE TOMAN DE LA ITERACIÓN ANTERIOR 371 vectoruextend=compruebasigno(vectoru,EVG,ii); %Aquí se le añadirá la segunda columna. Se revisan todos los valores Código 2.14 Llamada a compruebasigno desde el algoritmo de triangulación. En el caso de que el punto de paso sea superior a 3 , se llamará a la función compruebasigno . Esta función llama a una subrutina que añade una segunda columna con los valores de los signos del cociente de los productos vectoriales que surgen de la siguiente manera: se obtiene el vector que va desde ii hasta el punto anterior de vectoru respecto al que se está analizando, análogamente se procesa el vector posterior; cada uno de estos vectores se proyecta sobre el vector que que va de ii hacia el punto sobre el que se quiere calcular el signo. Este cociente si se normaliza respecto a si mismo, da −1 si el punto de estudio no es tangente y 1 si lo es. Estas proyecciones pueden verse de forma visual en la figura2.16 y su implementación en la caja de código 2.15. jj signo=-1 prodvec<0: jj-1 prodvec>0: jj+1 prodvec<0: jj+1 jj: signo=1 prodvec<0: jj-1 jj: signo=1 Figura 2.16 Ejemplo de la lógica de compruebasigno. 2.3 Bucle en ii 33 1%% Función comprueba signo 2% Comprueba como son los signos de los diferentes puntos a recorrer 3 4function vectoruextend=compruebasigno(vectoru,EVG,nodoii) 5vectoruextend=[vectoru,zeros(length(vectoru),1)]; 6vectoru=[vectoru(end);vectoru;vectoru(1)]; 7 8v1=zeros(2,1);%Vector que lleva de N a p1(punto anterior) 9v2=zeros(2,1);%Vector que lleva de N a p2(punto a evaluar) 10 v3=zeros(2,1);%Vector que lleva de N a p3(punto posterior) 11 12 p1=zeros(2,1);%Punto anterior 13 p2=zeros(2,1);%Punto a evaluar 14 p3=zeros(2,1);%Punto posterior 15 p4=EVG.WPs(nodoii,:)’;%Punto asociado a N 16 17 for jj=2:(length(vectoru)-1) 18 %Así evitamos que tener que utilizar casos específicos para los valores 19 %de 1 y 2 20 p1=EVG.WPs(vectoru(jj-1),:)’; 21 p2=EVG.WPs(vectoru(jj),:)’; 22 p3=EVG.WPs(vectoru(jj+1),:)’; 23 24 v1=p1-p4; 25 v2=p2-p4; 26 v3=p3-p4; 27 28 prodveca=v1(1)*v2(2)-v1(2)*v2(1); 29 prodvecb=v3(1)*v2(2)-v3(2)*v2(1); 30 31 signo=prodveca*prodvecb/abs(prodveca*prodvecb); 32 33 vectoruextend(jj-1,2)=signo; 34 end 35 36 end Código 2.15 Función compruebasigno. Todos los puntos de interés estarán ubicados entre dos tangentes, y se evaluará si se procesaran sus familias de embudos en función de sus visibilidad. Una vez se ha distinguido cuales son los puntos que potencialmente podrían resultar interesantes para analizar según la metodología de Ghosh & Mount, [6], es necesario comprobar si existen interrupciones de la visibilidad con este subconjunto de valores que se encuentran entre dos tangentes. Para ello se identifica cual es el primer y el último valor igual a 1 en el vectoru , usando dos comandos find , uno en sentido directo y otro inverso. Gracias a esto se puede evitar consultar el vector entero, debido a que la búsqueda para cuando se encuentra el único nodo de interés. 373 % Se computa el vector subu, solo hace falta revisar entre el primer 1(pos) y 374 % el último 1 (pos) 375 376 %% Creación vectorsubu 377 378 %El vector subu es el segmento del vector que potencialmente contiene los 379 %elementos a recorrer (si no hubiera aristas que cortan la visión vectorsubu sería igual a vector recorre) 380 381 ind1=find(vectoruextend(:,2)==1,1);% Primera aparición de 1 382 ind2=length(vectoruextend(:,2))-find(vectoruextend(end:-1:1,2)==1,1)+1; %Última aparición de 1 383 384 vectorsubu=vectoru(ind1:ind2); 385 vectorsubuextend=[vectoruextend(ind1:ind2,:)]; 386 387 %% Comprobación visibilidad 388 % Estas visibilidades se calculan revisando los conjuntos de polígonos y la 34 Capítulo 2. Método de Ghosh & Mount 389 % envolvente 390 391 aristaintersec=NaN(2,1); 392 poligono2nodoextend=[EVG.poligono2nodo;vectoru’,zeros(1,N-length(vectoru))];%comprobamos sobre la envolvente completa 393 vectorpolextend=[length(poligono2nodoextend(:,1));vectorpol]; Código 2.16 Identificación de polígonos de interés en la triangulación. Para comprobar si existen cortes de visibilidad con los puntos concretos de la envolvente no convexa, se seguirá el siguiente procedimiento: 1. Definición de variables de utilidad en el bucle: •aristaintersec es una variable que almacena la última arista que cortó la visibilidad para el punto anterior, una vez se inicia el bucle esta información todavía no existe. •poligono2nodoextend : es una versión extendida de polgono2nodo que contiene además, información acerca del polígono que conforma la envolvente no convexa, que podría causar intersecciones en casos concretos, dado que su perímetro puede desarrollar dientes de sierra. •vectorpolextend : que contiene la misma información que vectorpol , acerca de que polígonos son susceptibles de interrumpir la visibilidad, pero añadiendo el polígono n+1 , que sería el de la envolvente. Como este cambia para cada iteración, es necesario actualizarlo. 2. Se recorren todos los valores de vectorsubu buscando cortes, comprobando primero si aristaintersec sigue bloqueando la visibilidad. inter funciona como una bandera, en el momento en el que se detecta una intersección se levanta, inter =1 , y se registra acordemente que no es visible el punto de estudio y se pasa al siguiente. Esta información acerca de la visibilidad se guarda en la respectiva fila de la tercera columna vectorsubuextend(j j,3) . Si por un casual, la conexión entre ii y el punto de estudio es un Camino_prohibido , véase la sección 2.2.3 para más detalle, se interpreta que hay una intersección dado que la visibilidad esta bloqueada pero no por una intersección, sino por producirse este intento de visión dentro de uno de los agujeros. a) Para cada valor de vectorsubu que no se haya podido descartar de forma preliminar, es necesario comprobar si hay intersecciones con los polígonos registrados en vectorpol , se recorre cada uno de ellos buscando una intersección. b) Las aristas a estudiar están registradas de manera conveniente en poligono2nodoextend , por lo que eliminados los valores nulos residuales, para cada pareja de vértices consecutivos, que conforman una línea, se busca si existe intersección con la recta que va desde ii hasta el nodo de estudio. c) Este corte con las aristas de los polígonos solo puede darse si se cumplen las siguientes condiciones: • Al menos uno de los vértices que compone la arista debe tener preceder a ii . Si están más allá de ii no puede producirse un corte. • No se considerará ninguna de las aristas que tengan a ii o al nodo de vectorsubu que se estudie en ese momento. Dado que solo puede existir una intersección entre dos rectas, y ese punto sería o bien ii o el nodo de estudio. d) Se guardan las coordenadas asociadas a los puntos de paso de la arista que se vaya a estudiar line1, y la recta que une ii y el punto de estudio, line2. 2.3 Bucle en ii 35 e) Se calcula si existe un punto de intersección usando geometría euclidiana básica, que será desarrolla y mostrada en el correspondiente Apéndice A, sección A.3.2. f) En el caso de existir se corta la revisión de aristas gracias al levantamiento de la variable bandera inter y se paran de procesar las aristas gracias a un break. 3. Por defecto, no existe corte hasta que se detecta uno, es por ello que vectorsubuextend(j j,3) = 1es establecido al principio del bucle y cuando se levanta la bandera se cambia a o. 395 for jj=1:length(vectorsubu) 396 % ’siguiente u’ 397 inter=0; 398 contadorpol=1; 399 vectoraristaextend=zeros(1,length(poligono2nodoextend(1,:))); 400 401 %Antes de empezar a calcular para cada polígono debería comprobarse si 402 %aristaintersec sigue bloqueando la visión del nuevo punto 403 404 %% Comprobamos visibilidad con arista intersec 405 406 if ~isnan(aristaintersec(1,1)) && aristaintersec(1)~=ii && aristaintersec(2)~=ii && aristaintersec(1)~=vectorsubu(jj) && aristaintersec(2)~=vectorsubu(jj) 407 %Si no es nan calculalo, si no lo es, compruebalo normal 408 409 %line1 es la línea que une N con el punto de la secuencia a comprobar 410 %line2 es la línea con la que se busca que interseque (cualquiera de las aristas de los polígonosconsiderados) 411 line1=[EVG.WPs(ii,1),EVG.WPs(vectorsubu(jj),1);... 412 EVG.WPs(ii,2),EVG.WPs(vectorsubu(jj),2)]; 413 line2=[EVG.WPs(aristaintersec(1),1),EVG.WPs(aristaintersec(2),1);... 414 EVG.WPs(aristaintersec(1),2),EVG.WPs(aristaintersec(2),2)]; 415 416 [inter,~]=cortapuntos(line1,line2); 417 418 if inter==1 419 %si hay intersección se registra y arista intersec 420 %se queda igual. 421 vectorsubuextend(jj,3)=0; 422 % ’CorteAristaintersec’ 423 end 424 end 425 426 if Caminos_prohibidos(ii,vectorsubu(jj))==1 427 inter=1; 428 vectorsubuextend(jj,3)=0; 429 % ’Camino prohibido’ 430 end 431 432 %% Se comprueban cortes con cada uno de los polígonos candidato 433 while inter==0 && contadorpol<=length(vectorpolextend(:,1)) % En el momento que se identifica un corte la línea N-vectorsubu(jj) es no visible 434 % ’comprobación pol’ 435 436 poligono=vectorpolextend(contadorpol); 437 vectoraristaextend=poligono2nodoextend(poligono,:);%Extend incluye para poligono2 nodo el vectoru, para arista los 0 sobrantes 438 naristaspol=find(vectoraristaextend(1:1:end)==0,1)-1; 439 vectorarista=[vectoraristaextend(naristaspol),vectoraristaextend(1:naristaspol), vectoraristaextend(1)]; 440 %Procesarlo así hace mucho más sencillo gestionar las parejas a 441 %procesar 442 vectorsubuextend(jj,3)=1; 443 444 % Para cada polígono se comprueban los cortes con cada arista 445 for kk=2:(naristaspol+2) 446 if (vectorarista(kk-1)<ii || vectorarista(kk)<ii) && vectorarista(kk-1)~=ii && vectorarista(kk)~=ii && vectorarista(kk-1)~=vectorsubu(jj) && vectorarista(kk)~= vectorsubu(jj)%Las dos últimas condiciones sirven para que no se compruebe si hay intersecci ón con las aristas que pasan por el punto que vamos a añadir. Esto es así porque daría un 36 Capítulo 2. Método de Ghosh & Mount falso positivo diciendo que existe intersección entre las aristas que conforman el punto que vamos añadir, pero esas no bloquean la visibilidad, serían otras del polígono. 447 %line1 es la línea que une N con el punto de la secuencia a comprobar 448 %line2 es la línea con la que se busca que interseque (cualquiera de las aristas de los polígonosconsiderados) 449 %Las lineas vendrán dadas en el formato 450 % [x1, x2 451 % ;y1,y2] 452 453 line1=[EVG.WPs(ii,1),EVG.WPs(vectorsubu(jj),1);... 454 EVG.WPs(ii,2),EVG.WPs(vectorsubu(jj),2)]; 455 % linea1=[ii,vectorsubu(jj)] 456 line2=[EVG.WPs(vectorarista(kk-1),1),EVG.WPs(vectorarista(kk),1);... 457 EVG.WPs(vectorarista(kk-1),2),EVG.WPs(vectorarista(kk),2)]; 458 % linea2=[vectorarista(kk-1),vectorarista(kk)] 459 [inter,~]=cortapuntos(line1,line2); 460 if inter==1 461 % ’Corte’ 462 % Guardamos la información de la arista que cortó 463 % para empezar ahí la siguiente iteración 464 aristaintersec=[vectorarista(kk-1),vectorarista(kk)]; 465 vectorsubuextend(jj,3)=0; 466 break 467 end 468 end 469 end 470 contadorpol=contadorpol+1;% cambio de polígono 471 end 472 inter=0;%reseteo para comprobar el siguiente valor de u 473 end Código 2.17 Cálculo de puntos de corte en el área visible de la envolvente. Una vez obtenida la información acerca de que puntos están entre dos tangentes y además no encuentran su visión obstaculizada por los polígonos potenciales ni la propia envolvente no convexa, se computa el vectorrecorre. 477 vectorrecorre=vectorsubuextend(vectorsubuextend(:,3)==1); Código 2.18 Identificación de vectorrecorre. Este vectorrecorre tiene su contraparte en el documento original de Ghosh & Mount en el Lema 2.1, [6], del que deriva la siguiente casuística: •vectorrecorre contiene dos puntos: Es decir, solo esa línea de visión puede servir como cono para conectar el punto ii al dominio triangulado extendido Pk. •vectorrecorre contiene más de dos puntos: El dominio se puede ampliar añadiendo más de un triángulo, pero existen ahora diferentes conos de visión a comprobar. La información de estos conos de visión deberá luego sintetizarse. Este caso y el de dos puntos son los más comunes. •vectorrecorre solo tiene un valor: Esto significa que al no existir cono de visión, el único contacto posible es el de ii con el único valor perteneciente al vector. Esto sucede para cuando el punto de paso está a la vuelta de la esquina, y su única conexión es otro punto anterior, por lo general del mismo polígono, que bloquea la visibilidad. •vectorrecorre está vacío: Esto sólo sucede en recintos no convexos, cuando existen dos puntos de paso que conectan a un tercero entre los dos, que tiene menor coordenada x que los anteriores. Estas aristas laterales actúan como una pantalla y hacen imposible cualquier tipo de conexión con la envolvente. 2.3 Bucle en ii 37 Esta casuística puede verse representada para un mismo dominio para diferentes valores de punto de paso en la figura 2.17. Caso a): Dos puntos Caso b): Más de dos puntos Caso c): Un único punto Caso d): Ningún punto ii ii ii ii Figura 2.17 Diferentes casos para vectorrecorre. Cada caso influye a su vez con la manera de introducir ii en el vectoru de cara a la siguiente iteración: •vectorrecorre contiene dos puntos o más: Se localiza en vectoru la ubicación de los dos nodos extremos de vectorrecorre , se introduce ii sobrescribiendo los nodos que hubiera entre medias. •vectorrecorre solo tiene un valor: Dado que no se tiene referencia de entre que valores hay que insertar el nuevo nodo ii , debido a que solo hay un valor, es necesario establecer algún tipo de regla. Se comprueba cuales son los nodos consecutivos asociados al nodo de vectorrecorre y se insertan en orden de y creciente. Distinguiendo dos casos: xv si ii tiene mayor coordenada yyvy en caso contrario. •vectorrecorre está vacío: Se recorre el vectoru , insertando ii cuando se detecte la primera posición en la que su posible predecesor y su posible predecesor tengan menor y mayor coordenada y respectivamente. Esto es debido a que se carece completamente de referencias de vectorrecorre acerca de donde debería conectar. 479 %% Actualización vectoru 480 % Añadiendo valores a vectoru 481 482 ind1=find(vectorsubuextend(:,3)==1,1); 483 if ~isempty(ind1) 484 nodo1=vectorsubuextend(ind1); 485 %Dónde se ubica el valor nodo1 en vectoru? 486 ind1=find(vectoru(:)==nodo1,1); 487 44 Capítulo 2. Método de Ghosh & Mount xy ii ab xy ii padre en el árbol inferior padre en el árbol superior xy ii Dentro del cono: Visible Dentro del cono: Visible Fuera del cono: no Visible Figura 2.20 Cálculo de la visibilidad mediante el empleo de las características de los embudos. 1. Embudos que se encuentran fuera del cono de visión formado por x−ii eii −y. 2. Embudos cuya visibilidad está bloqueada por un tercer nodo en el ala izquierda (aprovechando el conocimiento acerca del padre en el árbol inferior). 3. Embudos cuya visibilidad está bloqueada por un tercer nodo en el ala derecha (aprovechando el conocimiento acerca del padre en el árbol superior). En teoría, esto debería funcionar, y de forma general lo hace con resultados satisfactorios, el problema radica en que existen ciertas configuraciones en las que lo anteriormente explicado no aplica. Este hecho se ha descubierto durante las labores de prueba, que en el caso de que un embudo pudiera tener dos padres potenciales, existirían dos conos posibles, y probar ambas alternativas no sería factible dado que solo se puede identificar un padre. El ejemplo de la figura 2.21 lo muestra sin lugar a dudas, dos configuraciones diferentes de embudos dan lugar a dos rangos de visibilidad muy diferentes. Dado que en el documento original no se dan indicaciones acerca de cómo seleccionar o reorganizar embudos al cambiar de familias, se ha optado por resolver el problema por fuerza bruta y comprobar si existen cortes con las aristas de los polígonos asociados a los nodos presentes en la familia, y en última instancia a todos los nodos asociados a los polígonos que se visitan. Todo esto queda desarrollado en el Apéndice A, sección A.3.2. Nótese además, que por como están configurados los SPLIT , al basarse estos en un proceso de triangulación incremental, las visibilidades en el grafo se añaden desde el punto de estudio ii a todos aquellos que tengan menor coordenada x . Esto se hacía de forma similar en el método de Lee, que computaba los semiplanos superiores al punto de estudio, salvo que en este caso son los semiplanos izquierdos. La visibilidad hacia valores mayores de x se computará cuando se alcancen estos puntos y se detecte su visibilidad, dado que el grafo de visibilidad es simétrico. Cálculo de u Dentro del documento de Ghosh & Mount[6], se especifica que u es un nodo visible desde ii y que precede a otro nodo, sin dar detalles de su cálculo. Con tal de trazar cierta simetría en el problema se han tomado propiedades y características de t para poder usar de forma similar sus propiedades en el árbol izquierdo. Tampoco especifica directamente como identificar el nodou , y parte de la base de que ya está localizado. Con las propiedades asignadas de t, esto puede lograrse. • Se recorre de uno en uno cada embudo de la familia, computando si son visibles para ii o no. Y en el caso de no serlo, se revisa la forma en la que no son visibles, gracias a modosalida . A partir de aquí existen dos posibilidades: que sea visible o que no. 2.3 Bucle en ii 45 x ii y ab Figura 2.21 Ejemplo de fallo en la demostración geométrica de la visibilidad con las propiedades de los embudos. • Si el embudo es visible: Si las condiciones son apropiadas, se busca aplicar SPLIT al embudo actual de la operación y el anterior embudo que fuera visible, se almacenan los embudos obtenidos. • Si no es visible: En el caso de que no sea visible por ′f ueraconoi′ , es decir, por un corte de visibilidad, producido por un nodo del ala izquierda, se almacena el embudo en f amiliavy. • Se comprueba si el bucle debería parar, ya sea si es a causa de que se ha alcanzado nodoy o si es porque el embudo de estudio no cumple con la condición ′f ueraconoi′. 740 %% Calculo u 741 % u es el último nodo que puede bloquear la visibilidad de nodos que se encuentren en el ala izquierda de la familia 742 % ’Calculo u’ 743 744 745 indiceaux=1; 746 nodoaux=familiaembudo(indiceaux,1); 747 748 EVG.Visibilidad(nodoaux,ii)=1; 749 EVG.Visibilidad(ii,nodoaux)=1; 750 751 familiaxv=[familiaxv;familiaembudo(indiceaux,:)]; 752 familiavy=[familiavy;familiaembudo(indiceaux,1),ii]; %Al partir el embudo queda la parte de xv igual y la de vy, se convierte en un nodo cuya base es N 753 754 modosalidaextend=strings(length(familiaembudo(:,1)),1); 755 modosalidaextend(1)=’VIS’; 756 757 %El primer nodo est añadido por defecto 758 flag=0; 759 indiceaux=2; Código 2.24 Implementación del cálculo de u: Arranque. Antes de empezar el bucle es necesario realizar cierta preparación, se añade x−ii al diagrama de visibilidad como visible y el correspondiente embudo a f amiliaxv y f amiliavy . Se define modosalidaextend , que guarda los out puts de texto del algoritmo de detección de visibilidad, 46 Capítulo 2. Método de Ghosh & Mount visibleDEF , en el vector para realizar de forma sencilla operaciones dentro del bucle, como pararlo o realizar SPLIT de forma recursiva. Se define la bandera f lag =0 que mantendrá el bucle ejecutándose hasta alcanzar la condición de parada, ya sea alcanzar el final de la lista, o encontrar el último embudo del ala izquierda que bloquea la visibilidad de los siguientes, es decir, el nodou. 761 while flag==0 762 %Arrancamos el bucle 763 [VIS,modosalida,EVG]=visible_DEF(EVG,familiaembudo,indiceaux,ii,nodox,nodoy); 764 modosalidaextend(indiceaux)=modosalida; 765 766 if VIS==true 767 nodoaux=familiaembudo(indiceaux,1); 768 EVG.Visibilidad(nodoaux,ii)=1; 769 EVG.Visibilidad(ii,nodoaux)=1; 770 %Si es visible se comprueba que se pueda splitear 771 772 for jj=(indiceaux-1):-1:1 773 if strcmp(’VIS’,modosalidaextend(jj)) 774 nodoa=familiaembudo(jj,1); 775 break 776 end 777 end 778 779 nodob=familiaembudo(indiceaux,1); 780 781 %Una vez se tiene la base de SPLIT, se comprueba que esta sea 782 %congruente 783 784 %Se hace una necesidad técnica el comprobar que la forma del SPLIT 785 %es adecuada, nodob debe preceder a nodoa CW 786 p1=EVG.WPs(nodoa,:)’; 787 p2=EVG.WPs(nodob,:)’; 788 p3=EVG.WPs(ii,:)’; 789 v1=p1-p3; 790 v2=p2-p3; 791 prodvec=v1(1)*v2(2)-v1(2)*v2(1); 792 793 if prodvec<0 794 llamadaCW=true; 795 else 796 llamadaCW=false; 797 end 798 799 % CONDVACIO 800 %Para evitar que nodos que se postcedan en la familia no estén 801 %entorpeciendo y cumplir la condición de que los nodos deban 802 %encontrarse vacíos, introducimos una condición adicional (4,6,caso1; check nodo34) 803 804 %para ello se organizan todos los nodos alrededor de ii partiendo 805 %desde nodoy, se revisan si entre los dos nodos propuestos hay uno o 806 %más nodos 807 808 posnodoy=EVG.Mat_ind_pos_enMatrizorden(ii,nodoy); % Se saca la posición de N en la secuencia de nodoorden 809 nodoorden=EVG.Matrizorden(ii,:); % Se extrae el vector nodoorden que contiene los nodos ordenados en sentido CCW 810 nodoorden=[nodoorden(posnodoy+1:end),nodoorden(1:posnodoy)]; % Se ordena en CCW dejando último a N 811 nodoorden=nodoorden(end:-1:1); %A hora N es primero en orden CW 812 813 [~, indicefiltrado]= ismember(nodoorden, familiaembudo(:,1)); % Se filtran los valores posibles a partir de los nodos que pertenencian a subfamxv 814 familiaord=nodoorden(indicefiltrado>0); % En este caso, se tienen todos los nodos ordenados respecto a ii 815 816 ind1=find(familiaord==nodoa); 817 ind2=find(familiaord==nodob); 818 819 condvacio=[]; 2.3 Bucle en ii 47 820 if ind2-ind1>1 821 for jj=(ind1+1):(ind2-1) 822 %Se comprueba si se encuentran DENTRO del cono 823 nodocheck=familiaord(jj); 824 isInside = checkIfInsideTriangle(ii, nodoa, nodob, nodocheck, EVG); 825 if isInside==true 826 condvacio=false; 827 break 828 end 829 end 830 if isempty(condvacio) 831 condvacio=true; 832 end 833 834 else 835 condvacio=true; 836 end 837 838 839 840 if ~(nodox==nodoa && nodoy==nodob) && loop<=maxloop && llamadaCW && condvacio 841 %Para evitar recursividades, es necesario meter esta 842 %condición, tmbn es necesario integrar en familiaguia estás familias de embudos, si no se pierden las posibles visibilidades 843 844 [EVG,familiaxv_sub,familiavy_sub]=SPLIT(nodoa,nodob,EVG,ii,loop,’interno’); 845 846 %El primer valor ya está computado en familiaxv, el ultimo es 847 %ii por lo que no se añade 848 vecaux=familiaxv_sub(2:end-1,:); 849 850 [familiaxv]=add_familiaxv(EVG,familiaembudo,familiaxv,vecaux,nodox); 851 852 vecaux=familiavy_sub(3:end,:); % el primer valor es ii, no se añade, el segundo es el nodoy del anterior, por lo que tampoco 853 [familiavy]=add_familiavy(EVG,familiaembudo,familiavy,vecaux,ii); 854 else 855 [familiaxv]=add_familiaxv(EVG,familiaembudo,familiaxv,familiaembudo(indiceaux,:), nodox); 856 [familiavy]=add_familiavy(EVG,familiaembudo,familiavy,familiaembudo(indiceaux,:),ii); 857 end 858 end Código 2.25 Implementación del cálculo de u:nodoaux es visible. Cómo se indicó en el resumen inicial, lo primero que se hace es revisar la visibilidad y almacenar el modosalida . Dentro del i f que comprueba si la variable VIS contiene información acerca de si es visible o no, se guarda en el diagrama de visibilidad. A continuación se recorre hacia atrás los nodos ya visitados buscando el último nodo que fue visible. Ese nodo se identificará como el nodox de la siguiente generación, que se llamará para evitar confusiones nodoa , de forma análoga, el nodo de estudio que sería el nodoy en la recursividad se llamará nodob . Se adjunta la imagen 2.22, donde se puede observar como a lo largo de la búsqueda de nodou se van visitando diferentes nodos del ala izquierda y se realizan operaciones de SPLIT acorde, esto actúa como un seguro en caso de que, por ejemplo, no se hayan guardado debidamente embudos en la familia, y pudiera resultar altamente beneficioso en el caso del SPLIT 3de la imagen. Antes de ejecutar este SPLIT es necesario garantizar que nodoa preceda a nodob tanto en orden de embudo, cosa que está inherente al concepto de familia, como en sentido horario pivotando desde ii . Esto es algo que no se desarrolla en el documento original, pero existen casos concretos en los que puedan darse dos nodos visibles atravesando túneles formados entre polígonos, que provocarían fallos. Esto además actúa como un seguro en el caso de que un embudo tuviera algún fallo en su asignación en la familia. Todo esto se recoge en la definición de la variable llamadaCW. Además se comprueba si entre nodoa y nodob pudiera haber algún nodo dentro del cono, que provocaría falsos positivos, dado que se trabaja suponiendo que estos conos de visión están vacíos. 48 Capítulo 2. Método de Ghosh & Mount ii xy u a b 1 2 34 Figura 2.22 Uso de recursividad SPLIT para detección de embudos. Para ello se extrae de Matrizorden , véase 2.2.4, se busca la secuencia de nodos pivotados en torno a ii a partir de nodoy en el sentido de las agujas del reloj, y se verifica que entre los dos no existan nodos, y de existir, se comprueba si quedan dentro del cono de visión gracias a la función checkI f InsideTriangle . Esta lo que hace es una operación de cálculo de áreas y comprueba como se distinguen entre estas áreas entre sí, si el total de las secciones suma igual o menor que el cono, el punto se encuentra dentro. Una explicación esquemática puede encontrarse en 2.23. Si ningún nodo se encuentra dentro, no se ha recurrido más allá del número permitido máximo de bucles maxloop y no se está intentando partir el mismo SPLIT que el actual, se procede a ejercer la recursividad. a b x ii a b ii ii y Obstáculo Área total a b ii Dentro a b x ii a b ii a b a b ii y Obstáculo Área total Fuera nodo nodo nodo nodo nodo Figura 2.23 Uso de checkI f InsideTriangle para comprobar la posición relativa de un nodo externo. 2.3 Bucle en ii 49 Si se llega a realizar la llamada a SPLIT se extrae la familia partida, como se esperaba, y esas familias se añaden a f amiliaxv y a f amiliavy usando las funciones análogas add_f amiliaxv y add_f amiliavy . En el código 2.26, se puede identificar una estructura sencilla consistente en identificar nodos repetidos en vecaux que estén a su vez en f amiliaxv , a partir de ahí los nodos no repetidos se añaden. Se organizarán fuera de la llamada de SPLIT principal. Las funciones para f amiliaxv yf amiliavy son completamente iguales por lo que por compacidad solo se añade una. En el caso de que no se cumplan todas las condiciones que permitan hacer la llamada a SPLIT solo se añaden los nodos a y b a las correspondientes familias usando las ya mencionadas funciones. 1function [familiaxv]=add_familiaxv(EVG,familiaembudo,familiaxv,vecaux,nodox) 2 3% Obtener los índices de los elementos de vecaux(:,1) que NO están en familiaxv(:,1) 4[ismem, ~] = ismember(vecaux(:,1), familiaxv(:,1)); 5 6% Invertir el valor lógico para obtener los índices NO repetidos 7indicesnorep = find(~ismem); % Aquí mantenemos el orden de vecaux 8 9if ~isempty(indicesnorep) 10 % Mantener solo los elementos de vecaux que no están en familiaxv 11 vecaux = vecaux(indicesnorep, :); 12 else 13 vecaux = []; 14 end 15 16 % Mantener solo los valores únicos en familiaxv respetando el orden original 17 18 % ’Inicio’ 19 if ~isempty(vecaux) 20 familiaxv=[familiaxv;vecaux]; 21 end 22 23 end Código 2.26 Implementación de add_f amiliaxv. La última parte del código presenta una comprobación del modo de salida, para ubicar correctamente el nodo en cuestión en la correspondiente familia. En el caso ′f ueraconoi1′ , es necesario añadir el embudo a la familia vy , si no fue visible para ii por un obstáculo que interrumpía su visión por la izquierda, tampoco lo será para la base xv . Esto queda reflejado en la figura 2.24, dónde se puede observar la estrategia de determinación de visibilidad mostrada en la figura 2.20, quedando constancia visual de lo anteriormente explicado. Es por ello que hace una revisión rápida del caso y se añade si procede. En el caso de que el modo de salida fuera otro, no sería necesario añadirlo ya que en etapas posteriores esto será contemplado, por lo que no se perderá esta información. Para decidir si continuar con el bucle o no se hace de nuevo una distinción de casos: 1. Nodo de estudio visible o ′f ueraconoi1′ : Se revisa si se ha alcanzado el último embudo de la familia, si es el caso se corta el bucle. Si no, se incrementa indiceaux y se continua con el escrutinio. Se levanta una bandera llamada resuelta f amilia , para indicarle al resto del código si es necesario seguir aplicando la metodología implementada. 2. Nodo de estudio no visible en casos diferentes a ′f ueraconoi1′ : El bucle se corta, el último nodo registrado como visible se identifica como nodou . La familia no está entera resuelta, por lo que es necesario continuar el procesado. 860 if strcmp(’fueraconoi1’,modosalida) 861 [familiavy]=add_familiavy(EVG,familiaembudo,familiavy,familiaembudo(indiceaux,:),ii); 862 end 863 50 Capítulo 2. Método de Ghosh & Mount xy ii nodo NO: familiaxv SI: familiavy Figura 2.24 Justificación de añadir a f amiliavy pero no a f amiliaxv: Caso ′f ueraconoi1′. 864 865 if strcmp(’VIS’,modosalida) || strcmp(’fueraconoi1’,modosalida) 866 %Cualquiera de esos 2 casos continua corriendo el bucle en 867 %búsqueda de $u 868 if indiceaux<length(familiaembudo(:,1)) 869 indiceaux=indiceaux+1; 870 nodoaux=familiaembudo(indiceaux,1); 871 else 872 flag=1; %Con esto lo que se hace es salir del bucle una vez se han computado todos los nodos 873 indicesol=indiceaux; 874 end 875 else 876 flag=1;%se rompe el bucle 877 for jj=(indiceaux-1):-1:1 878 if strcmp(modosalidaextend(jj),’VIS’) 879 indicesol=jj; 880 break 881 end 882 end 883 end 884 end 885 886 nodou=familiaembudo(indicesol,1); 887 indiceu=indicesol; 888 889 resueltafamilia=false; 890 891 if indiceu==length(familiaembudo(:,2)) % este es el caso particular de que se han computado en el calculo de u todos los embudos 892 resueltafamilia=true; 893 end Código 2.27 Implementación del cálculo de u:nodoaux no visible y cierre. Calculo de u′ Una vez identificado u , si la familia no está resuelta, el código se bifurca gracias a una sentencia i f en base a la variable resuelta f amilia . El objetivo de esta parte del código es identificar u′ , el hijo de u ubicado más en el sentido de las agujas del reloj cuya visibilidad de ii está interrumpida por u , para que sirva como ayuda para localizar r , el primer nodo en orden de embudo del ala derecha, 2.3 Bucle en ii 51 siendo el más lejano de esta ala. Tampoco se proporcionan indicaciones en el artículo original acerca de su cálculo pero gracias a la claridad de la definición su localización es sencilla. En la primeras primeras líneas se observa una vez más la preparación para el arranque del bucle, y el inicio de la sentencia i f que realiza la derivación conforme a la naturaleza de f amiliaresuelta . Se recorren desde nodoy hasta nodou todos los nodos, aumentando contadorindice y restándolo al índice auxiliar que realiza el seguimiento de la posición en el bucle indicecandidatouprima . Si el padre es nodou , comprobación trivial si se comprueba la segunda componente del vector fila correspondiente que guarda el embudo; se comprueba si la componente qz resultante del cálculo de v1×v2 , donde v1 es el vector de origen ii y destino nodou , mientras que v2 comparte origen pero su destino es el nodo candidato a ser uprima , es mayor o igual a 0 . En el caso de ser mayor o igual a 0 , nodou se encontraría más adelantado en el sentido de las agujas del reloj, por lo que bloquearía su visión, y por ende, dado que se está recorriendo la familia de fin a principio, se obtendría el primer nodo cuyo padre es nodou y su línea de visión con ii se ve interrumpida por este, figura 2.25. x ii y u u' v1 v2 Figura 2.25 Ejemplo de identificación uprima. 895 %% Cálculo u’ 896 % ’Cálculo uprima’ 897 % uprima es el hijo de u más clockwise (por definición de u, no es visible uprima desde v), puede no existir 898 % Se recorrerá usando un bucle que vaya desde y hasta u, si no encuentra un 899 % solo hijo de u entonces saldremos del bucle 900 901 contadorindice=0; 902 903 indicecandidatouprima=(length(familiaembudo(:,2))-contadorindice); 904 nodouprima=NaN; 905 indiceuprima=NaN; 906 907 if resueltafamilia==false 908 909 910 while indicecandidatouprima>indiceu 911 % Solo realizamos este bucle cuando u sea diferente de y 912 if familiaembudo(indicecandidatouprima,2)==nodou 913 p1=EVG.WPs(familiaembudo(indiceu,1),:); 914 p2=EVG.WPs(familiaembudo(indicecandidatouprima,1),:); 52 Capítulo 2. Método de Ghosh & Mount 915 p3=EVG.WPs(ii,:); 916 917 v1=p1-p3; 918 v2=p2-p3; 919 920 prodvec=v1(1)*v2(2)-v1(2)*v2(1); 921 if prodvec>0 922 %El padre del candidato nodouprima es u por lo que no es candidato, es 923 %el nodo 924 indiceuprima=indicecandidatouprima; 925 nodouprima=familiaembudo(indiceuprima,1); 926 break 927 end 928 end 929 contadorindice=contadorindice+1; 930 indicecandidatouprima=(length(familiaembudo(:,2))-contadorindice); 931 end Código 2.28 Cálculo de u′. Cálculo de qyr La identificación de estos nodos si se encuentra documentada en el documento original de Ghosh & Mount, [6], y para ello es necesario distinguir dos casos: 1. nodouprima no existe: No hay entonces en la familia ningún embudo cuya visibilidad es bloqueada por nodou , lo que implica que escondido por nodou para ii no hay ningún nodo, por lo que la definición de q (el último nodo del ala izquierda) queda invalidada, y se realiza la siguiente asignación nodoq =nodou . El siguiente nodo entonces en orden de embudo es nodor (el último nodo del ala derecha). 2. nodouprima existe: Que será el caso que se explicará con más detalle a continuación, dado que es necesario aplicar las propiedades de la familia de embudos para encontrarla. En el caso de existir uprima , hay que considerar que este mismo embudo, todos sus hermanos CCW , es decir, los hijos de nodou que se encuentran posicionados más en el sentido de las agujas del reloj respecto a nodouprima , y los descendientes de todos los nodos anteriormente citados, debido a la convexidad de las familias de embudos (regla de la goma elástica, véase sección 2.2.7), no son visibles para ii por la presencia de u . El embudo posterior al último descendiente de uprima , se encontrará en el ala derecha y será nodor , y que será hijo o bien de nodou o bien de cualquiera de sus ascendientes en el árbol superior. Es por ello que se usará la siguiente estrategia: 1. Se parte de uprima y se busca si tiene algún hermano posterior a este en la familia de embudos. Si lo tiene, ese hermano se identifica como nodor y el nodo que precede al mismo se identifica como nodoq. 2. Si uprima no tiene hermanos CW , se busca si su ancestro directo tiene hermanos entre nodouprima e y . En caso afirmativo se para el bucle y se realiza la identificación ya mencionada de nodor ynodoq. 3. En caso contrario, se repite la búsqueda entre los antecesores hasta llegar a nodox , el embudo primigenio. Se adjunta para dos casos diferentes ejemplos del procedimiento, figura 2.26, en el caso A) se puede observar que nodouprima si tiene un hermano CW , por lo que la identificación de nodor , y consecuentemente la de nodoq es directa. Mientras que en el caso B) , es necesario remontar hacia el padre de upara encontrar un nodo válido. 2.3 Bucle en ii 53 x ii y x ii y u u' r u'' Antecesor en la familia Antecesor en la familia q ut u' r q Caso A): u' tiene hermano CW Caso B): Búsqueda en los antecesores de u Figura 2.26 Ejemplo de identificación nodor ynodoq. En el desarrollo teórico del documento original, [6], se plantea la deducción de estas propiedades a partir de lo que los autores denominan el reloj de arena,que consiste en la forma que presentarían los nodos de la familia si se representaran para una familia dada x−u−q−r−t−y , siempre y cuando nodoq y nodou sean embudos diferenciados. Tal y como se presenta la familia de embudos y sus propiedades, existe un cuello estrecho vacío de nodos conformado por nodou y nodot ; de forma adicional justifica que nodoq y nodor sean consecutivos en orden de nodo. Esto de forma teórica es muy interesante, pero en la práctica ciertas configuraciones pudieran no ser reducibles a estos parámetros, especialmente cuando existen pasillos estrechos después de haber procesado muchos polígonos, dándose varios relojes de arena. Se presentan los relojes de arena asociados a las familias propuestas para el cálculo de nodor en la figura 2.27. x ii y x ii y u r u'' q ut t u' u' r q Figura 2.27 Identificación de relojes de arena para los casos de 2.26. Nótese que al principio del código se considera el caso que el nodouprima no existe, por lo que el código genera dos ramificaciones que se vuelven a unir al final de la sección, con ambos valores de nodor ynodoq ya identificados según el caso. 933 %% Calculo q y r + Adición nodos solo visibles para (v,y) 934 %Se irá desde u hasta x viendo si hay hermanos 935 % ’Calculo q y r + Adición nodos solo visibles para (v,y)’ 936 nodoq=NaN; 937 nodor=NaN; 60 Capítulo 2. Método de Ghosh & Mount 1167 %Se busca el primer nodo de familiaxv_sub que se encuentre 1168 %en familia xv, -2 porque los últimos valores son nodoy e ii. 1169 1170 vecaux=familiaxv_sub(1:end-2,:); 1171 1172 %Este sigue el nodo propuesto como padre 1173 1174 [familiaxv]=add_familiaxv(EVG,familiaembudo,familiaxv,vecaux,nodox); 1175 1176 %El primer nodo es ii, el último es nodob, se añadirá a la salida del 1177 %bucle si es la última ejecución, o como nodoa en la siguiente 1178 vecaux=familiavy_sub(2:end-1,:); 1179 1180 [familiavy]=add_familiavy(EVG,familiaembudo,familiavy,vecaux,ii); 1181 1182 indiceaux=tkindrec(kk); 1183 nodoaux=familiaembudo(indiceaux,1); 1184 1185 EVG.Visibilidad(nodoaux,ii)=1; 1186 EVG.Visibilidad(ii,nodoaux)=1; 1187 1188 if kk==length(tkindrec) 1189 if familiaembudo(indiceaux,1)==familiaembudo(indiceaux,2) 1190 familiaxv=[familiaxv;familiaembudo(indiceaux,1),familiaxv(1,1)]; 1191 familiavy=[familiavy;familiaembudo(indiceaux,1),ii]; 1192 else 1193 [familiaxv]=add_familiaxv(EVG,familiaembudo,familiaxv,familiaembudo(indiceaux ,:),nodox); 1194 [familiavy]=add_familiavy(EVG,familiaembudo,familiavy,familiaembudo(indiceaux ,:),ii); 1195 end 1196 end 1197 %Fin bucle jj 1198 end 1199 else 1200 1201 indiceaux=tkindrec(1); 1202 nodoaux=familiaembudo(indiceaux,1); 1203 1204 EVG.Visibilidad(nodoaux,ii)=1; 1205 EVG.Visibilidad(ii,nodoaux)=1; 1206 1207 % El único punto que hay, es necesario añadirlo 1208 1209 if familiaembudo(indiceaux,1)==familiaembudo(indiceaux,2) 1210 familiaxv=[familiaxv;familiaembudo(indiceaux,1),familiaxv(1,1)]; 1211 familiavy=[familiavy;familiaembudo(indiceaux,1),ii]; 1212 else 1213 [familiaxv]=add_familiaxv(EVG,familiaembudo,familiaxv,familiaembudo(indiceaux,:), nodox); 1214 [familiavy]=add_familiavy(EVG,familiaembudo,familiavy,familiaembudo(indiceaux,:),ii); 1215 end 1216 end Código 2.33 Recurrencia tk:Parte2. Últimas comprobaciones Esta sección actúa como cajón de sastre, guardando ciertas acciones que se realizan en el código para garantizar ciertos estándares en la salida. La primera de estas acciones es comprobar que se estén conservando todas las familias que había en f amiliaembudo , en el caso de que la concurrencia no hubiera podido tener lugar. Para ello se revisan para las parejas de nodos de tkindrec , todos los nodos intermedios en vectorangulo , y se añaden con el padre que tenían asociado a f amiliaembudo. 1218 %% Se añaden nodos 0 entre t+1 e y 1219 2.3 Bucle en ii 61 1220 nodoscheckrecorrido=familiaembudo(tkindrec,1); 1221 1222 %Si por la razón que fuera no se hubieran añadido todos los nodos, se 1223 %buscan entre las parejas de tkindrec, los embudos que no estén y se 1224 %añaden. Se filtrarán y ordenarán fuera del bucle. 1225 1226 for jj=1:length(nodoscheckrecorrido)-1 1227 nodoseg=nodoscheckrecorrido(jj); 1228 nodosegpost=nodoscheckrecorrido(jj+1); 1229 1230 indnodoseg_vectorangulo=find(vectorangulo(:,1)==nodoseg); 1231 indnodosegpost_vectorangulo=find(vectorangulo(:,1)==nodosegpost); 1232 1233 vectorvalores_no1=vectorangulo((indnodoseg_vectorangulo+1):(indnodosegpost_vectorangulo -1),1); 1234 1235 if ~isempty(vectorvalores_no1) 1236 for kk=1:length(vectorvalores_no1) 1237 nodo1=vectorvalores_no1(kk); 1238 padrefamiliaembudo=familiaembudo(familiaembudo(:,1)==nodo1,2); 1239 vectorvalores_no1(kk,2)=padrefamiliaembudo; 1240 end 1241 indnodosegpost_famx=find(familiaxv(:,1)==nodosegpost,1); 1242 familiaxv=[familiaxv(1:(indnodosegpost_famx-1),:);vectorvalores_no1;familiaxv( indnodosegpost_famx:end,:)]; 1243 end 1244 end 1245 %Fin del if de familia resuelta 1246 end Código 2.34 Manteniendo consistencia de familias. Es necesario hacer notar que esto podría romper el orden de embudo de las familias, que antes podría haberse mantenido considerando la metodología de introducción, aunque con padres no actualizados, pero eso ya no sucede así. De todas maneras, no vale la pena realizar un esfuerzo adicional de implementación en esta parte del bucle debido a que una vez terminadas las instancias más externas de SPLIT se realizará un refactorizado completo, en el caso de que se hayan añadido. Desde el final de la última instancia de código presentada, se vuelven a considerar ahora aquellas ejecuciones en los que durante el cálculo de u , se computaron todos los nodos de la familia, es decir, pasó el end del if de resuelta f amilia == true . En la siguiente instancia de código se usará el comando unique en modo ′stable′ para eliminar embudos que pudieran estar duplicados, tanto para f amiliaxv como para f amiliavy . De forma adicional, se revisa si nodoy está en la última posición de f amiliavy dado que si no lo estuviera el algoritmo de reordenación de familias no funcionaría. 1248 %% Reasignación familiaxv,vy 1249 1250 [~,indiceunico]=unique(familiaxv(:,1),’stable’); 1251 familiaxv=familiaxv(indiceunico,:); 1252 1253 [~,indiceunico]=unique(familiavy(:,1),’stable’); 1254 familiavy=familiavy(indiceunico,:); 1255 1256 1257 indy=find(familiavy(:,1)==nodoy,1); 1258 if indy~=length(familiavy(:,1)) 1259 familiavy(indy,:)=[]; 1260 familiavy=[familiavy;nodoy,ii]; 1261 end Código 2.35 Limpieza de nodos duplicados. 62 Capítulo 2. Método de Ghosh & Mount Para cerrar la función es necesario añadir a f amiliaxv el nodo ii , que dado que la base de f amiliaxv es la conformada por x−ii , sería el último de la familia. Se añade con su padre siendo nodox y con esto concluye la ejecución de SPLIT. 1263 %% Completando familiaxv y familiavy 1264 %Hay que añadir ii a familia xv y ii se añadió a familiavy al principio de 1265 %la ejecución de SPLIT 1266 1267 familiaxv=[familiaxv;ii,nodox]; 1268 %FIN SPLIT 1269 end Código 2.36 Fin de SPLIT. 2.3.4 Procesado de las familias Una vez se han dividido las familias y acumulados todos los nodos visibles para las bases xv e vy , se recuerda que v≡ii , se realiza un proceso de almacenamiento en dos matrices llamadas TOT_f amxv y TOT_f amvy que guardan la información de las familias de par en par de columnas, una para los nodos hijos y otra para los padres. El número de pares de columnas es igual a las parejas de valores consecutivos que se encuentran en vectorrecorre , que conforman las diferentes bases de estas familias divididas, es decir, la longitud del mencionado vectorrecorre menos uno. Para poder introducir estas familias se comprueba si el número de filas de la susodicha es menor, igual o mayor a la que tiene la matriz. En caso de que sea menor, se amplia el número de filas utilizando zeros(l,2) y si es mayor se amplia el resto de la matriz con matrices de ceros para igualar. Estas familias se almacenan en una estructura de datos almacenada en el grafo de visibilidad extendido. 567 [EVG,familiaxv,familiavy]=SPLIT(nodox,nodoy,EVG,ii,loop,’interno’); 568 569 if jj==2 570 TOT_famxv=familiaxv; 571 TOT_famvy=familiavy; 572 else 573 if length(familiaxv(:,1))>=length(TOT_famxv(:,1)) 574 TOT_famxv=[[TOT_famxv;zeros(-length(TOT_famxv(:,1))+length(familiaxv(:,1)), length(TOT_famxv(1,:)))],familiaxv] ; 575 else 576 TOT_famxv=[TOT_famxv,[familiaxv;zeros(+length(TOT_famxv(:,1))-length( familiaxv(:,1)),2)]]; 577 end 578 579 if length(familiavy(:,1))>=length(TOT_famvy(:,1)) 580 TOT_famvy=[[TOT_famvy;zeros(-length(TOT_famvy(:,1))+length(familiavy(:,1)), length(TOT_famvy(1,:)))],familiavy] ; 581 else 582 TOT_famvy=[TOT_famvy,[familiavy;zeros(+length(TOT_famvy(:,1))-length( familiavy(:,1)),2)]]; 583 end 584 end 585 end 586 EVG.TOT_famxv_glob{ii}=TOT_famxv; 587 EVG.TOT_famvy_glob{ii}=TOT_famvy; Código 2.37 Almacenaje de los resultados obtenidos. En versiones anteriores del código sólo se trabajaba con las familias que resultaban de forma directa del SPLIT , lo que implicaba que solo las familias a las que se podían acceder era a las conformadas por xj j −vj j y vj j −yj j , siendo j j el índice que recorre las parejas de vectorrecorre . 2.3 Bucle en ii 63 En geometrías sencillas funcionaba bien, pero al aumentar el número de obstáculos y la complejidad de los mismos faltaba información por lo que el grafo quedaba incompleto. Procedimiento: Longitud de vectorrecorre es inferior a 2 Cuando el vectorrecorre solo tiene un valor, no se ha podido partir ninguna familia, por lo que se crean dos familias triviales, es decir, que solo contienen los correspondientes nodox y nodoy . Por otro lado se añade la línea de visión entre el valor presente de vectorrecorre , que se le llamará a por comodidad para representar a continuación la familia de embudos, e ii. Las familias que se añaden al diagrama de visibilidad extendida son las siguientes: f amiliaa,ii =a a ii a,f amiliavy =ii ii a ii(2.8) Cabe destacar que una de las familias apunta hacia dentro del polígono que comparten, y la otra hacia el lado externo. Pero ambas sólo contienen a su base, dado que las familias guardan información acerca de lo que se dirige hacia los valores de coordenadas xmenores. 588 elseif length(vectorrecorre)==1 589 590 EVG.TOT_famxv_glob{ii}=[vectorrecorre,vectorrecorre;ii,vectorrecorre]; 591 EVG.TOT_famvy_glob{ii}=[ii,ii;vectorrecorre,ii]; 592 593 EVG.Visibilidad(ii,vectorrecorre)=1; 594 EVG.Visibilidad(vectorrecorre,ii)=1; Código 2.38 Caso vectorrecorre contiene solo un valor. En el caso de que vectorrecorre esté vacío, situación que se da en recintos cóncavos, dado que en estos existe la posibilidad de que para tres puntos consecutivos del perímetro de un polígono, los extremos tengan mayor coordenada x que el de en medio. No es necesario añadir nada al grafo de visibilidad dado que no existe visibilidad ninguna de cara a valores con menor coordenada x , por lo que se añade un embudo trivial. 595 else 596 %Es necesario introducir esta información, esta debe ser 597 %visible para ser consultada. No conectar con ningún valor es un caso 598 %que puede darse, y se da, en recintos cóncavos 599 EVG.TOT_famvy_glob{ii}=[ii,ii]; 600 EVG.TOT_famxv_glob{ii}=[ii,ii]; 601 602 end Código 2.39 Caso vectorrecorre no contiene ningún valor. Adición de los nodos asociados a las familias al diagrama de visibilidad extendido Para poder computar las familias empleando add_f amilia2 , es necesario identificar todos los nodos asociados a las familias de tipo xv y de tipo vy obtenidas para ii. Para ello se emplean las estructuras de datos EVG.nodosTOT _f amxv_glob{ii}y EVG.nodosTOT_f amvy_glob{ii} , a las que se le añade el vector columna con los nodos únicos asociados a las familias, según si tienen orientación xv ovy. 604 %De esta manera se evita el cálculo repetido de estos valores 605 EVG.nodosTOT_famvy_glob{ii}=extraenodos(EVG.TOT_famvy_glob{ii}); 606 EVG.nodosTOT_famxv_glob{ii}=extraenodos(EVG.TOT_famxv_glob{ii}); 607 608 end Código 2.40 Extracción de los nodos de las familias asociadas a la base completa. 64 Capítulo 2. Método de Ghosh & Mount Por su parte extraenodos es una subrutina que agrupa todos los nodos hijos asociados a la ejecución, que se encuentran en las columnas impares de EV G.TOT_f amvy_glob{ii}y EVG.TOT_f amxv_glob{ii} , y se les aplica el comando unique , dado que al integrar diferentes familias podrían darse duplicados. 1function nodos_TOT=extraenodos(TOT_fam) 2if length(TOT_fam(1,:))==2 3%Esta información ya ha pasado por un unique 4nodos_TOT=[TOT_fam(:,1)]; 5else 6[~,d]=size(TOT_fam); 7nodos_TOT=[]; 8 9for jj=1:d/2%Siempre será par, por lo que no hay problema 10 subTOT=TOT_fam(:,2*jj-1); 11 indcero=find(subTOT==0,1); 12 13 if ~isempty(indcero) 14 subTOT=subTOT(1:indcero-1); 15 end 16 nodos_TOT=[nodos_TOT;subTOT]; 17 end 18 [~,indiceunico]=unique(nodos_TOT(:,1),’stable’); 19 nodos_TOT=nodos_TOT(indiceunico,:); 20 end 21 %Fin Bucle ii 22 end Código 2.41 Función extraenodos. Dependiendo de las preferencias del usuario, pueden comentarse o descomentarse las dos siguientes llamadas a figuras. Siendo la primera la que se emplea para pintar el diagrama al completo, y la segunda la que muestra de forma escalonada para cada nodo, como conecta con los anteriores, hasta llegar a N. 611 %% Postprocesado + Representación resultados RES 612 Descomentar si se quiere pintar la gráfica de el grafo de visibilidad al 613 completo 614 %%%%%%%%%%%%%%%%%% 615 fig=figure(2) 616 hold on 617 mapshow(V(:,1),V(:,2),’LineWidth’,1) 618 for i = 1:size(WPs, 1) 619 % Identificador basado en el orden 620 identificador = num2str(i); 621 % Mostrar el identificador junto al punto en el gráfico 622 text(WPs(i, 1), WPs(i, 2), identificador, ... 623 ’VerticalAlignment’, ’bottom’, ’HorizontalAlignment’, ’right’, ... 624 ’FontSize’, 6, ’Color’, ’r’); 625 end 626 627 % Personalizar el gráfico 628 xlim([min(V(:,1)),max(V(:,1))]) 629 ylim([min(V(:,2)),max(V(:,2))]) 630 ylabel(’Latitud’); 631 xlabel(’Longitud’); 632 title(’Diagrama de Visibilidad’); 633 grid on; 634 axis equal; 635 636 % for i = 5:N-4 637 for i = 1:N 638 % for j = i+1:N-4 639 for j = i+1:N 640 % for j = i+1:271 641 642 if EVG.Visibilidad(i, j) == 1 2.3 Bucle en ii 65 643 plot([WPs(i, 1), WPs(j, 1)], [WPs(i, 2), WPs(j, 2)], ’k--’); 644 end 645 end 646 647 end 648 649 set(fig, ’WindowState’, ’maximized’); 650 651 Descomentar si se quiere usar el modo manual para comprobar punto a punto 652 grafo de visibilidad. 653 654 %%%%%%%%%%%%%%%%%%%%% 655 figure(2) 656 hold on 657 mapshow(V(:,1),V(:,2),’LineWidth’,1) 658 659 for i = 1:size(WPs, 1) 660 % Identificador basado en el orden 661 identificador = num2str(i); 662 % Mostrar el identificador junto al punto en el gráfico 663 text(WPs(i, 1), WPs(i, 2), identificador, ... 664 ’VerticalAlignment’, ’bottom’, ’HorizontalAlignment’, ’right’, ... 665 ’FontSize’, 8, ’Color’, ’r’); 666 end 667 % Personalizar el gráfico 668 xlim([min(V(:,1)), max(V(:,1))]) 669 ylim([min(V(:,2)), max(V(:,2))]) 670 ylabel(’Latitud’); 671 xlabel(’Longitud’); 672 title(’Diagrama de Visibilidad’); 673 grid on; 674 axis equal; 675 676 % Dibujar las aristas de los polígonos primero 677 hPolygons = mapshow(V(:,1), V(:,2), ’LineWidth’, 1); 678 679 % Array para almacenar los handles de las líneas dibujadas 680 visible_lines = []; 681 invisible_lines = []; 682 683 for i = 1:N 684 for j = 1:i 685 if i ~= j % Evitar dibujar líneas de un punto a sí mismo 686 if EVG.Visibilidad(i, j) == 1 687 % Dibujar la línea visible y almacenar el handle 688 visible_lines(end+1) = plot([WPs(i, 1), WPs(j, 1)], [WPs(i, 2), WPs(j, 2)], ’k--’, ’LineWidth’, 1.5); 689 else 690 % Dibujar la línea no visible y almacenar el handle 691 invisible_lines(end+1) = plot([WPs(i, 1), WPs(j, 1)], [WPs(i, 2), WPs(j, 2)], ’r --’, ’LineWidth’, 1.5); 692 end 693 end 694 end 695 696 % Pausar 697 pause; 698 699 % Eliminar todas las líneas dibujadas en esta iteración 700 if ~isempty(visible_lines) 701 delete(visible_lines); 702 visible_lines = []; % Reiniciar el array de handles 703 end 704 705 if ~isempty(invisible_lines) 706 delete(invisible_lines); 707 invisible_lines = []; % Reiniciar el array de handles 708 end 709 end Código 2.42 Final Main_v13. 66 Capítulo 2. Método de Ghosh & Mount Con esto finaliza el bucle en ii , el grafo de visibilidad queda completado y solo queda extraerlo del diagrama de visibilidad extendido y pasarlo como salida de la función, junto con el tiempo computado tras el toc. 711 t2=toc; 712 % Guardar la visibilidad final 713 Visibilidad = EVG.Visibilidad; Código 2.43 Final Main_v13. Esta es la versión final implementada del código, la teoría que la sustenta y las diferencias presentes con el método originalmente planteado. En el siguiente capítulo se procederá a mostrar los resultados obtenidos con el mismo y el rendimiento que proporciona. 3 Estudio de resultados U na vez desarrollado el algoritmo implementado, se va a proceder a ponerlo a prueba con otros métodos y comprobar empíricamente la complejidad algorítmica real. Con este propósito se proponen los siguientes puntos de análisis: 1. Comparativa con el método de Lee: Para casos de tormentas reales, se procederá a trazar una comparativa para probar los resultados ofrecidos por la implementación. 2. Análisis de recintos complejos: Simulaciones en diferentes entornos repetibles, tiempo necesario para ello y prueba empírica de la complejidad algorítmica del método. 3. Profile analysis: Se comprobará para un caso concreto de ejecución en que operaciones se ha invertido más tiempo, propuestas de desarrollo y limitaciones. 3.1 Comprobación del caso de estudio: Método de Lee en tormenta realista Debido a las diferencias en el planteamiento de ambos problemas, el que resuelve SPLIT y el método de Lee, resulta delicado medir ambos programas en igualdad de condiciones. Por un lado el método de Lee trabaja con un destino y un origen en un entorno abierto, mientras que el método SPLIT necesita un recinto cerrado pero no tiene la necesidad de plantear un origen y un destino, pero por otro lado no permite que dos puntos compartan coordenada x y no permite la existencia de 3 puntos colineares. Esto se ha tenido en mente desde el principio del proyecto, procurando hacer que ambos algoritmos puedan trabajar, si bien no con los mismos inputs, con cierto preprocesado. Gracias a esto pueden trazarse comparativas aproximadas entre ambos métodos. Para ello se tomará el caso de estudio utilizado en el análisis de resultados desarrollado por Narciso Valverde, [5], se explicarán qué cambios se han tenido que realizar para poder realizar esta adaptación y se mostrarán los resultados obtenidos, comprobando que el número de conexiones del grafo es igual para ambos casos. Para poder ejecutar un mapa de coordenadas preparado para ser usado con el método de Lee se plantea la siguiente metodología: 1. Se verifica que los polígonos estén dados en sentido contrario a las agujas del reloj en torno a su centroide, en caso contrario se reordenan los nodos asociados a los polígonos en la matriz de coordenadas. 2. Se identifican los valores menores y mayores para las coordenadas x e y . Si hay algún punto fuera del cuadrante x positivo e y positivo se desplazan todos los puntos y después se aplica un escalado. 67 68 Capítulo 3. Estudio de resultados 3. Finalmente, tras esta normalización, se añaden dos triángulos auxiliares: uno con un vértice en el punto de origen y otro con un vértice en el punto de destino. En este aspecto, el método SPLIT es más inflexible. El caso contrario, pasar de mapa preparado para ser usado con el método SPLIT a método de Lee, es bastante más sencillo, consistiendo en eliminar el recinto exterior e introducir dos coordenadas, origen y destino, al final de la matriz de coordenadas. Se muestran a continuación los resultados obtenidos de ejecutar el mapa de estudio del trabajo de Narciso Valverde, [5], con el método SPLIT, el grafo de visibilidad obtenido es el mismo una vez se descartan las lineas de visión asociadas al recinto exterior y a los puntos que no son de interés de los triángulos auxiliares de destino y origen. Figura 3.1 Ejemplo de región delimitada por meteorología adversa. El tiempo de ejecución de esta simulación es de t=210,57s en el modo de ejecución de SPLIT adaptado a este tipo de geometrías, presentando un rendimiento menor al del método de Lee, véase la tabla 3.1. Por otro lado el número de puntos de paso de la simulación asciende a N=1206 , 8 puntos más que para el método de Lee debido a la necesidad técnica de añadir el recinto exterior, 4 puntos, y los nodos auxiliares no utilizados de los triángulos, otros 4 puntos. Cabe destacar que la metodología empleada en el trabajo de Narciso Valverde [5] fue diferente: partiendo de un caso específico, se realizaron varias simulaciones para calcular el tiempo promedio de ejecución en distintas versiones reducidas del mismo caso, utilizando el algoritmo de RamerDouglas-Peucker. Este algoritmo permite reducir el número de puntos de los polígonos de un recinto dado, manteniendo aproximadamente sus geometrías para las tolerancias establecidas. Siguiendo esta metodología, Narciso Valverde, en su estudio [5], aplicó diferentes reducciones en el número 3.2 Ejecución en casos propuestos: Recintos convexos 69 de nodos del recinto mostrado en esta sección, registrando los resultados correspondientes para cada valor de Ngenerado por la simplificación. Comparando estos valores con las ejecuciones presentes en [5], mostrados en la tabla 3.1, los tiempos son superiores a los obtenidos mediante el método de Lee pero inferiores a los obtenidos mediante el método Ingenuo avanzado que se planteo en su investigación. Estos, además, se encuentran en orden de magnitud más cercanos a los obtenidos mediante el método Ingenuo, debido a las limitaciones inherentes de la implementación. Tabla 3.1 Comparación de tiempos de ejecución: Algoritmo basado en el método Ingenuo, método SPLIT y el método de Lee. Obtenida de [5]. N (-) Promedio de tIng,v2(s) Promedio de tSPLIT (s) Promedio de tLee (s) 1198 454.5 210.57 22.7 Se puede observar en la imagen 3.1 que SPLIT respeta de forma adecuada las conexiones del grafo de visibilidad en los puntos ubicados en las concavidades de los obstáculos, lo cual se considera un objetivo cumplido. Además las envolventes se encuentran correctamente computadas y las zonas de alta densidad de grafos corresponden con las zonas con mayor densidad de puntos expuestos a otros obstáculos. Más allá de esta comprobación cualitativa, se han comparado los grafos de visibilidad propuestos por los método Ingenuo y de Lee, arrojando el mismo número de conexiones, validando así el modelo. 3.2 Ejecución en casos propuestos: Recintos convexos Durante el transcurso de la implementación se han establecido procedimientos y herramientas para poder depurar la algoritmia, todos estos estarán en sus respectivas secciones del Apéndice A. Entre las herramientas más útiles generadas, de cara a encontrar fallos de forma eficaz, es lo que se ha llamado un generador de recintos determinista. Para un mismo número de puntos, la ejecución puede variar mucho dependiendo de diversos factores, como la posición relativa entre los obstáculos y el tamaño de los mismos. Es por ello que, buscando estas situaciones límite, se han creado dos generadores de recintos, uno de entornos complejos y otro de entornos sencillos, con los que se ha ejecutado el algoritmo en busca de fallos para poder depurarlos. Se ha tomado la decisión de emplear recintos convexos y no cóncavos debido a que, para un mismo número de puntos de paso, los recintos cóncavos requieren un menor esfuerzo computacional. Esto es causado por la naturaleza dependiente de las conexiones en el grafo. Dado que se conecta un punto con otros que sean visibles, los recintos cóncavos tienen regiones en las que sólo pueden ver su entorno, lo cual reduce significativamente el esfuerzo respecto a un caso convexo, en el que la mayor parte de los nodos se encuentran expuestos a muchos otros polígonos, siempre y cuando no haya obstáculos muy voluminosos. La matemática detrás del generador de recintos es sencilla, para poder generar entornos repetibles y comprobar de manera sencilla que las soluciones programadas tienen un impacto positivo en la implementación, se entrega una semilla seed para que todo lo generado aleatoriamente, usando los comandos rand , sea repetible consistentemente. Para poder definir la complejidad del recinto, se ha tomado la solución de indicar el número de polígonos que se necesitan al programa. De esta manera cada polígono se añade sobre una base ya computada, lo que evita que un entorno con una semilla concreta no varíe sustancialmente al querer añadir un nuevo punto. 76 Capítulo 3. Estudio de resultados Figura 3.6 Estudio del comportamiento del número de conexiones en el grafo, E , respecto al número de puntos de paso N. 3.3.3 Desglose de tiempos utilizando profiler Para ver como de eficiente es el algoritmo internamente, más allá de la complejidad algorítmica mostrada, se va a hacer uso de la herramienta de profiler. Esta permite para un intervalo de código concreto, al igual que se haría con un tic-toc, comprobar cuanto tiempo se ha invertido en cada sección del código pudiendo así optimizar de manera más eficiente. Esta herramienta ha sido ampliamente usada a lo largo del desarrollo. Dado que se estaba realizando una implementación no fiel a la estructura de datos originalmente propuesta, para conseguir resultados aceptables en cuanto a tiempos de ejecución, se ha estado vigilando las partes internas del programa. Buscando diversas estrategias de optimización minimizando el número de cálculos y maximizando el número de casos que pueden ser procesados con garantía de validez de resultados. 1profile clear 2profile on 3 4[t2,visibilidad]=Main_v13(V); 5 6profile off 7profile viewer Código 3.2 Ejemplo de uso de profiler. Los argumentos utilizados para sacar el máximo partido a la herramienta son los siguientes: 1. clear: Limpia el historial registrado. 2. on: Indica a Matlab que empiece a registrar los datos de la ejecución. 3. off : Análogo a on pero para que se detenga. 4. viewer: Abre el lector, que contiene la información de interés. 3.3 Resultados de la ejecución 77 Cabe destacar que el empleo de la misma empeora los tiempos de rendimiento, dado que se registran de forma activa los tiempos de ejecución y llamadas a las funciones. Por suerte, esto no resulta un gran problema, dado que es una herramienta de depuración. Para un caso de estudio como el de los recintos convexos planteado anteriormente, seed =1 y npol =100, se pueden observar a continuación los tiempos de computación obtenidos: Figura 3.7 Resultados obtenidos del profiler para seed =1ynpol =100. Los dos principales sumideros de carga computacional durante la ejecución son, precisamente, los derivados de las características de la implementación, inter_visible1 que calcula si existen algún corte entre las aristas vinculadas a los nodos visibles para la base de la familia de embudos que se esté SPLITeando, con la línea de visión entre ii y el nodo en cuestión. corta_puntos es una función hija de inter_visible que calcula si dos segmentos de recta intersecan entre sí. Por otro lado add_f amilia2 que es la función que calcula la familia de embudos que se solicite ordenando los nodos asociados a la base de la familia. add_f amilia e inter_visible1 ocupan aproximadamente un 90% del tiempo total de ejecución. La complejidad algorítmica propuesta por Ghosh & Mount resultaría factible entonces, de implementarse con los requisitos mencionados, pudiendo mejorarse así la complejidad algorítmica propuesta por el método de Lee, O(N2logN), [5]. Ambas ejecuciones son críticas y afectan en gran manera al resultado final obtenido, por ejemplo, si la familia obtenida está mal ordenada porque existiera algún nodo que se ha insertado en el padre equivocado, todo el fundamento teórico en el que se basa SPLIT queda completamente invalidado, pudiendo provocarse falsos positivos y falsos negativos, que además irían propagándose durante el resto de la ejecución completa del bucle debido a a naturaleza aditiva del método. Continuando esta línea de pensamiento, si debido a que la información de las familias no se encuentra completa a causa de un problema con la ordenación eficaz de los nodos, podrían darse falsos positivos al no encontrarse aristas que produzcan un corte. Aunque es posible, que debido a 78 Capítulo 3. Estudio de resultados la robustez artificial proporcionada al método, con varias capas de comprobaciones y condiciones, estos fallos no lleguen a propagarse, pero la complejidad computacional se ha visto lastrada por esto. Una de estas comprobaciones es, por ejemplo, el añadir a la lista de nodos de la familia de embudos, los vértices de los polígonos que están en la familia. Esto es necesario para ejemplos como el de la figura 3.8. En polígonos altamente complejos las estrategias de visibilidad necesitan ser adaptadas, es por ello que existe una versión a y b del código, cada una implementando estrategias de visibilidad distintas con este propósito. Todas estas técnicas, al final, se sirven como refuerzo al concepto de los embudos y sus familias. En última instancia ha sido necesario establecer un equilibrio entre garantizar familias completas y realizar comprobaciones de visibilidad mediante cortes. ii x y Figura 3.8 Caso de intersección no detectada sin comprobaciones adicionales de todos los vértices asociados a los polígonos de la familia de embudos. El método de cálculo de visibilidad presentado en la sección 2.3.3, la basada en la definición de los embudos, llegó a funcionar de forma sobresaliente en etapas tempranas de la implementación, pero acabó siendo dejada en segundo plano en aras de garantizar en mejor medida resultados fiables, debido a la casuística comentada en la sección mencionada anteriormente. El hecho de que la escalabilidad de la complejidad algorítmica calculada empíricamente presente resultados tan notables, aunque los tiempos de computación sean elevados, no es si no un indicador del potencial que tiene este método. Aunque se realizará una disertación en las conclusiones acerca de si la relación costo-riesgo-beneficio de trabajar con algoritmos así de complejos pudiera resultar a la larga una desventaja. Finalmente, los resultados obtenidos son mejores que los del algoritmo Ingenuo debido a que aunque se pueda llegar a hacer un número de comprobaciones considerables de cortes de segmentos, se hace siguiendo la estrategia que proporcionan las familias. Esto representa en sí, un éxito. 4 Conclusiones y futuras líneas de trabajo P resentado el método, la implementación y los resultados obtenidos se desarrollarán las lecciones aprendidas, diversas puntualizaciones que pudieran resultar de interés para el lector y diferentes propuestas de desarrollo. El primer punto a tratar será la complejidad de la implementación del algoritmo, ¿es eficiente desde un punto de vista ingenieril considerando el entorno que rodea a la resolución del problema? El principal atractivo de implementar un algoritmo así de complejo, pero que escala de manera tan eficiente, es poder computar de forma muy precisa rutas que tanto directamente como indirectamente, reportan un beneficio a diferentes escalas, tal y como se planteó en la introducción. Fuera del entorno teórico existen varias limitaciones y restricciones que podrían echar por tierra estos beneficios obtenidos. ¿Es necesario realmente considerar cada obstáculo por separado y respetando sus concavidades? Podría plantearse el caso de que considerar estas concavidades es imprescindible si, por ejemplo, se quisiera aterrizar de forma segura en un aeropuerto que se encontrara ubicada en una; pero el análisis realizado se ha hecho para un único nivel de vuelo, y lo que es más crítico, considerando los obstáculos fijos. Por supuesto, existen en la vida real obstáculos fijos, como áreas de vuelo restringidas delimitadas geográficamente, pero las rutas o sugerencias de ruta de un vuelo comercial están trazadas de forma que estas condiciones de contorno estén evitadas por defecto, dado que todos los actores que toman parte del mundo aeroespacial deben ser consciente de estas limitaciones a la hora de establecer rutas de acuerdo con los estándares de seguridad operacional y física esperadas en el sector. Es por ello, que a la hora de la verdad, los agentes meteorológicos juegan un papel muy importante a la hora de tomar decisiones en el día a día. No se plantean de la misma manera una operación de aterrizaje o despegue con condiciones óptimas de visibilidad y/o viento alineado con la pista, que días con niebla, nieve, lluvia intensa o fuertes vientos cruzados. En condiciones extremas, se encuentran contemplados procedimientos para aterrizar en aeropuertos cercanos e incluso que la torre de control no proporcione autorización para despegar. Las tormentas se desarrollan con el tiempo en la atmósfera, pudiendo cambiar su posición y forma en ordenes de magnitud temporales inferiores a la duración estándar de un vuelo comercial. Entonces cabe plantearse lo siguiente, ¿resulta interesante usar algoritmos muy complejos que puedan detectar caminos y pasadizos angostos en un contexto aeronáutico? ¿Existen garantías de que esa ruta exista una vez comenzada? Asimismo, el cálculo desde tierra de tormentas para un nivel de vuelo dado no es tan preciso como el algoritmo que se ha planteado, [11]. Todos estos factores, acumulados, desde el punto de vista del redactor, hacen más interesante de cara a líneas de desarrollo futuras, optar por estrategias que trabajen con geometrías más sencillas y entornos dinámicos, estableciendo metodologías de decisión de ruta basados en el control en bucle cerrado con diferentes horizontes temporales. 79 80 Capítulo 4. Conclusiones y futuras líneas de trabajo Por otro lado, la implementación, aunque haya satisfecho las expectativas propuestas, plantea una serie de dudas. Los resultados son satisfactorios en cuanto a la capacidad de escalado del método con N , pero la implementación es tan robusta como sensible. Están programados dentro del código secciones enteras de comprobaciones adicionales que han tenido que desestimarse en aras de mantener la complejidad en rangos aceptables, lo cual no ha representado un problema en la extensa variedad de casos de estudio depurados pero cabe la posibilidad de que se dieran inconsistencias, si se considera la casuística en geometrías increíblemente complejas. Continuando con esta línea de desarrollo, establecer la complejidad algorítmica en base al número de conexiones visibles podría considerarse delicado, dado que esa información no se conoce de antemano. Es por ello que el artículo original se traduce como Un algoritmo sensible a los resultados para resolver grafos de visibilidad. Todo esto deriva a que los tiempos de computación no sean consistentes entre sí, dado que diferentes geometrías pueden llevar a comportamientos del algoritmo diametralmente opuestos. En el capítulo anterior se pudo observar, que para un número similar de puntos, existen casos que se resuelven en la mitad de tiempo debido a sus peculiaridades. Después de analizar los resultados, se ha estimado que E , el número de conexiones en un grafo de visibilidad en el contexto de los casos de estudio es aproximadamente: E≈N√2,(4.1) que es sin lugar a dudas, superior a la complejidad algorítmica del término que si era dependiente de N, O(NlogN+E≈N√2);(4.2) que es similar a la pendiente de la curva de regresión potencial propuesta en el apartado anterior. Aún así, en el propio documento original no se especifica directamente la complejidad algorítmica del algoritmo de triangulación empleado, sólo la complejidad del de Mehlhorn, [9]. Este algoritmo no está planteado para funcionar en recintos con agujeros, pero se indica que este se generaliza fácilmente, [6], sin indicar ninguna estrategia, ni aportar consideraciones adicionales acerca de como esto afecta a la complejidad algorítmica del método de triangulación original. De la misma manera, no se han especificado condiciones de parada en la recurrencia propuesta, pero la implementación realizada consigue sacar partido de ella, siendo esta consistente con los tiempos de computación. Cabe recalcar que, existe una variabilidad de casos que no llega a ser completamente cubierta en el documento original, habiéndose considerado imprescindible aplicar varias capas de parches para cubrir los casos límite, véase la referencia [12], dónde se comentan ciertas dificultadas similares a las documentadas en esta memoria. Estas situaciones, originalmente planteadas como límite y no cubiertas del todo por la casuística, en entornos muy complejos y con muchos obstáculos, podrían darse de forma recurrente. Otra desventaja a comentar sería la necesidad técnica de realizar un preprocesado más extenso que para otros métodos, debido a que se pide como dato, todos los nodos ordenados unos alrededor de otros, tal y como se especifica en la sección 7 de [6], o que los puntos de paso deban estar ordenados en orden de x creciente, esto introduce cuellos de botella independientemente de la estructura de datos que se emplee. En última instancia, de considerarse apropiado continuar esta línea de investigación a la hora de resolver el grafo de visibilidad, se recomendaría implementar la estructura de datos según lo planteado en el documento original, la propuesta por Gabow & Tarjan, [10]. Las decisiones adoptadas en este trabajo se han visto influenciadas por consideraciones relacionadas con el lenguaje de programación utilizado, Matlab™, los conocimientos específicos en el ámbito de la ciencia de la computación y el horizonte temporal del mismo. Aún así, bajo estas circunstancias se han podido alcanzar de forma satisfactoria los objetivos originalmente propuestos. Apéndice A Códigos de Matlab En este apéndice se recopilarán todos los códigos que no hayan sido expuestos a lo largo del documento. Se hará un esquema y resolverán posibles dudas que pudieran surgir acerca de las relaciones de los mismos. A.1 Relaciones entre los códigos A lo largo de los meses invertidos en este trabajo, se ha generado una gran cantidad de código que recoge varios programas y subrutinas que sirven de apoyo a las secciones principales comentadas a lo largo del presente documento. Existen 5 archivos con extensión .m que recogen la última versión de los algoritmos empleados: 1. Programas parent: Sirven de hub y gestionan las llamadas a las funciones child. En estas se realiza gestiona el entorno que envolverá a la ejecución. • Ejecutaprograma: Diseñado para las ejecuciones únicas del programa: Gestiona la creación de un recinto, llama a Main_v13_C y al profiler. • Estudio_parametrico: Usando el módulo de Parallel Computing, realiza un barrido para tantas semillas como se le indique y para valores de polígonos generados por un logspace. Toda esta información se almacena en un archivo de extensión .xlsx, para procesar los datos de forma rápida con las herramientas manuales que proporciona. • Compatibiliza: Importa los datos empleados en [5]. Adapta la matriz de coordenadas, debido a que los nodos en los polígonos están dados en sentido horario respecto al centroide y llama a la función Main_v13_C para obtener el grafo de visibilidad. 2. Programas children: su ejecución está autocontenida y proporcionan una salida concreta. •calcula_recinto: Devuelve un recinto parametrizado en función de seed ynpol. •Main_v13_C: Devuelve el grafo de Visibilidad y el tiempo necesario para ejecutarlo. Se proporcionarán los códigos en el orden en el que se han descrito y todas las subfunciones que pudieran tener serán añadidas con una breve contextualización que desarrolle los comentarios adjuntos en el código. A.2 Programas parent A.2.1 Ejecutaprograma 81 82 Capítulo A. Códigos de Matlab 1%% Código de arranque del programa 2clc, clear all 3 4seed=2; 5npol=50; 6 7V=calcula_recinto(seed,npol); 8 9profile clear 10 profile on 11 12 [t2,visibilidad]=Main_v13_C(V); 13 t2 14 15 profile off 16 profile viewer Código A.1 Ejecutaprograma.m. A.2.2 Estudio_parametrico 1%% Itera_ejecuciones en paralelo 2tic 3reiniciar = true; % Cambia a ’false’ si no se desea borrar el contenido 4 5% Si se indica reiniciar, borra todo el contenido del archivo 6if reiniciar 7delete(’datos_ejecucion_.xlsx’); 8end 9 10 % Inicializa un grupo de trabajadores si no está activo 11 if isempty(gcp(’nocreate’)) 12 parpool; % Inicia un pool de procesamiento en paralelo 13 end 14 15 % Número de iteraciones para cada bucle 16 nseed = 50; 17 18 % Utiliza parfor para el bucle exterior en paralelo 19 parfor seed = 1:nseed 20 for npol = round(logspace(1,2.5,10)) 21 22 V = calcula_recinto(seed, npol); 23 [t2, visibilidad] = Main_v13(V); 24 25 % Guardar los resultados en el archivo Excel 26 guarda_excel(t2, visibilidad, seed, npol); 27 end 28 strcat(num2str(seed)) 29 end 30 31 t=toc 32 33 function guarda_excel(t2,visibilidad,seed,npol) 34 35 % Nombre del archivo Excel 36 archivo = ’datos_ejecucion.xlsx’; 37 38 % Definir los valores que se quieren registrar 39 40 visibilidad(isnan(visibilidad))=0; 41 E=sum(sum(visibilidad)); 42 [N,~]=size(visibilidad); 43 version=13; 44 45 % Cargar los datos previos, si existen, para conocer la próxima fila vacía 46 if isfile(archivo) A.2 Programas parent 83 47 % Leer el contenido del archivo para determinar la próxima fila vacía 48 datosPrevios = readcell(archivo); 49 filaInicio = size(datosPrevios, 1) + 1; % Fila vacía 50 else 51 % Si el archivo no existe, se creará y se iniciará en la primera fila 52 filaInicio = 1; 53 end 54 55 % Escribe el encabezado si es la primera vez 56 if filaInicio == 1 57 encabezado = {’seed’, ’npol’, ’N’, ’t2’,’E’,’version’}; 58 writecell(encabezado, archivo, ’Sheet’, 1, ’Range’, ’A1’); 59 filaInicio = filaInicio + 1; % La siguiente fila vacía 60 end 61 62 % Escribe los datos en la fila disponible 63 nuevaFila = {seed, npol,N, t2,E,version}; 64 writecell(nuevaFila, archivo, ’Sheet’, 1, ’Range’, [’A’ num2str(filaInicio)]); 65 66 end Código A.2 Estudio_parametrico. A.2.3 Compatibiliza 1clc,clear all, close all 2% Compatibiliza la matriz de coordenadas para el caso de estudio del método 3% de Lee 4load(’TFG_Narciso.mat’) 5 6V=[LON,LAT]; 7%% Transformación de V 8% En este vector de coordenadas los datos de los nodos asociados a los 9% polígonos están dado en orden antihorario respecto al centroide, por lo 10 % que se realizará a continuación una adaptación para que sea compatible 11 % con el código Main_v13 12 13 poligonos = {}; % Celdas para almacenar los polígonos 14 current_poly = []; % Almacena los nodos actuales 15 16 % Separar los polígonos utilizando NaN 17 for i = 1:size(V, 1) 18 if any(isnan(V(i, :))) % Si es un NaN 19 if ~isempty(current_poly) % Si hay nodos acumulados 20 poligonos{end + 1} = current_poly; % Guardar el polígono 21 current_poly = []; % Reiniciar 22 end 23 else 24 current_poly = [current_poly; V(i, :)]; % Añadir nodo 25 end 26 end 27 28 % Agregar el último polígono si existe 29 if ~isempty(current_poly) 30 poligonos{end + 1} = current_poly; 31 end 32 33 % Reordenar los nodos en cada polígono 34 for i = 1:length(poligonos) 35 if ~isempty(poligonos{i}) 36 % Cambiar el orden; aquí lo hacemos de la manera deseada 37 poligonos{i} = poligonos{i}(end:-1:1, :); % Invertir el orden 38 end 39 end 40 41 % Combinar los polígonos reordenados en un solo vector 42 V_reordenado = []; 43 for i = 1:length(poligonos) 84 Capítulo A. Códigos de Matlab 44 V_reordenado = [V_reordenado; poligonos{i}]; % Añadir el polígono 45 V_reordenado = [V_reordenado; NaN(1, 2)]; % Añadir el NaN separador 46 end 47 48 % Eliminar el último NaN si es necesario 49 if all(isnan(V_reordenado(end, :))) 50 V_reordenado(end, :) = []; 51 end 52 53 % Mostrar el vector reordenado 54 disp(V_reordenado); 55 V=V_reordenado; 56 57 V400=V; 58 59 vminx=min(V400(:,1)); 60 vminy=min(V400(:,2)); 61 62 border=10; 63 64 Vmod=zeros(size(V400)); 65 66 if vminx<=0 67 Vmod(:,1)=V400(:,1)+vminx+border; 68 Vmod(:,2)=V400(:,2); 69 % vminx=border; 70 71 elseif vminy<=0 72 Vmod(:,2)=V400(:,2)+vminy+border; 73 % vminx=border; 74 Vmod(:,1)=V400(:,1); 75 76 else 77 Vmod=V400; 78 end 79 80 vminx=min(Vmod(:,1))*0.9; 81 vminy=min(Vmod(:,2))*0.9; 82 83 vmaxx=max(Vmod(:,1))*1.1; 84 vmaxy=max(Vmod(:,2))*1.1; 85 86 Vmod(:,1)=(Vmod(:,1)-vminx)./((vmaxx-vminx)); 87 Vmod(:,2)=(Vmod(:,2)-vminy)./((vmaxy-vminy)); 88 89 Vmod=[0,0;1,0;1-1e-4,1;1e-4,1,;0,0;NaN(1,2);Vmod;NaN(1,2)]; 90 91 indprimernan=find(isnan(Vmod(:,1)),1); 92 93 vminpostx=min(Vmod(indprimernan+1:end,1)); 94 vmaxpostx=min(Vmod(indprimernan+1:end,1)); 95 96 D1=vminpostx-Vmod(4,1); 97 D2=abs(vmaxpostx-Vmod(3,1)); 98 99 origen=[Vmod(4,1)+0.25*(D1),0.5]; 100 101 destino=[Vmod(3,1)-0.1*(D2),0.5]; 102 103 poligonoorigen=[origen;Vmod(4,1)+0.1*(D1),0.5*1.1;Vmod(4,1)+0.1*(D1),0.5*0.9;origen]; 104 poligonodestino=[destino;Vmod(3,1)-0.01*(D2),0.5*0.9;Vmod(3,1)-0.01*(D2),0.5*1.1;destino]; 105 106 Vmod=[Vmod;poligonoorigen;NaN(1,2);poligonodestino]; 107 108 [t2,visibilidad]=Main_v13_C(Vmod); 109 t2 Código A.3 Compatibiliza. A.3 Programas Children 85 A.3 Programas Children A.3.1 calcula_recinto Esta función realiza la algoritmia desarrollada en 3.2, genera recintos basados en una semilla dada, configurada mediante rng(seed)y normalizados para xeyentre [0,100]. Modificando los siguientes parámetros puede cambiarse la naturaleza de los recintos: •maxaristaspol: Número máximo de aristas permitidos por polígono. •minrad : Tamaño mínimo del radio en relación a la posición frente a la arista del recinto externo más cercano. •f actormargen : Define la mínima distancia a mantener entre el centroide del polígono y el margen. Esta función tiene tres funciones children asociadas: • distancia: Calcula la distancia entre dos puntos dado su vector director. Su usa para comprobar que la distancia entre centroides sea suficiente para evitar interferencias. • distance_point_to_polygon: Calcula la mínima distancia del centroide a añadir respecto al recinto exterior. •point_to_segment_distance: Calcula la mínima distancia entre un punto y una recta. 1function V=calcula_recinto(seed,npol) 2%% Inicio Programa 3 4%Se va a generar un código que reciba como input el número de polígonos, 5%máximo de lados por cada uno, y diferentes parámetros de forma y tamaño y 6%calcule los diferentes dominios y los de como input al algoritmo. 7 8 9%% Marco inicial 10 V0=[0,0;100,0;99.99,100;0.01,100;0,0;NaN(1,2)]; 11 V=[0,0;100,0;99.99,100;0.01,100;0,0;NaN(1,2)]; 12 13 rng(seed) 14 maxaristaspol=6; 15 factormargen=10;%10 16 17 minrad=20;%20 %Tamaño mínimo en relación al margen 18 19 margenx=(max(V(:,1))-min(V(:,1)))/factormargen; 20 margeny=(max(V(:,2))-min(V(:,2)))/factormargen; 21 22 margen=max([margenx,margeny]); 23 vectorcentro=[]; 24 vectordist=[]; 25 26 for ii=1:npol 27 28 xcentro=margen+(100-2*margen)*rand(1); 29 ycentro=margen+(100-2*margen)*rand(1); 30 R=min([distance_point_to_polygon(V0, xcentro, ycentro)*rand(1),minrad]); 31 contador=1; 32 33 if ~isempty(vectorcentro) 34 while contador<=length(vectorcentro(:,1)) 35 36 if dist([xcentro,ycentro],vectorcentro(contador,:))<(R+vectordist(contador))*1.1 37 xcentro=margen+(100-2*margen)*rand(1); 38 ycentro=margen+(100-2*margen)*rand(1); 92 Capítulo A. Códigos de Matlab 106 107 108 vectoraristas=[seg_aristas(indunos:end,:);seg_aristas(1:indunos-1,:)]; 109 110 indunos=find(vectoraristas(:,2)==1); 111 112 p1=EVG.WPs(vectoraristas(indunos(1),1),:)’; 113 p2=EVG.WPs(vectoraristas(indunos(2),1),:)’; 114 115 v1=p1-p4; 116 v2=p2-p4; 117 118 prodvec1=v1(1)*v2(2)-v1(2)*v2(1); 119 120 %por definición, uno va a quedar más CW que el otro, buscaremos cual es 121 %cual y a partir de ahí se buscan los valores factibles para nodos padres 122 123 if prodvec1>0 124 %p2 está por encima de p1 (lead de indunos(1)) 125 padresposiblespol=vectoraristas(indunos(1):(indunos(2)-1),1); 126 padresnoposiblespol=vectoraristas(indunos(2):end,1); 127 128 for jj=1:length(padresposiblespol) 129 nodopos=padresposiblespol(jj); 130 EVG.MATposiblepadre(nodox,nodopos)=1; 131 end 132 133 for jj=1:length(padresnoposiblespol) 134 nodonopos=padresnoposiblespol(jj); 135 EVG.MATposiblepadre(nodox,nodonopos)=0; 136 end 137 else 138 %p1 está por encima de p2 (lead de indunos(2)) 139 padresposiblespol=vectoraristas(indunos(2):end,1); 140 padresnoposiblespol=vectoraristas(indunos(1):(indunos(2)-1),1); 141 142 for jj=1:length(padresposiblespol) 143 nodopos=padresposiblespol(jj); 144 EVG.MATposiblepadre(nodox,nodopos)=1; 145 end 146 147 for jj=1:length(padresnoposiblespol) 148 nodonopos=padresnoposiblespol(jj); 149 EVG.MATposiblepadre(nodox,nodonopos)=0; 150 end 151 end 152 end 153 end Código A.6 Función test_de_paternidad_MAIN. cortapuntos Una forma más de calcular intersecciones, esta vez resolviendo un sistema de ecuaciones con el método de Cramer empleando las formulación pendiente-punto de las rectas. Este método presenta el inconveniente de que calcula el punto de intersección independientemente de si se encuentra dentro o fuera de los segmentos de las rectas, por lo que se verifica que el punto de intersección se encuentra entre ambas con una sentencia if. 1%% Función Comprueba visibilidad aristas (Intersección) 2 3function [inter,xinter]=cortapuntos(line1,line2) 4%Las lineas vendrán dadas en el formato 5% [x1, x2 6% ;y1,y2] 7 8%line1 es la línea que une N con el punto de la secuencia a comprobar A.3 Programas Children 93 9%line2 es la línea con la que se busca que interseque (cualquiera de las aristas de los polí gonosconsiderados) 10 11 x1=line1(1,1); 12 y1=line1(2,1); 13 x2=line1(1,2); 14 y2=line1(2,2); 15 %Por simplicidad vamos a forzar 16 17 if x1>x2 18 xaux=x2; 19 yaux=y2; 20 21 x2=x1; 22 y2=y1; 23 24 x1=xaux; 25 y1=yaux; 26 27 end 28 29 x3=line2(1,1); 30 y3=line2(2,1); 31 x4=line2(1,2); 32 y4=line2(2,2); 33 34 m1=(y2-y1)/(x2-x1); 35 36 37 if x3>x4 38 xaux=x4;yaux=y4; 39 x4=x3;y4=y3; 40 41 x3=xaux; 42 y3=yaux; 43 44 end 45 46 m2=(y4-y3)/(x4-x3); 47 48 if m1~=m2 49 xinter=(y3-y1+m1*x1-m2*x3)/(m1-m2); 50 51 %Esta intersección podría provocarse tanto antes como después del 52 %segmento que nos interesa comprobar su visibilidad por lo que habría 53 %que comprobar si la intersección corta la visibilidad. 54 55 if isnan(xinter) 56 warning(’NaN en intersección’) 57 end 58 59 if xinter>=x1 && xinter<=x2 && xinter>=x3 && xinter<=x4%Para que sea una intersección real tiene que dentro de als dos rectas, cualquiera de las condiciones saca del rango 60 % El punto está dentro del margen de las rectas, es que la 61 % intersección está produciendose en el rango de interés de las 62 % rectas 63 inter=1; %Intersección 64 else 65 %En caso de quese vaya del límite de alguna, no se ha producido. 66 inter=0; 67 end 68 else 69 %Si son iguales, es decir que son paralelas por lo que no puede haber punto de corte 70 inter=0; 71 xinter=nan(1); 72 % ’Bandera2’ 73 end 74 75 end Código A.7 Función cortapuntos. 94 Capítulo A. Códigos de Matlab Visible_DEF Una de las rutinas más importantes de todo el código, su propósito es el de establecer si la visión entre un punto del embudo y el nodo ii son visibles o ven su visión interrumpida. Para ello se establece un sistema de switch-case que va comprobando cada vez con mayor nivel de profundidad si el nodo es potencialmente visible o no. De esta manera pueden descartarse pasos dentro de esta rutina con mayores costes computaciones, como por ejemplo inter_visible_1. Para ello se hacen las siguientes comprobaciones en orden, siempre y cuando el caso se mantenga en ′proceso′: 1. Se revisa que el cono se encuentre dentro del cono de visión, si no lo está, no es visible dentro de ese cono. en el caso de que si lo fuera, por como se ha planteado la triangulación, en las anteriores o siguientes iteraciones de las bases de SPLIT definidas por vectorrecorre se encontraría. 2. Se revisa si el nodo se encuentra detrás de su padre en el árbol inferior, de ser así, es necesariamente no visible y además el modo de salida sería ′f ueraconoi1′ , como se mencionó en la sección del capítulo 2 correspondiente. Si no bloquea su visión potencialmente, entonces no es su padre en el árbol inferior, debido a la definición de embudo presentada. 3. Se revisa con inter_visible_1 si existen cortes con las aristas asociadas a los nodos con el resto de familias del embudo, o bien con las aristas asociadas a los polígonos cuyos embudos pertenezcan a la familia. Este paso es el más costoso computacionalmente pero el más efectivo a la hora de localizar cortes que pudieran encontrarse a causa de ausencia de embudos en la familia. 4. Se revisa si el padre en el árbol superior bloquea la visibilidad. Se identifica este nodo usando la función sucesores. Este método no es infalible, por lo que no se usa como criterio directo para definir la visibilidad, sino como soporte a inter_visible_1 . Esto se debe a que no existen garantías de que si el polígono al que pertenece el padre en el árbol superior es pequeño, el nodo de estudio no sea visible a la derecha del ya mencionado polígono, se adjunta la figura A.2 como soporte visual, que presenta un caso de ambigüedad de padre de embudo, que provocaría un falso negativo. En el artículo original se propone que la recurrencia de SPLIT debería encontrarlos al aplicar SPLIT a los nodos del ala derecha, pero en las pruebas del algoritmo implementado en geometrías complejas esta llamada pudiera no llegar a darse si faltara algún nodo. Con objeto de mantener baja la complejidad, se ha observado que era menos costoso computacionalmente comprobar los cortes que garantizar una familia completa, hecho comentado en las conclusiones del capítulo 3. Si todas estas condiciones se cumplen satisfactoriamente el nodo se valida como visible. 1function [VIS,modosalida,EVG]=visible_DEF(EVG,familiaembudo,indicevis,ii,nodox,nodoy) 2caso=’proceso’;%Variable que gestiona el switch case 3 4% Se comprueba que el nodo no pertenezca a la base del embudo 5if indicevis==1 || indicevis==length(familiaembudo) %|| EVG.Visibilidad(ii,familiaembudo( indicevis,1))==1 6caso=’nodobase’; 7end 8 9% Se comprueba que se encuentre dentro del cono de Visibilidad 10 if strcmp(caso,’proceso’) 11 %Comprobamos que se encuentre dentro del arco que forman x-v-y, si no está 12 %dentro, es necesariamente no visible. 13 14 p1=EVG.WPs(familiaembudo(1,1),:)’; %Punto x 15 p2=EVG.WPs(familiaembudo(indicevis,1),:)’; %Punto vis 16 p3=EVG.WPs(familiaembudo(length(familiaembudo(:,1)),1),:)’; %Punto y 17 p4=EVG.WPs(ii,:)’; %Punto v/ii A.3 Programas Children 95 Visión permitida ii y x Figura A.2 Ejemplo de caso de ambigüedad de embudo. 18 19 v1=p1-p4; 20 v2=p2-p4; 21 v3=p3-p4; 22 23 prodveca=v1(1)*v2(2)-v1(2)*v2(1); 24 prodvecb=v3(1)*v2(2)-v3(2)*v2(1); 25 26 if prodveca*prodvecb>0 27 caso=’fuera’; 28 if prodveca>0 29 %Es demostrable visualmente que está fuera del cono por la 30 %izquierda, siempre y cuando no se haya comprobado en una 31 %iteración anterior que si era visible 32 caso=’fueraconoi1’; 33 end 34 end 35 end 36 37 % Se comprueba que el padre en el árbol inferior de nodoVIS no bloquea su 38 % visibilidad 39 if strcmp(caso,’proceso’) 40 p1=EVG.WPs(familiaembudo(indicevis,2),:)’; %Punto padre vis 41 p2=EVG.WPs(familiaembudo(indicevis,1),:)’; %Punto vis 42 p4=EVG.WPs(ii,:)’; %Punto v/ii 43 44 v1=p1-p4; 45 v2=p2-p4; 46 47 prodvec=v2(1)*v1(2)-v2(2)*v1(1); 48 if prodvec<0 49 %Es decir, si v2 tiene que rotar CW, es que padrevis bloquea su 50 %visibilidad 51 caso=’fueraconoi1’; 52 end 53 end 54 55 % Se comprueba que no haya cortes con las aristas pertenecientes a los 56 % polígonos que tomen parte de la familia 57 if strcmp(caso,’proceso’) 58 if EVG.Caminos_prohibidos(familiaembudo(indicevis,1),ii)~=1%si no es un camino prohibido se computa directamente 59 [VIS,~,EVG]=inter_visible_1(indicevis,familiaembudo,EVG,ii); 60 else 61 VIS=false; 62 caso=’inter’; 63 end 96 Capítulo A. Códigos de Matlab 64 65 if VIS==true 66 caso=’nocorte’; 67 else 68 %Se revisa de forma adicional si está oculto detrás de un nodo que le postcede 69 posnodobloq=calculaCW_CCW(’CCW’, EVG.Matrizorden(familiaembudo(indicevis,1),:), EVG.Mat_ ind_pos_enMatrizorden(familiaembudo(indicevis,1),:),familiaembudo(indicevis,2),ii, familiaembudo(1:end,:)); 70 71 p1=EVG.WPs(familiaembudo(indicevis,1),:)’; 72 p2=EVG.WPs(posnodobloq,:)’; 73 p3=EVG.WPs(ii,:)’; %Punto v/ii 74 75 v1=p1-p3; 76 v2=p2-p3; 77 prodvec=v1(1)*v2(2)-v1(2)*v2(1); 78 79 if prodvec>0 80 p2=perpporpunto(EVG,nodox,nodoy,ii); 81 82 v1=p1-p3; 83 v2=p2-p3; 84 prodvec=v1(1)*v2(2)-v1(2)*v2(1); 85 86 if prodvec>0 87 %lo que corta la visibilidad 88 p1=EVG.WPs(familiaembudo(indicevis,1),:)’; 89 p2=EVG.WPs(familiaembudo(indicevis,2),:)’; 90 p3=EVG.WPs(familiaembudo(familiaembudo(:,1)==familiaembudo(indicevis,2),2),:)’; %Punto v/ii 91 92 v1=p1-p3; 93 v2=p2-p3; 94 prodvec=v1(1)*v2(2)-v1(2)*v2(1); 95 96 if prodvec<0 %Para que se de el caso de igualdad, tendría que haber saltado la condición de fueraconoi1 del principio 97 caso=’fueraconoi1’; 98 else 99 caso=’inter’; 100 end 101 else 102 % caso=’fueraconoi2’; 103 caso=’inter’; 104 end 105 else 106 caso=’fueraconoi1’; 107 end 108 end 109 end 110 111 switch caso 112 113 case ’nodobase’ 114 VIS=true; 115 modosalida=’VIS’; 116 case ’fuera’ 117 VIS=false; 118 modosalida=’fueraconod’; 119 case ’fueraconoi1’ 120 VIS=false; 121 modosalida=’fueraconoi1’; 122 case ’fueraconoi2’ 123 VIS=false; 124 modosalida=’fueraconoi1’; 125 case ’nocorte’ 126 VIS=true; 127 modosalida=’VIS’; 128 case ’inter’ 129 VIS=false; 130 modosalida=’inter’; A.3 Programas Children 97 131 otherwise 132 error(’Fallo Visibilidad’) 133 end 134 end Código A.8 Función Visible_DEF. Se adjunta la función perpporpunto, que como su nombre indica para una recta dada, calcula la intersección entre la recta y la línea que pasa por un punto externo con pendiente perpendicular a la recta. 1function [p_interseccion]=perpporpunto(EVG,nodox,nodoy,ii) 2 3p1=EVG.WPs(nodox,:)’; 4p2=EVG.WPs(nodoy,:)’; 5v=EVG.WPs(ii,:)’; 6% Vector director de la línea l1 7d = p2 - p1; 8 9% Vector normal a la línea l1 (perpendicular) 10 n = [-d(2), d(1)]; % Este es el vector normal 11 12 % Ecuaciones paramétricas de la línea l1 y de la normal 13 % l1: r(t) = p1 + t * d 14 % Normal: r(s) = v + s * n 15 16 % Resolver el sistema de ecuaciones: 17 % p1(1) + t*d(1) = vx + s*n(1) 18 % p1(2) + t*d(2) = vy + s*n(2) 19 20 A = [d(1), -n(1); d(2), -n(2)]; 21 b = [v(1) - p1(1); v(2) - p1(2)]; 22 23 % Resolver para t y s 24 sol = A\numberstyleb; 25 26 % Encontrar el punto de intersección usando t 27 t = sol(1); 28 p_interseccion = p1 + t * d; 29 30 end Código A.9 Función perpporpunto. inter_visible_1 inter_visible_1 es una subrutina dependiente de Visible_DEF que comprueba para los nodos de la familia, si existen cortes con las aristas que pertenecen a ella, usando la misma estrategia de detección de cortes empleada para la triangulación. Existen dos versiones posibles de esta función, uno que aplica esta definición directamente, y otra que extiende esta familia empleando todos los nodos asociados a los polígonos asociados a los nodos de la familia. La primera estrategia funcionó perfectamente para el barrido de los recintos convexos pero para el caso de estudio del método de Lee, [5], fue necesario aplicar la segunda. Ambas mantienen la complejidad en valores similares al ser aplicadas según la naturaleza de la geometría, por lo que se adjuntan ambas. En la versión estándar del código se añaden todos los vértices asociados al polígono del que ii forma parte, a los ya presentes en la familia. 1function [VIS,modosalida,EVG]=inter_visible_1(indicevis,familiaembudo,EVG,ii) 2nodo2poligono=EVG.nodo2poligono; 3poligono2nodo=EVG.poligono2nodo; 4 5nodox=familiaembudo(1,1); 6nodoy=familiaembudo(end,1); 98 Capítulo A. Códigos de Matlab 7nodoVIS=familiaembudo(indicevis,1); 8if nodox>nodoy 9nodo=nodox; 10 familiaembudo=EVG.nodosTOT_famvy_glob{nodo}; 11 else 12 nodo=nodoy; 13 familiaembudo=EVG.nodosTOT_famxv_glob{nodo}; 14 15 end 16 17 %Se añaden el resto de nodos del polígono a la lista de comprobación 18 19 listapol=[]; 20 21 poljj=nodo2poligono(ii); 22 if poljj==0 23 N=length(nodo2poligono(:,1)); 24 aristavector=[1,2,N-1,N,0]; %se saca el vector que contiene todas las aristas de EVG, atención contiene 0s 25 else 26 aristavector=poligono2nodo(poljj,:); %se saca el vector que contiene todas las aristas de EVG, atención contiene 0s 27 end 28 naristas=find(aristavector==0,1)-1; %Se identifica el números de aristas diferentes de 0 29 aristas=aristavector(1:naristas); %se extrae un vector de aristas limpio 30 familiaembudo=[familiaembudo;aristas’]; 31 32 33 34 [~,indiceunico]=unique(familiaembudo(:,1),’stable’); 35 familiaembudo=familiaembudo(indiceunico,:); 36 indicevis=find(familiaembudo(:,1)==nodoVIS,1); 37 38 VIS=true; 39 modosalida=’’; 40 41 for jj=[1:(indicevis-1),(indicevis+1):length(familiaembudo)] 42 poligono=EVG.nodo2poligono(familiaembudo(jj,1)); 43 inter1=0; 44 inter2=0; 45 46 if poligono==0 47 N=length(EVG.nodo2poligono(:,1)); 48 aristavector=[1,2,N-1,N,0]; %Se saca el vector que contiene todas las aristas de EVG, atención contiene 0s 49 else 50 aristavector=EVG.poligono2nodo(poligono,:); %se saca el vector que contiene todas las aristas de EVG, atención contiene 0s 51 end 52 53 naristas=find(aristavector==0,1)-1; %Identificamos el números de aristas diferentes de 0 54 aristas=aristavector(1:naristas); %se extrae un vector de aristas limpio 55 indpoljj=find(aristas==familiaembudo(jj,1),1); %buscamos en que posición se ubica el valor del nodo jj en las aristas 56 57 58 if indpoljj~=1 && indpoljj~=naristas 59 indarista1=[indpoljj+1,indpoljj]; 60 indarista2=[indpoljj-1,indpoljj]; 61 62 elseif indpoljj==1 63 indarista1=[indpoljj+1,indpoljj]; 64 indarista2=[naristas,indpoljj]; 65 66 elseif indpoljj==naristas 67 indarista1=[indpoljj-1,indpoljj]; 68 indarista2=[1,indpoljj]; 69 end A.3 Programas Children 99 70 71 aristaincidente1=aristas(indarista1); %Estos valores están en nodos 72 aristaincidente2=aristas(indarista2); %Estos valores están en nodos 73 74 line1=[EVG.WPs(ii,1),EVG.WPs(familiaembudo(indicevis,1),1);... 75 EVG.WPs(ii,2),EVG.WPs(familiaembudo(indicevis,1),2)]; 76 line2=[EVG.WPs(aristaincidente1(1),1),EVG.WPs(aristaincidente1(2),1);... 77 EVG.WPs(aristaincidente1(1),2),EVG.WPs(aristaincidente1(2),2)]; 78 79 [inter1,xinter1]=cortapuntos(line1,line2); 80 81 line2=[EVG.WPs(aristaincidente2(1),1),EVG.WPs(aristaincidente2(2),1);... 82 EVG.WPs(aristaincidente2(1),2),EVG.WPs(aristaincidente2(2),2)]; 83 [inter2,xinter2]=cortapuntos(line1,line2); 84 85 %inter vale 0 cuando no hay intersección 86 87 if ~(inter1==0 && inter2==0) 88 if inter1==1 && ~(abs(xinter1-EVG.WPs(ii,1))<1e-3) && ~(abs(xinter1-EVG.WPs( familiaembudo(indicevis,1),1))<1e-3)%Comprobamos que el punto de intersección no sea el nodo en si mismo 89 %si existe una intersección que no sea ii o nodo vis es necesariamente no visible 90 VIS=false; 91 % line1 92 % line2 93 % xinter1 94 modosalida=’inter1_intervisible’; 95 break 96 elseif inter2==1 && ~(abs(xinter2-EVG.WPs(ii,1))<1e-3)&& ~(abs(xinter2-EVG.WPs( familiaembudo(indicevis,1),1))<1e-3) 97 VIS=false; 98 modosalida=’inter2_intervisible’; 99 break 100 end 101 end 102 103 if jj==length(familiaembudo) 104 VIS=true; 105 modosalida=’’; 106 end 107 end 108 109 end Código A.10 Función inter_visible_1. A continuación se encuentra la versión alternativa de la extensión de la familia de embudos, nótese el bucle que recorre la familia comprobando si los nodos asociados al polígono han sido ya añadidos, comprobándolo con listapol. 1function [VIS,modosalida,EVG]=inter_visible_1(indicevis,familiaembudo,EVG,ii) 2nodo2poligono=EVG.nodo2poligono; 3poligono2nodo=EVG.poligono2nodo; 4 5nodox=familiaembudo(1,1); 6nodoy=familiaembudo(end,1); 7nodoVIS=familiaembudo(indicevis,1); 8if nodox>nodoy 9nodo=nodox; 10 familiaembudo=EVG.nodosTOT_famvy_glob{nodo}; 11 % familiaembudo2=EVG.nodosTOT_famxv_glob{nodo}; 12 13 else 14 nodo=nodoy; 15 familiaembudo=EVG.nodosTOT_famxv_glob{nodo}; 16 % familiaembudo2=EVG.nodosTOT_famvy_glob{nodo}; 17 18 end 19 100 Capítulo A. Códigos de Matlab 20 %Se añaden el resto de nodos del polígono a la lista de comprobación 21 22 listapol=[]; 23 for jj=1:length(familiaembudo) 24 poljj=nodo2poligono(familiaembudo(jj)); 25 indjj=find(listapol==poljj,1); 26 27 if isempty(indjj) 28 listapol=[listapol;poljj]; 29 30 31 if poljj==0 32 N=length(nodo2poligono(:,1)); 33 aristavector=[1,2,N-1,N,0]; %se saca el vector que contiene todas las aristas de EVG, atención contiene 0s 34 else 35 aristavector=poligono2nodo(poljj,:); %se saca el vector que contiene todas las aristas de EVG, atención contiene 0s 36 end 37 naristas=find(aristavector==0,1)-1; %Identificamos el números de aristas diferentes de 0 38 aristas=aristavector(1:naristas); %se extrae un vector de aristas limpio 39 familiaembudo=[familiaembudo;aristas’]; 40 41 end 42 43 end 44 [~,indiceunico]=unique(familiaembudo(:,1),’stable’); 45 familiaembudo=familiaembudo(indiceunico,:); 46 indicevis=find(familiaembudo(:,1)==nodoVIS,1); Código A.11 Versión alternativa de la ampliación familiaembudo para método de Lee. calculaCW_CCW:sucesores Se utilizan los nodos ordenados unos alrededor de otros que se obtuvieron en el preprocesado para convertir este cálculo en una consulta a una matriz y una búsqueda del siguiente nodo de la familia, filtrando fuera del subvector ordenado los nodos no presentes en la familia de embudos. 1%% Funcion Sucesores CW/CCW 2 3% nodoorden es el vector fila asociado al nodo i respecto al cual rotaremos 4% para buscar el sucesor, tanto CW como CCW del nodo j incidente (nodoinc) en el. 5% invnodoorden es el vector fila asociado al nodo i de Mat_ind_pos_enMatrizorden. 6% maxnodo es la manera de filtrar los nodos que todavía no se han 7% incorporado a la región P, estos no participarían en el cálculo de los 8% embudos. Sería el nodo que incorporamos a la región triangulada. 9 10 11 function sucesor=calculaCW_CCW(modo, nodoorden, invnodoorden, nodoinc,maxnodo,familiaembudo) 12 % indicenodoinc=invnodoorden(nodoinc); ha quedado obsoleto por culpa de que 13 % solo hay que usar los valores de pertenencientes a familiaembudo que son 14 % todos menores a maxnodo pero no son todos los menores a maxnodo 15 % nodoorden 16 [~, indicefiltrado]= ismember(nodoorden, familiaembudo(:,1)); 17 18 nodoorden=nodoorden(indicefiltrado>0); 19 indicenodoinc=find(nodoorden==nodoinc,1); 20 21 flag=0; 22 23 if ~isnan(indicenodoinc) 24 if strcmp(modo,’CCW’) 25 while flag==0 26 if indicenodoinc==length(nodoorden) 27 indicenodoinc=1; 28 if nodoorden(indicenodoinc)<maxnodo A.3 Programas Children 101 29 flag=1; 30 end 31 else 32 indicenodoinc=indicenodoinc+1; 33 % [nodoorden(indicenodoinc),maxnodo] 34 if nodoorden(indicenodoinc)<maxnodo 35 flag=1; 36 end 37 end 38 end 39 elseif strcmp(modo,’CW’) 40 while flag==0 41 if indicenodoinc==1 42 indicenodoinc=length(nodoorden); 43 if nodoorden(indicenodoinc)<maxnodo 44 flag=1; 45 end 46 else 47 indicenodoinc=indicenodoinc-1; 48 if nodoorden(indicenodoinc)<maxnodo 49 flag=1; 50 end 51 end 52 end 53 else 54 error(’Wrong input, please check modo variable input in calculaCW_CCW.’) 55 end 56 sucesor=nodoorden(indicenodoinc); 57 58 elseif strcmp(’CCW’,modo) 59 %Cualquier nodo cuyo padre en el lower tree sea si mismo, su padre en 60 %el arbol superior debe ser y 61 sucesor=familiaembudo(end,1); 62 else 63 sucesor=nodoinc; 64 end 65 66 end Código A.12 Función calculaCW_CCW. checkIfInsideTriangle Esa rutina implementa el cálculo para comprobar si para un linea de visión a la que se le va a aplicar una operación de SPLIT , el cono asociado está vacío. Para ello implementa la estrategia desarrollada en la figura 2.23. Para comprobar si existe un punto entre el triángulo formado por una linea de visión y un tercer punto externo o interno, se calculan las áreas utilizando la fórmula del determinante. A continuación se compara si las áreas que emplean el punto externo sumadas proporcionan sumadas valores mayores al triángulo original planteado, el punto se encontraría fuera del cono. De lo contrario de ser las áreas iguales, o ligeramente menores debido a errores numéricos, el punto se encontraría dentro del cono. 1function isInside = checkIfInsideTriangle(ii, nodoa, nodob, nodocheck, EVG) 2 3% Calcular las áreas 4% Área total del triángulo ii-nodoa-nodob 5areaTotal = triangleArea(ii, nodoa, nodob, EVG); 6 7% Área del subtriángulo formado por ii, nodoa y nodocheck 8area1 = triangleArea(ii, nodoa, nodocheck, EVG); 9 10 % Área del subtriángulo formado por ii, nodob y nodocheck 11 area2 = triangleArea(ii, nodob, nodocheck, EVG); 12 13 % Área del subtriángulo formado por nodoa, nodob y nodocheck 14 area3 = triangleArea(nodoa, nodob, nodocheck, EVG); 108 Capítulo A. Códigos de Matlab 241 end 242 243 % Dado que la condición de salida no es que se asignen todos los nodos, es 244 % posible que queden algunos que pudieron haber sido asignados pero quedaron 245 % pendientes al ser el mejor padre propuesto no válido, o que este nunca se 246 % llegara a asignar, por esto se plantea una segunda ronda. 247 248 flag_pruebagracia=1; 249 flag_padres=0; 250 padresposiblesfuturo=[]; 251 % seguimientonodos(isnan(seguimientonodos(:,2)),2)=0; 252 253 MATfuturasgeneraciones=zeros(length(vecaux(:,1)),1); 254 flagmodoejecucionfinal=1; 255 256 p3=WPs(vecaux(1,1),:)’;p4=WPs(vecaux(end,1),:)’; 257 258 while flag_pruebagracia==1||flagmodoejecucionfinal==1 259 flag_pruebagracia=0; 260 nodospendientes=seguimientonodos(seguimientonodos(:,2)==0,1); 261 262 %En esta ocasión el planteamiento será diferente, se buscará un padre para 263 %cada nodo 264 if flagmodoejecucionfinal==0 265 if flag_padres==0 266 padresposibles=familia(familia(:,4)==1,1);%Se extraen los nodos que pueden ser padres 267 flag_padres=1; 268 else 269 padresposibles=padresposiblesfuturo; 270 padresposiblesfuturo=[]; 271 end 272 else 273 padresposibles=familia(familia(:,4)==1,1);%Modo de ejecución final, todos los padres vá lidos son ahora posibles 274 end 275 276 for jj=1:length(nodospendientes) 277 nodoVIS=nodospendientes(jj); 278 p1=WPs(nodoVIS,:)’; 279 for kk=1:length(padresposibles) 280 nodopropuestopadre=padresposibles(kk); 281 p2=WPs(nodopropuestopadre,:)’; 282 cond=Visibilidad(nodoVIS,nodopropuestopadre); 283 284 if cond 285 v1=p2-p1;v2=p3-p1; 286 v3=p1-p2;v4=p3-p2; 287 prodvec=v1(1)*v2(2)-v1(2)*v2(1); 288 prodvec2=v3(1)*v4(2)-v3(2)*v4(1); 289 v5=p2-p4;v6=p1-p4; 290 prodvec3=v5(1)*v6(2)-v5(2)*v6(1); 291 292 if prodvec<=0 && prodvec2>=0 && prodvec3<=0 293 if nodopropuestopadre~=vecaux(1,1) 294 %Tecnicamente es una condición más desarollada que la 295 %de prodvec, más genérica pero solo valida para casos 296 %que padre y abuelo no son iguales 297 % p1=EVG.WPs(nodoVIS,:)’; %Nodo 298 % p2=EVG.WPs(nodopropuestopadre,:)’; %Padre Actual 299 nodoabuelo=familia(familia(:,1)==nodopropuestopadre,2); 300 pabuelo=WPs(nodoabuelo,:)’; %Nodo abuelo (Padre del padre) 301 302 v1=p1-p2; 303 v2=pabuelo-p2; 304 prodvecx=v1(1)*v2(2)-v1(2)*v2(1); 305 306 if prodvecx>0 307 condconvex=true; 308 else 309 condconvex=false; 310 end A.3 Programas Children 109 311 else 312 %Esta condición no tendría sentido si se analiza nodox como padre por lo que por 313 %defecto es verdadera 314 condconvex=true; 315 end 316 317 %Se va a meter una condición adicional, un nodo solo se 318 %puede añadir si queda CW de la arista 319 %nodopropuestopadre-siguiente aristapol 320 if condconvex==true 321 nodosconsecpadre=Consec_aristasrec(nodopropuestopadre,:); 322 323 listasec=nodoVIS; 324 325 if ~isempty(nodosconsecpadre) 326 for mm=1:length(nodosconsecpadre) 327 if nodosconsecpadre(mm)~=nodoVIS 328 listasec=[listasec,nodosconsecpadre(mm)]; 329 end 330 end 331 end 332 333 % %%%%%%%%%%% %ordenarlistasec 334 % pA=EVG.WPs(nodopropuestopadre,:); %Actua nodox 335 pA=p2; %Actua nodox 336 % pB=EVG.WPs(vecaux(end,1),:); %Actua nodoy 337 pB=p4; %Actua nodoy 338 339 vA=pB-pA; 340 listaauxsec=[listasec’,zeros(length(listasec’),1)]; 341 342 for mm=1:length(listasec) 343 nodomm=listasec(mm); 344 pC=WPs(nodomm,:)’; %Actua nodo a ordenar 345 vB=pC-pA; 346 ang=acos(dot(vB,vA)/(norm(vB)*norm(vA))); 347 prodveck=vB(1)*vA(2)-vB(2)*vA(1); 348 349 if prodveck<0 350 %Se busca ordenar rotando desde x y 351 %empezando desde y 352 ang=2*pi-ang; 353 end 354 listaauxsec(mm,:)=[nodomm,ang]; 355 end 356 listaauxsec=ordenalistasec(listaauxsec); 357 listasec=listaauxsec(:,1); 358 359 if length(listasec)==2 360 if nodoVIS==listasec(end) 361 %Caso D 362 condarista=true; 363 else %nodoVIS==listasec(1) 364 %Caso E 365 condarista=false; 366 end 367 elseif length(listasec)==3 368 if nodoVIS==listasec(1) 369 %CasoA 370 condarista=false; 371 elseif nodoVIS==listasec(2) 372 %Caso B: Este podría verse afectado en 373 %recintos no convexos 374 condarista=true; 375 elseif nodoVIS==listasec(end) 376 %caso C 377 condarista=true; 378 end 379 else 380 listasec 110 Capítulo A. Códigos de Matlab 381 error(’check condarista’) 382 end 383 else 384 condarista=false; 385 end 386 387 polnodopropuesto=EVG.nodo2poligono(nodopropuestopadre); 388 polnodox=EVG.nodo2poligono(vecaux(1,1)); 389 390 indfutgen=MAT_posnodosvecaux(nodoaux); 391 if flagmodoejecucionfinal==0 392 if condarista&&condconvex && nodopropuestopadre ~=vecaux(1,1)&& polnodopropuesto~=polnodox&& MATfuturasgeneraciones(indfutgen,1)==0 393 %Esto es un último check, para evitar que se 394 %introduzcan nodos antes de tiempo 395 [condpadrefuturo,posiblespadres]=futurasgeneraciones(EVG, nodospendientes,nodopropuestopadre,nodoVIS,vecaux(1,1)); 396 MATfuturasgeneraciones=extiendematriz(MATfuturasgeneraciones, indfutgen,posiblespadres); 397 398 elseif condarista&&condconvex && nodopropuestopadre ~=vecaux(1,1)&& polnodopropuesto~=polnodox&& MATfuturasgeneraciones(indfutgen,1)~=0 399 %Para el caso de que haya computados nodos, se 400 %revisará si el nodo que se propone en esta 401 %iteración es uno de los que se propuso como 402 %mejor padre 403 flag=0; 404 contador=1; 405 condpadrefuturo=false; 406 while flag==0 407 nodofuturo=MATfuturasgeneraciones(indfutgen,contador); 408 if nodofuturo==nodopropuestopadre 409 %Si el padre coincide con uno de los 410 %que se propuso en su momento se da 411 %como válido y corta la ejecución 412 [condpadrefuturo,posiblespadres]=futurasgeneraciones(EVG, nodospendientes,nodopropuestopadre,nodoVIS,vecaux(1,1)); 413 if ~isempty(posiblespadres) 414 MATfuturasgeneraciones=extiendematriz( MATfuturasgeneraciones,indfutgen,posiblespadres); 415 end 416 417 flag=1; 418 end 419 420 contador=contador+1; 421 if contador<=length(MATfuturasgeneraciones(1,:)) 422 if MATfuturasgeneraciones(indfutgen,contador)==0 423 flag=1; 424 end 425 else 426 %si se ha llegado al final de la 427 %matriz, se corta la ejecución del 428 %bucle 429 flag=1; 430 end 431 end 432 else 433 condpadrefuturo=true; 434 end 435 else 436 condpadrefuturo=false; 437 end 438 439 if cond && prodvec<=0 && prodvec2>=0 && prodvec3<=0 && condconvex && condarista && (condpadrefuturo||flagmodoejecucionfinal==1) 440 %Se busca en familia los hermanos 441 hermanosVIS=familia(familia(:,2)==nodopropuestopadre,1);%Nodos 442 if ~isempty(hermanosVIS) 443 hermanosVIS=[hermanosVIS;nodoVIS]; A.3 Programas Children 111 444 posnodoy=EVG.Mat_ind_pos_enMatrizorden(nodopropuestopadre,vecaux(end ,1)); 445 nodoorden=EVG.Matrizorden(nodopropuestopadre,:); %Se saca el vector de orden 446 nodoorden=[nodoorden(posnodoy+1:end),nodoorden(1:posnodoy)]; %Se pone a N último en orden CCW 447 nodoorden=nodoorden(end:-1:1); % Ahora N es primero en orden CW 448 [~, indicefiltrado]= ismember(nodoorden, hermanosVIS); 449 hermanosVIS=nodoorden(indicefiltrado>0); 450 451 indVIS=find(hermanosVIS==nodoVIS,1); %Una vez ordenados se buscan en que posición ha quedado nodoVIS 452 indhermanosVIS=find(familia(:,2)==nodopropuestopadre);%Indices 453 454 if indVIS==length(hermanosVIS) 455 %Si es el último de la lista se introduce después de 456 %los descendientes del último hermano 457 indiniciobusqueda=indhermanosVIS(end);%Esto es valido por que indhermanosVIS(end) se basa en familia 458 padresdescendencia=familia(indiniciobusqueda,1); 459 460 flag=0; 461 indseguimiento=indiniciobusqueda+1; 462 463 if indseguimiento>length(familia(:,1)) 464 flag=1; 465 indinserta=indiniciobusqueda; 466 end 467 468 while flag==0 469 if ~ismember(familia(indseguimiento,2),padresdescendencia)%si el padre del nodoindseguimiento no es de los permitidos entonces es que no es descendiente 470 flag=1; 471 %se rompe la ejecución, se ha identificado 472 %indiceinserta 473 indinserta=indseguimiento-1;%Se resta uno porque se ha pasado ya, se inserta sobre el anterior 474 break 475 end 476 477 if flag==0 478 padresdescendencia=[padresdescendencia;familia( indseguimiento,1)];%Si el nodo es hijo de los padres propuestos, entonces se añade a la lista de posibles padres futuros 479 indseguimiento=indseguimiento+1; 480 end 481 482 if indseguimiento==length(familia(:,1))+1 && flag==0 483 % Si se ha llegado al final de la lista se 484 % inserta en el último valor registrado 485 flag=1; 486 indinserta=indseguimiento-1; 487 end 488 end 489 490 familia=[familia(1:indinserta,:);nodoVIS nodopropuestopadre NaN (1,1) 0;familia(indinserta+1:end,:)]; 491 else 492 %En este caso solo es necesario insertarlo justo antes 493 %del siguiente hermano 494 nodohermanopos=hermanosVIS(indVIS+1);%hermanosVIS es una lista basada en nodos 495 indinserta=find(familia(:,1)==nodohermanopos,1)-1;%La posición del hermano posterior en familia -1 496 familia=[familia(1:indinserta,:);nodoVIS nodopropuestopadre NaN (1,1) 0;familia(indinserta+1:end,:)]; 497 end 498 499 else 500 %Si nodopropuestopadre no tiene hijos ya, se inserta el 112 Capítulo A. Códigos de Matlab 501 %nodo inmediatamente después del padre 502 indinserta=find(familia(:,1)==nodopropuestopadre,1); 503 familia=[familia(1:indinserta,:);nodoVIS nodopropuestopadre NaN(1,1) 0;familia(indinserta+1:end,:)]; 504 %FIN isempty 505 end 506 507 %Se hace una rápida comprobación de si nodoVIS se puede añadir a la lista de padres 508 509 padrevalido=EVG.MATposiblepadre(vecaux(1,1),nodoVIS); 510 511 if padrevalido==true 512 padresposiblesfuturo=[padresposiblesfuturo;nodoVIS]; %De esta manera en la siguiente iteración solo se trabajará con padres asignados en esta, mejorando el rendimiento 513 familia(familia(:,1)==nodoVIS,4)=1; 514 end 515 516 seguimientonodos(seguimientonodos(:,1)==nodoVIS,2)=1; 517 518 if flagmodoejecucionfinal~=1 519 flag_pruebagracia=1; 520 end 521 522 break % el nodo está asignado, se dejan de revisar padres 523 %Fin if es válido 524 end 525 %fin segundo if 526 end 527 %fin VIS 528 end 529 %Fin bucle padres 530 end 531 %Fin bucle nodos 532 end 533 534 if flag_padres==0 535 padresposibles=familia(familia(:,4)==1,1);%Se extraen los nodos que pueden ser padres 536 flag_padres=1; 537 else 538 padresposibles=padresposiblesfuturo; 539 padresposiblesfuturo=[]; 540 end 541 542 if isempty(padresposibles)&& flagmodoejecucionfinal~=1 543 flagmodoejecucionfinal=1; 544 elseif flagmodoejecucionfinal==1 545 %Si se llega con la bandera extendida, se corta la ejecución 546 flagmodoejecucionfinal=0; 547 end 548 end 549 550 familia=familia(:,1:2); 551 familia=[familia;vecaux(end,1),vecaux(1,1)]; 552 553 end Código A.14 Función add_familia2. A continuación se mostrarán las funciones auxiliares que emplea add_familia2 para cumplir sus funciones. futurasgeneraciones futurasgeneraciones pone en práctica parte de las ideas presentes para validar si un nodo puede presentar mejores padres o no. Por un lado, sólo considera aquellos nodos pendientes de ser asignados, ubicados en un sentido angular, desde ii entre nodopropuestopadre y nodox . Todas estas reglas aplican solo para los nodos que no pertenezcan al polígono de nodox , dado que se A.3 Programas Children 113 considerará que son los mejores candidatos al no poderse trazar rutas más cortas hacia nodox que por la frontera del polígono al que pertenece. Una vez con los nodos que se encuentren entre nodox y nodopropuestopadre localizados, se verifican las siguientes condiciones de forma secuenciada, descartando aquellos nodos que no pasen por el filtro, siendo la aplicación de estas condiciones cada vez más complejas. 1. Si existen varios nodos asociados a un mismo polígono en la lista de nodos, solo se considera aquel que forme parte de la cadena de tangentes que dirige a nodox , que por como está construida add_familia2 va a ser siempre nodopropuestopadre , por lo que se eliminan el resto de nodos que pertenezcan a este mismo polígono. De manera similar, no se va a llamar a la función futurasgeneraciones si nodopadrenextgen pertenece al polígono del nodox por lo que se eliminan también estos nodos. 2. Se comprueba que exista línea de visión directa entre los nodos restantes en la lista y nodoV IS , los que no cumplan esta condición se eliminan. De forma análoga a como se hacia en add_familia2, se comprueba que los nodos propuestos como padre bloqueen la visibilidad de nodoV IS respecto a nodox. 3. Se descartan aquellos nodos que no sean fatibles para ser padres usando la información del diagrama de visibilidad extendido, EV G.MAT posiblepadre(nodox,nodopadrenextgen). Los nodos que queden restantes son susceptibles a ser una mejor opción de padre para nodoV IS , se pasan como salida de la función y en add_familia2 se almacenará esta información en la matriz MAT f uturasgeneraciones, usando la función extiendematriz. 1function [condpadrefuturo,posiblespadres]=futurasgeneraciones(EVG,listanodos,nodopropuestopadre, nodoVIS,nodox) 2%listanodos en un vector columna 3caso=’proceso’; 4 5posnodopivote=EVG.Mat_ind_pos_enMatrizorden(nodoVIS,nodox); %Se rota desde nodoVIS a partir de nodox 6 7nodoorden=EVG.Matrizorden(nodoVIS,:); %Sacamos el vector deorden 8nodoorden=[nodoorden(posnodopivote:end),nodoorden(1:posnodopivote-1)]; %Se pone a nodox primero en ordenccw 9 10 %Para reducir el número de comprobaciones innecesarias se recortanlos 11 %potenciales padres futuros a revisar a aquells entre nodo x y 12 %nodopropuesto padre, esos no estarán en listanodos por lo que este 13 %recorte debe hacerse desde nodoorden 14 indnodopadre=find(nodoorden==nodopropuestopadre); 15 16 nodoorden=nodoorden(1:indnodopadre); 17 18 %Se filtra a partir de el vector recortado 19 [~, indicefiltrado]= ismember(nodoorden, listanodos); 20 listanodos=nodoorden(indicefiltrado>0); 21 22 if isempty(listanodos) 23 caso=’vacio’; 24 end 25 26 if strcmp(caso,’proceso’) 27 %Comprobación nodos mismo polígono-nodopropuestopadre 28 [a,b]=size(listanodos); 29 30 if b>a 31 listanodos=listanodos’; 32 end 33 34 listanodosextend=[listanodos,zeros(length(listanodos),1)]; 35 114 Capítulo A. Códigos de Matlab 36 for jj=1:length(listanodosextend(:,1)) 37 % nodoanalisis=listanodosextend(jj,1) 38 % EVG.nodo2poligono(nodoanalisis) 39 listanodosextend(jj,2)=EVG.nodo2poligono(listanodosextend(jj,1)); 40 end 41 polnodopropuesto=EVG.nodo2poligono(nodopropuestopadre); 42 polnodox=EVG.nodo2poligono(nodox); 43 44 listanodosextend(listanodosextend(:,2)==polnodopropuesto,:)=[];%El primer intento de engarce que funcione con un nodo de un polígono es el válido, el resto de nodos del polígono son sub óptimos por lo que hay que eliminarlos 45 listanodosextend(listanodosextend(:,2)==polnodox,:)=[];%No se consideran los nodos asociados a los poligonos que pertenezcan a nodox 46 47 if isempty(listanodos) 48 caso=’vacio’; 49 end 50 51 end 52 53 54 if strcmp(caso,’proceso’) 55 % Se revisa si los padres restantes son compatibles con nodoVIS 56 %Dado que falta información se usarán sólo dos criterios: Visibilidad y 57 %convexidad con nodox 58 listanodosextend(:,2)=0; %Se borra la información de los polígonos, ya no es relevante 59 60 for jj=1:length(listanodosextend(:,1)) 61 nodopadrenextgen=listanodosextend(jj,1); 62 63 condVIS=EVG.Visibilidad(nodoVIS,nodopropuestopadre); 64 65 p1=EVG.WPs(nodoVIS,:); %NodoVIS 66 p2=EVG.WPs(nodox,:); %Nodox: Servirá como referente para comprobar la convexidad 67 p3=EVG.WPs(nodopadrenextgen,:); %Nodopadrenextgen: se valorará como pivote 68 69 v1=p1-p3; 70 v2=p2-p3; 71 prodvec=v1(1)*v2(2)-v1(2)*v2(1); 72 73 %Para mantener la convexidad es necesario que nodovis encuentre el 74 %camino más corto pivotando por nodonextgen hacia nodox en sentido 75 %CCW 76 if prodvec>0 && condVIS 77 listanodosextend(jj,2)=1; 78 end 79 end 80 81 listanodosextend(listanodosextend(:,2)==0,:)=[];%Todos los valores que sean cero no son de interés, son padres predescartados 82 83 if isempty(listanodos) 84 caso=’vacio’; 85 end 86 end 87 88 if strcmp(caso,’proceso’) 89 % Comprobación test de paternidad: Esos nodos restantes pueden ser 90 % padres? 91 listanodosextend(:,2)=0; %Se borra la información de los polígonos, ya no es relevante 92 flag=0; 93 94 for jj=1:length(listanodosextend(:,1)) 95 nodopadrenextgen=listanodosextend(jj,1); 96 padrevalido=EVG.MATposiblepadre(nodox,nodopadrenextgen); 97 98 if padrevalido==true 99 %Se registra de los nodos restantes si hay alguno que pueda 100 %ser padre 101 % nodopadrenextgen A.3 Programas Children 115 102 listanodosextend(jj,2)=1; %este padre es válido 103 flag=1; 104 end 105 end 106 107 if flag==0 108 caso=’vacio’; 109 else 110 listanodosextend(listanodosextend(:,2)==0,:)=[];%Todos los valores que sean cero no son de interés, son padres predescartados 111 end 112 end 113 114 switch caso 115 case ’proceso’ 116 condpadrefuturo=false; 117 posiblespadres=listanodosextend(listanodosextend(:,2)==1,1); 118 case ’vacio’ 119 %si no hay padres posibles esta condición debe ser true 120 condpadrefuturo=true; 121 posiblespadres=[]; 122 end 123 124 end Código A.15 Función futurasgeneraciones. generalistanodos generalistanodos genera la lista de nodos por asignar ordenados respecto a nodopropuestopadre . Para ello toma el conjunto de nodos ordenados al completo y filtra respecto a la lista de nodos no asignada. Esta lista se reordena de manera que comienza en N , que su posición relativa es siempre atrás a la izquierda del nodox , por como se plantean los conos de visión siempre hacia occidente, y así garantizar que la secuencia es parte de una referencia correcta. 1function listanodos=generalistanodos(nodopropuestopadre,nodopivote,listanodos,EVG) 2%Se mira como se ordenan los nodos alrededor del padre propuesto 3 4posnodopivote=EVG.Mat_ind_pos_enMatrizorden(nodopropuestopadre,nodopivote); 5 6nodoorden=EVG.Matrizorden(nodopropuestopadre,:); %Se saca el vector deorden 7nodoorden=[nodoorden(posnodopivote+1:end),nodoorden(1:posnodopivote)]; %Se pone a N último en orden CCW 8nodoorden=nodoorden(end:-1:1); %Ahora N es primero en orden CW 9 10 [~, indicefiltrado]= ismember(nodoorden, listanodos); 11 listanodos=nodoorden(indicefiltrado>0); 12 end Código A.16 Función generalistanodos. extiende_matriz incorpora a la matriz, la lista de nodos que podrían conformarse como padres óptimos para un nodoV IS dado, este se identifica en la matriz con el indicador de posición pasado como argumento. Rellena con ceros para mantener consistente la concatenación de la matriz, dependiendo si el vector a introducir tiene más columnas que esta o no. 1function MAT=extiendematriz(MATBASE,indice,vector) 2 3[a,b]=size(MATBASE); 4[c,d]=size(vector);%Debe ser fila, por si acaso se añade revisión para hacerlo columna 5 6if c>d 7vector=vector’; 8[~,d]=size(vector); 116 Capítulo A. Códigos de Matlab 9end 10 11 %indice indica donde hay que insertar el valor 12 13 if b>d 14 %Si la matriz es más ancha que el vector se introduce ampliando el vector 15 MAT=[MATBASE(1:indice-1,:);vector,zeros(1,b-d);MATBASE(indice+1:end,:)]; 16 else 17 %En caso contrario hay que expandir la matriz 18 MAT=[MATBASE(1:indice-1,:),zeros(length(1:indice-1),d-b);vector;MATBASE(indice+1:end,:),zeros (length(indice+1:a),d-b)]; 19 end 20 21 end Código A.17 Función extiendematriz. ordenalistasec ordenalistasec es una función auxiliar del cálculo de listasec , que sirve para calcular condarista . Los argumentos del programa son la lista de nodos consecutivos en el polígono de nodopadrepropuesto y nodoV IS , junto con los ángulos dados respecto a la recta nodopropuestopadre −nodoy , para poder así extraer la información con los nodos dados en sentido creciente de las agujas del reloj. Como el número de puntos a ordenar va a ser siempre igual a 2 o a 3, esta ordenación se hace a mano, dado que la función sortrows consumía mucho tiempo en su ejecución. Es por ello que se ordena a mano a través de un sistema de sentencias if, siguiendo una lógica similar a de Bubble-sort, cambiando posiciones en los vectores con ciertas hipótesis adicionales que pueden brindarse debido a tener, en el peor de los casos, sólo tres componentes. 1function listaaordenar=ordenalistasec(listaaordenar) 2if length(listaaordenar(:,1))==2 3 4if listaaordenar(1,2)>listaaordenar(2,2) 5listaaordenar=[listaaordenar(2,:);listaaordenar(1,:)]; 6else 7listaaordenar=[listaaordenar(1,:);listaaordenar(2,:)]; 8end 9elseif length(listaaordenar(:,1))==3 10 if listaaordenar(1,2)>listaaordenar(2,2) 11 listaaordenar=[listaaordenar(2,:);listaaordenar(1,:);listaaordenar(3,:)]; 12 else 13 listaaordenar=[listaaordenar(1,:);listaaordenar(2,:);listaaordenar(3,:)]; 14 end 15 16 if listaaordenar(3,2)<listaaordenar(2,2) 17 listaaordenar=[listaaordenar(1,:);listaaordenar(3,:);listaaordenar(2,:)]; 18 19 if listaaordenar(1,2)>listaaordenar(2,2) 20 listaaordenar=[listaaordenar(2,:);listaaordenar(1,:);listaaordenar(3,:)]; 21 end 22 23 end 24 else 25 error(’Revisar listasec, longitud no aceptada.’) 26 end 27 28 end Código A.18 Función ordenalistasec. Con ello concluye el conjunto de códigos elaborados a lo largo de este trabajo. Índice de Figuras 1.1 Uso de productos vectoriales para determinar si una línea cruza por dentro de un polígono 4 1.2 Ejecución método de Lee para un punto de paso 5 2.1 Ejemplo práctico 8 2.2 Ejemplo de dominio con poligonos 10 2.3 Ejemplo práctico de triangulación respetando obstáculos 16 2.4 Sucesores CW y CCW 19 2.5 Extensiones CX y CCX 20 2.6 Dominio triangulado compuesto 21 2.7 Dominio triangulado compuesto: ii =622 2.8 Ejemplos de embudo en un dominio cualquiera 23 2.9 Regla de la goma elástica 23 2.10 Ambigüedad de padres factibles y muestra de padres no válidos 24 2.11 Familia de embudos para (x,y)25 2.12 Recorrido en preorden en sentido horario de la familia presente en la figura 2.11 26 2.13 Casos ejemplificados de los criterios aplicados 28 2.14 Reordenación de la familia partiendo de una base dada 28 2.15 Ejemplo de procesamiento de vectorpol, bloqueo de visión respecto a la envolvente 31 2.16 Ejemplo de la lógica de compruebasigno 32 2.17 Diferentes casos para vectorrecorre 37 2.18 vectorrecorre contiene solo un nodo: Caso xv y caso vy 38 2.19 Ejemplo de aplicación de SPLIT 41 2.20 Cálculo de la visibilidad mediante el empleo de las características de los embudos 44 2.21 Ejemplo de fallo en la demostración geométrica de la visibilidad con las propiedades de los embudos 45 2.22 Uso de recursividad SPLIT para detección de embudos 48 2.23 Uso de checkI f InsideTriangle para comprobar la posición relativa de un nodo externo 48 2.24 Justificación de añadir a f amiliavy pero no a f amiliaxv: Caso ′f ueraconoi1′50 2.25 Ejemplo de identificación uprima 51 2.26 Ejemplo de identificación nodor ynodoq 53 2.27 Identificación de relojes de arena para los casos de 2.26 53 2.28 Relevancia del uso de maxangulo 57 3.1 Ejemplo de región delimitada por meteorología adversa 68 3.2 Ejemplo esquemático de creación de recintos 71 3.3 Ejemplos de recintos complejos y sencillos 71 117 124 Bibliografía [11] Ligthart,L. P.,Yanovsky,F. J. yProkopenko,I. G. «Adaptive algorithms for radar detection of turbulent zones in clouds and precipitation». En: IEEE Transactions on Aerospace and Electronic Systems 39 (1 ene. de 2003), págs. 357 - 367. doi:10.1109/ TAES.2003. 1188918. [12] Kitzinger,J. «The Visibility Graph Among Polygonal Obstacles: a Comparison of Algorithms». En: Article. 2003. url:https:// api.semanticscholar.org/ CorpusID:9518166.