Día 131028_1 : Simulación : Detección de colisiones

DETECCIÓN DE COLISIONES

La detección de colisiones es útil para:

• Determinar la visibilidad de los objetos, haciendo la colisión de la visual del observador con los objetos que pueda haber delante de lo visible u oculto.
• Para determinar la velocidad y las posiciones de los objetos después de las colisiones.

En geometría computacional al usarse la coma flotante, se debe usar un mínimo de distancia de colisión:

• ε = 10E-4

Conceptos previos

Distancia de un punto a un segmento

• Se analiza si el punto esta dentro o fuera del segmento, según el producto dot.
• Si esta fuera del segemento, se calcula el módulo de la distancia punto a extremo del segmento.
• Si esta dentro del segmento, se hace igual que distancia punto a recta.

Distancia entre líneas

• La distancia mínima entre dos rectas que se cruzan es el módulo del vector w.
• El vector w pasa por los puntos P y Q de cada recta.
• La líneas están parametrizada por el punto por donde pasa y el vector dirección:
• P=P0 + ut
• Q=Q0 + vs
• D0 = Q0 - P0 
• Si son paralelas, la distancia es: D = |Do-Dot(Do,u)u|
• Si no son paralelas, la distancia es: D = Dot(Do,C)/|C| siendo C=Cross(U,V)

Intersección rayo-triángulo

• Muy usado en Ray-Tracing

1. Ver la orientación del triángulo, suponemos que la normal N mira hacia la fuente del rayo.
2. P es el punto de intersección del rayo con el triángulo: P =  Lo + L< (Po-Lo), N >
3. P puede estar dentro o fuera del triángulo -> Ver si esta dentro de todos los lados del triángulo.
4. Se dice que un punto "q" pertenece a un plano si: < q - Po, N > = 0, siendo Po un punto del plano.
5. Sustituyendo q por P: < Lo + L·d -Po, N > = 0
6. Con lo cual la distancia d = < (Po-Lo), N > / < L, N >

Código para la detección de colisiones:

//Codigo C para interseccion rayo-triangulo

//Punto a la izquierda

bool  pointLeft( vec3D q, vec3D p0, vec3D p1, vec3D N)

{

vec3D n=cross((p1-p0),N);

return( dot(q-p0,n) < 0.0);

}

//Interseccion rayo-plano

Vec3D  findRayIntersectPlane(Vec3D l0, Vec3D l,  Vec3D n, Vec3D p0, float &d) {

 // assuming n,l are normalized 

float denom = dot(n, l); d = dot(p0 - l0, n);

if (denom < 1e-6)  return l0;

return (l0 +(d/ denom )*l) 

}

//Interseccion con triangulo

bool  rayIntersectTriangle(Vec3D l0, Vec3D l, Vec3D t0, Vec3D t1, Vec3D t2){

float d=0.0;

Vec3D  n=normalize(cross(t1-t0,t2-t1));

Vec3D  p= findRayIntersectPlane( l0, l, n,t0,d);

return( (pointLeft(p,t0,t1,n) && pointLeft(p,t1,t2,n) pointLeft(p,t2,t0,n));

}

Colisión entre esferas

• La colisión se detecta si la distancia entre los centros es menor que la suma de los radios 
• distancia entre centros: d =||c2-c1||
• Colisionan si d < r1 + r2
• Se optimiza el cálculo si se usan los cuadrados: d^2 < (r1 + r2)^2
• Es más rápido, pero menos preciso
• No confundir con d^2 < r1^2 + r2^2, faltaría el termino 2r1·r2

Colisión entre capsulas

• La colisión se detecta si la distancia entre los segmentos es menor que los Radios de las esferas de los extremos: D < ( R1 + R2 )
• Es análogo a las distancia entre esferas.
• Se usa para colisiones entre objetos redondeados que deben rodar.

Colisión entre cajas (Bouding boxes)

• Dos cajas colisionan si existe una recta (2D) o si existe un plano (3D) de separación.
• Para resolver esto se usa el teorema del eje de separación.
• Dos poliedros convexos están separados si existe un plano de separación que sea:
• paralelo a alguna de sus caras
• paralelo al plano determinado por cada par de aristas eligiendo una en cada poliedro.
• En cuanto se de una de estas condiciones se considera que no se tocan.
• En consecuencia, basta con proyectar los vértices de los poliedros en lineas perpendiculares a los planos de separación para determinar separación.
• Según esto hay sólo 15 casos posibles que comprobar:
• Los 6 Planos paralelos a caras de ambas cajas:
• X1,  Y1, Z1
• X2, Y2, Z2
• Los 9 Planos paralelos a producto de aristas: (productor vectorial de todas las aristas)
• X1 x X2, X1 x Y2, X1 x Z2
• Y1 x X2, Y1 x Y2, Y1 x Z2
• Z1 x X2, Z1 x Y2, Z1 x Z2
• En cuanto uno estos casos se cumpla: NO HAY COLISIÓN.

Colisión entre Convex Hull

• Se emplea un método generalizado del plano de separación aplicado a las caras y aristas de los poliedros convexos. 
• No hay semiejes como en las cajas, que facilitan el cálculo.
• Hay que hacer la comprobación de los planos paralelos y del productos vectorial de aristas.
• El número de comprobaciones es n1+ n2+ ( m1 x m2 )
• Siendo ni = número de caras y mi = número de aristas.
• No se llegan a procesar todas, en cuanto una comprobación da positivo: NO HAY COLISIÓN.

Día 131023: Simulación: Voronoi y Delaunay

Diagrama de Voronoi

Dado un conjunto de puntos en 2D o 3D, se reparte el espacio en celdas de forma:

• Cada celda tiene una sola arista en común con la celda vecina
• Se emplea para rotura de sólidos rígidos (sin deformaciones), en materiales quebradizos.
• Las celdas son siempre convexas,
• Las celdas no acotadas (las de fuera) forman el Convex Hull.
• Cada vértice de la partición es el centro de un circulo determinado por 3 puntos distintos de P en regiones vecinas.

Algorirmo de division de una malla poligonal en celdas de Voronoi:

Tomo un punto "p" y recorro el resto de puntos "qi" , haciendo los segmentos que los unen.
Trazo la mediatriz, quedándome con el área del lado del punto "p".
Repito con otro punto "qi", pero sólo dividiendo el área que me quedó de la etapa anterior.

For each p, //Para cada punto de la malla

 PolyMesh* mesh = mesh0; //Define una malla

 For each q<>p do: //Para el resto de puntos distintos de p

  r= (q+p)/2   //toma el punto medio

  n= (q-p)/2  //calcula la normal

  splitPlane( r,n, mesh) //divide el plano

V(p)=mesh; //Guarda la malla

Triangulacion de Delaunay (fonéticamente «Deloné»)

Red de triángulos que cumple la condición de Delaunay:  La circunferencia circunscrita de cada triángulo de la red no debe contener ningún vértice de otro triángulo

Las aristas de los triangulos de Delaunay unen dos regiones vecinas de Voronoi.

Propiedades

• Cada nodo de Delaunay tiene asociada una celda de Voronoi
• El borde de Delaunay es el Convex Hull de la nube de puntos.
• En el interior del reciento de Delaunay no existen puntos
• Estos triángulos maximizan los ángulos de todas las soluciones posibles. Son triángulos muy buenos para cálculos robustos.
• En 3D Voronoy y Delaunay son tetraedros de ángulos maximizados.

Día 131021: Simulación: Envolturas convexas I

ENVOLTURAS CONVEXAS

• Los algoritmos geométricos son mas rápidos y sencillos con mallas convexas, debemos usarlas siempre que sea posible o como solución aproximada de un problema.
• Un polígono es convexo si y solo si para cada par de puntos en el polígono, el segmento que los une esta incluido en el polígono.
• El polígono mas sencillo y convexo es la L.

Distancia entre dos puntos

Para los cálculos es más eficiente usar la distancia al cuadrado, para simplificar.

A los lados de una recta

Pasa determinar si se esta dentro o fuera, a un lado o a otro de una recta se usa el producto escalar normalizado.

Dada una recta PoP1, calculo su normal con el producto Cross del vector y lo normalizo.
El punto "q" estará fuera o a la izquierda si el producto vectorial es menor que cero.
Es decir que el ángulo de los dos vectores P0-P1 con P0-q es negativo si esta fuera.

Distancia punto - recta

Es la menor de las distancias posibles.

La distancia es la proyección de w sobre v: w//

Día 131021: Simulación: Envolturas convexas Métodos Convex Hull

ENVOLTURAS CONVEXAS DE COLISION

Bounding Sphere:

• La esfera mas pequeña que envuelve al objeto
• Centro = punto medio de los vertices extremos: Sumatorio de coordenadas /  Nºvertices
• Radio = distancia maxima del centro a los vertices: (Distancia vertices)^2 / Nºvertices
• La colisión entre esferas es lo más fácil.

Bounding Box:

• Es la caja mas pequeña que envuelve el objeto
• Largo = distancia entre los vertices mas alejados en el eje X
• Ancho = distancia entre los vertices mas alejados en el eje Y
• Alto = distancia entre los vertices mas alejados en el eje Z
• Centro = punto medio de los vertices mas alejados

Bounding Capsule:

• La capsula mas pequeña que envuelve al objeto
• Lado, eje mas largo del bounding box
• Radio, segundo semi-eje mas largo del bounding box
• Se usa mucho en videojuegos aunque sea más lento, para rodar por escalera mejor.
• Se define como el conjunto de puntos distantes de un segmento.
• Es el menos preciso de aproxima, ya que tiene un sólo grado de ajuste.

Convex Hull

• Poliedro convexo mas pequeño que envuelve al objeto.
• Muy rápido en detectar colisiones.
• Construcción complicada
• Tipos de convex hull:
• GiftWrap - O(n*n)
• Incremental - O(nlogn)
• Quickhull - O(nlogn)
• Divide & Conquer - O(nlogn)

Proceso GiftWrap(2D) 

1. Se toma el punto más bajo
2. Se toma el siguiente punto con menor angulo con la horizontal
3. Se recorren los demás puntos y se toman los que tengan menor angulo con el segmento de los dos puntos anteriores.
4. Final cuando se llega al primer punto.

Proceso Quick Hull

1. Se toma el punto más a la derecha: h0
2. Se toma el punto más a la izquierda: h1
3. Se toma el conjunto de puntos a la derecha de h0h1
4. Se toma el conjunto de puntos a la izquierda de h0h1
5. Recorrer los puntos de la derecha y me quedo con el más alejado. h2
6. Recorrer los puntos de la izquierda y me quedo con el más alejado. h3
7. Repetir punto 5 y 6 con h0h2  h2h1 h0h3 y h3h1.
8. Final los conjuntos obtenidos están vacíos.

Proceso Incremental(2D) (procedimiento más usado, ya que es fácil de actualizar si cambia la nube)

1. Se parte de un triangulo inicial
2. Se toma el punto más abajo: h0
3. Se toma el punto más alejado de h0: h1
4. Se toma el punto más alejado del segmento [h0,h1]: h2
5. Se toma el subconjunto de puntos H0 dentro de [h0,h1,h2]
6. Se toma el subconjunto de puntos Q0 fuera de H0.
7. Ampliamos H0 con el punto qi =hi+1 más alejado de los hi: Hi+1=[h0, h1,...,hi,hi+1]
8. Rehago Qi con los puntos externos a Hi+1
9. Qi queda vacio.

Proceso Divide and Conquer Es un método peor y no se usa.

Hull( [p1,p2,.. Pn])
ordenar los puntos por su coordenada X,
q1,q2,.. qn
set Q1={qi, 1 to n/2}
set Q2={qi, n/2+1 to n}
Si size(Qi) > 3 -> Hi=Hull(Qi)  para i=1,2
H = merge( H1, H2)

Día 131016: Simulación: Mallas poligonales

MALLAS POLIGONALES

• ¿Por qué usar mallas poligonales? 
• Porque las tarjetas gráficas pintan únicamente triángulos y puntos.
• Son muy manejables, para calcular áreas y volúmenes.
• Con mallas poligonales se pueden representar cualquier sólido en 2D ó 3D.

Las mallas poligonales se miden por el LOD (Level Of Detail), representan sólidos a distintas resoluciones:

• Los polígonos son un conjunto finito de segmentos formados por una curva cerrada simple, sin nudos ni cruces.
• Cualquier polígono puede dividirse en triángulos, que es el polígono más sencillo.
• Para representar volúmenes se usan tetraedros, pero se usan poco.

Componentes: vértices, aristas, caras

Una malla se considera SOLIDA si es CERRADA:

• Sin aristas al aire. Cada arista debe pertenecer a dos caras y con todas las Mallas solidas: 
• Con todas las caras: Cada cara comparte sus aristas con otra cara.

Malla Poligonal en en C++ 

struct Vec3D //Definición estructura del vector en 3D

{

   double x,y,z; //Coordenadas de cada punto

   double w; //Coordenada proyectiva

};

struct PolyEdge //Definición estructura de la arista

{

   int  i0,i1; //Coordenada de sus vértices

   int  ti0,ti1; //Coordenadas de sus vértices de textura

   int  f0,f1; //Caras a las que pertenece

};

struct PolyFace //Definición estructura de cara

{

    int nVertex; //Número de vértices (puede tener >3)

    int iv[MAX_FACE_VERTEX]; //Array de índices de vértices

    int nEdges; //Número de aristas

    int ie[MAX_FACE_VERTEX]; //Array de Índices de aristas

    int iNormal; //Índice al array de normales

    long int flag; //Flag de uso general, activo on/off

    short imat; //Índice al material de la cara

};

struct PolyMesh //Definición estructura de malla
{
int nFaces; //Número de caras = Número de normales
PolyFace* faces; //Array estática de caras
int nEdges; //Número de aristas
PolyEdge* edges; //Array estática de aristas
int nVertex; //Número de vértices
Vec3D*  vertex; //Array estática de vértices.
int nTexVertex; //Nºcoord. de textura = Nºvertices
UVWTex* texVertex; //Array estática de coord.Textura
Vec3D*  Normals; //Array de normales = Nº de caras
long int flag; //Flag de uso general
}

Para hacer cálculos en 3D o hacer visualizaciones en 2D se usan coordenadas homogéneas:

| 1 0 0 0 | | x |   | x' |

| 0 1 0 0 | | y |   | y' |

| 0 0 1 0 | | z | = | z' |

| 0 0 0 1 | | w |   | w' |

Pasan de 3D a 2D directamente.

Adyacencias: Es necesario conocer:

• Caras a las que pertenece un vertice dado
• Aristas a las que pertenece un vertice dado
• Caras vecinas de una cara dada

Estas características para cada elemento se pre-calcula una vez definida la geometría.

Adyacencia en C++

Para obtener las caras adyacentes a un vértices:

void getVertexFaces( //No devuelve nada

 int iv, //Indice del vértice

 PolyMesh* mesh, //Array estático a la malla

 List<int>& faces //Lista de caras

 ) 

{ //Recorre las caras de la malla

 for( int i=0; i <mesh-> nFaces; i++) 

   { //Recorre los vertices de cada cara

  for( int j=0; j < faces[i]-> nVertex; j++)

  { //Si coincide nuestro vértice con el de la cara

   if( iv== faces[i]-> vertex[j])

   { //Lo agrega a la lista de caras adyacentes y va a otra cara

    faces.append(i);break;

   }

  }

 }

}

Para obtener las aristas adyacentes a un vértice:

void getFaceFaces( 
 int iface,
 PolyMesh* mesh,
 List<int>& faces
 )
 {//Recorre las caras de la malla
   for( int i=0; i <mesh-> nFaces; i++)
  {
      bool bfound=false; //Define comprobador y fija en false
   //Recorre los vértices de cada cara  
      for( int j=0; j < faces[i]-> nVertex; j++)
   {
    for( int k=0; k < mesh->faces[iface]->nVertex; k++)
    {
      if(faces[i]-> vertex[j] == mesh->faces[iface]-> vertex[k])
      {
                faces.append(i);bfound=true; break;
      }
    }
    if(bfound) break;
   }
 }

Para obtener las caras adyacentes a una cara:

void getFaceFaces(
 int iface,
 PolyMesh* mesh,
 List<int>& faces
 )
 {
  for( int i=0; i <mesh-> nFaces; i++)
  {
   bool bfound=false;
   for( int j=0; j < faces[i]-> nVertex; j++)
   {
    for( int k=0; k < mesh->faces[iface]->nVertex; k++)
    {
     if(faces[i]-> vertex[j] == mesh->faces[iface]-> vertex[k])
     {
      faces.append(i);bfound=true; break;
     }
    }
   if(bfound) break;
   }
  }

Cálculo con mallas poligonales
En simulación se define un valor infinitesimal equivalente al 0: epsilón = 10^-5 para

Cálculo del área

• Área de la malla = suma de las áreas de todos los triángulos.
• Área de un  triangulo con vértices ( a,b,c) es:
• V1= b-a
• V2 = c-a
• Área = ½* ||CrossProduct(V1, V2)||
• Funciona mejor con triángulos de ángulos abiertos.
• En C++:

//Devuelve un valor double a partir de una malla

double getArea( PolyMesh* mesh)

{

  double area=0.0; //Pone a cero el área al iniciar

 //Recorre todos las caras de la malla

 for( int i=0; i <mesh-> nFaces; i++)

 {

  PolyFace* face = mesh->faces[i]; //Toma una cara

  Vec3D a= mesh-> vertex[face->iv[0]]; //Vertice a

  Vec3D b= mesh-> vertex[face->iv[1]]; //Vertice b

  Vec3D c= mesh-> vertex[face->iv[2]]; //Vertice c

  area+=.5f*lenght(crossProduct(b-a, c-a)); //Area

 }

 return area; //Devuelve el valor del área

} 

Cálculo del Volumen

Problema: los triángulos tienen volumen cero por dentro. Se rellenan de dos formas:

• Metodo NAIF: tetraelizar la malla, y sumar los volúmenes. Coste computacional  muy alto.
• 
  
• Teorema de la DIVERGENCIA: Iguala una integral de volumen a una integral de superficie:

en nuestro caso:

queda resolver la integral en x de la superficie de todas las caras.

Partimos de una cara F definida por tres puntos P0,P1 y P2 y por las aristas Ei=Pi-P0 con i=1,2.

Parametrizando en u y v obtenemos x en función de u y v. 

Definimos la normal a la cara Nf en función de las aristas E1 y E2.

Haciendo una transformación de coordenadas obtenemos la función final.

En código C++:

//Devuelve un valor double a partir de una malla

double ComputeVolume( const PolyMesh* malla)

{

 const double fOneDiv6 = (double)(1.0/6.0); //Un sexto

 double volume=0.0; //Pone volumen a cero

 for(int i=0; i < malla->nfaces; i++) //Recorre las caras

 {

  //Trianguliza las caras con el metodo fan

  Vec3D p0 = malla->vertex[malla->faces[i].iv[0]];

  for(int j=0; j <malla->faces[j].nVertex – 2; j++)

  {

   Vec3D p1 = malla->vertex[malla->faces[i].iv[j+1]];

   Vec3D p2 = malla->vertex[malla->faces[i].iv[j+2]];

   //Obtiene la normal con el producto Cross de las aristas.

   Vec3D e1= p1 - p0;

   Vec3D e2= p2 - p0;

   Vec3D N = Cross( e1,e2);

   //Calcula los terminos F y N de la integral

  double  F1x = p0.x + p1.x+ p2.x;

  volume += N.x*F1x;

  }

 }

 volume *= fOneDiv6; //divide por seis

 return volume; //devuelve el volumen

}

Día 131014: Simulación por ordenador

Primera clase del Máster con Carlos Pegar:

• Copropietario de ThinKinetic junto con Pedro Ivan, producto: Pulldownit.
• En España se valora más la simulación para Ingenieria y para Medicina, fuera videojuegos.

SIMULACION: Reproducir la realidad con un ordenador basándose en leyes físicas.
EMULACIÓN: Reproduce la realidad sin basarse en leyes, es más difícil de parametrizar.

Areas de aplicación:

• Videojuegos y simuladores: la velocidad de calculo es lo mas importante(25 fps), la precisión del calculo es secundaria, se admiten defectos.
• VFX en cine y publicidad, velocidad y precisión son igualmente importantes, tanbien el control.
• Ingenieria, programas de CAD, la precision es lo mas importante.

El SOLVER es el encargado de calcular la simulación mediante la INTEGRACIÓN de las ecuaciones de estado y según la naturaleza de la simulación se utiliza:

• En sólidos rígidos las ecuaciones de estado son la ecuación de Newton F=m·a y las ecuaciones de acción-reacción.
• En sólidos que se deforman se emplean la ley de Hooke.
• En fluidos se usan las ecuaciones de Navier-Stokes.
• La luz se calcula con la ecuación de radiosidad:

Se resuelven mediante métodos numéricos que dan soluciones aproximadas.

Partición del espacio:  Hay que acotar el espacio y parcelarlo de forma que los cálculos sean mas sistemáticos y rápidos:

• Solidos: cajas, esfera,  convex-hull y ABB-trees (para formas complejas con agujeros).
• Fluidos: Grids (cuadrícula con coordenadas)
• Luz: KD- trees (como BS-tree en videojuegos)

Validación: Acotar el error, verificar corrección del modelo, es la fase mas difícil.

• La coma flotante limita la precisión del modelo.
• Los métodos numéricos dan una solución aproximada.
• En efectos especiales, se busca que se parezca a la realidad, no que sea exacto.
• Física: Conservación de la energía y el momento (péndulo convergencia a pararse)
• Matemáticas: Métodos numéricos, métodos estadísticos.
• Se verifica con el teorema de la energía o del momento.

Fundamentos Matemáticos:

• Geometría proyectiva.
• Geometria computacional.
• Métodos numéricos.
• Estadística y probabilidad.

Geometria proyectiva:

• Representar objetos en el espacio 3D, mover, rotar y escalar los objetos.
• Proyecciones 2D, ortogonal, perspectiva, cámaras.
• Esta materia la da Álvaro.

Geometria computacional: (a continuación de la geometría proyectiva)

• Recorrer la superficie de los objetos.
• Las tarjetas gráficas sólo pintan triángulos y puntos, por eso se usan mallas de triángulos.
• Las mallas son manejables, para el cálculo de áreas y volúmenes.
• Las mallas permiten representar cualquier sólido en 2D ó 3D.
• Calcular envolventes para los objetos : Conejo de Stanford.
• Particionar objetos.
• Detectar contacto entre distintos objetos.

Métodos numéricos: Es lo más difícil.

• Resolver ecuaciones diferenciales.
• Aproximar funciones.

Estadistica y probabilidad:

• Simplificar grandes volúmenes de datos.
• Acotar el error.
• Aproximar funciones complejas.

Principios:

• LEY DE HOOKE: 
• Para un resorte:
•  
• Relaciona la fuerza F ejercida en el resorte con el alargamiento "delta" producido.
•   Es la constante elástica.
• 
  
• 
• Establece que la deformación de un material elástico "ε" es proporcional a la fuerza aplicada F.
• Siendo :
• ΔL : Alargamiento longitudinal
• L : Longitud original 
• E : Modulo de Young o de elasticidad 
• A : Sección transversal de la pieza estirada.
• Se aplica a los sólidos elásticos hasta un límite denominado límite de elasticidad.
•  
• Esta es la ecuación diferencial del muelle. Si se integra para todo x, se obtiene como ecuación de onda unidimensional que describe los fenómenos ondulatorios:
• La velocidad de propagación de las vibraciones en un resorte se calcula como:
• 
• Para un sólido elástico:
• En la mecánica de sólidos deformables elásticos la distribución de tensiones es mucho más complicada que en un resorte o una barra estirada sólo según su eje.
•  La deformación en el caso más general necesita ser descrita mediante un tensor de deformaciones mientras que los esfuerzos internos en el material necesitan ser representados por un tensor de tensiones. Estos dos tensores están relacionados por ecuaciones lineales conocidas por ecuaciones de Hooke generalizadas o ecuaciones de Lamé-Hooke:
• 
• Caso unidimensional: las deformaciones o tensiones en direcciones perpendiculares a una dirección dada son irrelevantes o se pueden ignorar:     
• La ecuación anterior se reduce a: , donde E es el módulo de Young.
• Caso tridimensional: Para caracterizar el comportamiento de un sólido elástico lineal e isótropo se requieren además del módulo de Young (E) otra constante elástica: el coeficiente de Poisson ().
• Por otro lado, las ecuaciones de Lamé-Hooke para un sólido elástico lineal e isótropo pueden ser deducidas del teorema de Rivlin-Ericksen, que pueden escribirse en la forma:
• 
• 
• 
• En forma matricial, en términos del módulo de Young y el coeficiente de Poisson como:
• 
• Las relaciones inversas vienen dadas por:
•