UN NUEVO ALGORITMO EN LA TÉCNICA DE VELOCIMETRÍA POR ...

INSTITUTO POLITÉCNICO NACIONAL

CENTRO DE INVESTIGACIÓN EN COMPUTACIÓN

UN NUEVO ALGORITMO EN LA TÉCNICA DE VELOCIMETRÍA POR IMÁGENES DE PARTÍCULAS

T E S I S

QUE PARA OBTENER EL GRADO DE

MAESTRO EN CIENCIAS DE LA COMPUTACIÓN

P R E S E N T A

ING. MARISOL LÓPEZ HINOJOSA

DIRECTOR DE TESIS

DR. GERMÁN GONZÁLEZ SANTOS

CODIRECTOR DE TESIS DR. ADOLFO GUZMÁN ARENAS

MÉXICO D.F. 2006

Agradecimientos.

A los siguientes laboratorios del IPN: Reología de la Escuela Superior de Física y Matemáticas,

Velocimetría de la Escuela Superior de Ingeniería Química e Industrias Extractivas y Velocimetría de

la Escuela Superior de Ingeniería Mecánica y Eléctrica, por las imágenes otorgadas para las pruebas

de esta tesis.

Al Consejo Nacional de Ciencia y Tecnología (CONACyT) por la beca otorgada para realizar mis

estudios de maestría.

ÍNDICE GENERAL

Índice de figuras IX

Índice de tablas XIII

1. Introducción a la Técnica de PIV. 1

1.1. Antecedentes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2. Planteamiento del problema. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.3. Objetivos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.4. Solución propuesta. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.5. Aportaciones. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.6. Estructura de la tesis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2. La Técnica de Velocimetría por Imágenes de Partículas 9

2.1. Partículas trazadoras. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.2. Fuente de luz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.3. Captura y almacenamiento de imágenes. . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

Índice general

2.4. Evaluación de las imágenes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.5. Procesamiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3. Diseño e Implementación 29

3.1. Etapa de pre-procesamiento. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.2. Etapa de procesamiento. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

3.3. Etapa de post-procesamiento. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

4. Pruebas y Resultados 51

4.1. Descripción de los programas de PIV evaluados. . . . . . . . . . . . . . . . . . . . . . . . 52

4.2. Resultados. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

5. Conclusiones y Trabajo Futuro. 65

Referencias 67

A. Características de los componentes de PIV. 71

A.1. Partículas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

A.2. Láser . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

A.3. Cámara digital . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

A.4. Sistemas de PIV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

B. Diagramas 79

B.1. Diagramas de flujo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

B.2. Diagrama de estados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

VI Un nuevo algoritmo en la técnica de velocimetría por imágenes de partículas

Índice general

B.3. Pantallas de la interfaz gráfica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

C. Código 89

C.1. Programa principal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

C.2. Pre-procesamiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

C.3. Procesamiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

C.4. Post-procesamiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

C.5. Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

D. Instalación del sistema 119

D.1. Matlab . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

D.2. Requerimientos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

D.3. Formatos para imágenes de entrada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

D.4. Funciones utilizadas incluidas en Matlab. . . . . . . . . . . . . . . . . . . . . . . . . . . 121

E. Manual de usuario 125

E.1. Paso 1: Datos generales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

E.2. Paso 2: Datos de pre-procesamiento. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126

E.3. Paso 3: Datos de procesamiento. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127

E.4. Paso 4: Datos de post-Procesamiento. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129

F. Pruebas Realizadas. 131

Un nuevo algoritmo en la técnica de velocimetría por imágenes de partículas VII

ÍNDICE DE FIGURAS

1.1. Sistema de velocimetría por imágenes de partículas . . . . . . . . . . . . . . . . . . . . . 2

1.2. (a) (b): Imágenes de partículas, (c) (d): Mapa de vectores de velocidad . . . . . . . . . 4

2.1. (a) Imagen con una exposición, (b) Imagen con doble exposición, (c) Imagen con múltiple

exposición. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.2. (a) Dos imágenes con una exposición, (b) Múltiples imágenes con una exposición. . . . . 13

2.3. Representación de las imágenes de PIV. . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.4. Representación digital de una imagen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.5. Los tres tipos de densidad de partículas en las imágnes: (a) Baja (PTV), (b) Media

(PIV), y (c) Alta (LSV). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.6. (a) Imagen de partículas tomada en el tiempo t0, (b) Imagen de partículas tomada en

el tiempo t1, (c) Desplazamiento de las partículas en el intervalo t0-t1. . . . . . . . . . . 17

2.7. Imagen de partículas de PIV con la amplificación de una parte de la imagen. . . . . . . 18

2.8. Dos imágenes de PIV divididas en zonas de interrogación. . . . . . . . . . . . . . . . . . 18

2.9. Correlación de dos zonas de interrogación. . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.10. (a) y (b) Imágenes de partículas, (c) y (d) Plano de correlación cruzada. . . . . . . . . . 21

Índice de figuras

2.11. Proceso para obtener el mapa de vectores de desplazamiento. . . . . . . . . . . . . . . . 23

2.12. (a) (b) Imágenes de partículas 512 por 512 píxeles, (c) mapa de vectores con zonas de

256 por 256 píxeles (d) mapa de vectores con zonas de 32 por 32 píxeles. . . . . . . . . . 23

2.13. (a) Imagen de partículas tomada en el tiempo t0 , (b) Imagen de partículas tomada en

el tiempo t1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.14. Mapa de vectores con zonas de (a) 128 x 128 píxeles (b) 64 x 64 píxeles. . . . . . . . . . 24

2.15. Mapa de vectores con zonas de (a) 32 x 32 píxeles (b) 16 x 16 píxeles. . . . . . . . . . . 26

2.16. Mapa de vectores con zonas de (a) 8 x 8 píxeles (b) 4 x 4 píxeles. . . . . . . . . . . . . . 26

3.1. Diagrama de estados de datos del sistema. . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.2. Imagen original con áreas a recortar. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3.3. Imagen recortada. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

3.4. Imagen con objetos no deseados estáticos. . . . . . . . . . . . . . . . . . . . . . . . . . . 32

3.5. Imagen con el contraste ecualizado. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.6. Imágenes de PIV con objetos fijos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

3.7. Máscara de objetos estáticos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3.8. Imagen de partículas sin objetos estáticos. . . . . . . . . . . . . . . . . . . . . . . . . . . 36

3.9. Imagen de partículas con objetos no deseados (burbujas). . . . . . . . . . . . . . . . . . 37

3.10. Filtro de Mediana. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

3.11. Imagen con elementos ajenos sin partículas. . . . . . . . . . . . . . . . . . . . . . . . . . 39

3.12. Imagen de partículas pre-procesada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

3.13. (a)Imagen sin traslape, (b) imagen con traslape del 50% . . . . . . . . . . . . . . . . . . 41

3.14. Proceso iterativo con multiresolución. Reducción del 50% del tamaño de las zonas de

interrogación con cada iteración. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

X Un nuevo algoritmo en la técnica de velocimetría por imágenes de partículas

Índice de figuras

3.15. Mapa de vectores sin pre-procesamiento. . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

3.16. Mapa de vectores con pre-procesamiento. . . . . . . . . . . . . . . . . . . . . . . . . . . 44

3.17. (a) Mapa con vectores incorrectos (b) mapa sin vectores incorrectos (c) mapa con vec-

tores sustituidos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

3.18. Matriz de vectores. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

3.19. Mapa con varios vectores a reemplazar. . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

3.20. Mapa de vectores sin pre-procesamiento, con validación. . . . . . . . . . . . . . . . . . . 49

3.21. Mapa de vectores con pre-procesamiento, con validación. . . . . . . . . . . . . . . . . . . 50

3.22. Mapa de vectores sin pre-procesamiento, con validación y sustitución. . . . . . . . . . . 50

3.23. Mapa de vectores con pre-procesamiento, con validación y sustitución. . . . . . . . . . . 50

4.1. Imagen de partículas 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

4.2. Imagen de partículas 2, tomada 25 ms después de la primera imagen . . . . . . . . . . . 53

4.3. Mapa de vectores resultante del modelo. . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

4.4. Parábola promedio del modelo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

A.1. (a) Foto tomada en el 2004 de un experimento en el laboratorio de Velocimetría, (b)

Buza de acrílico por donde sale el agua. . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

A.2. Óptica para generar el plano de luz láser. . . . . . . . . . . . . . . . . . . . . . . . . . . 74

B.1. Algoritmo de pre-procesamiento. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

B.2. Algoritmo de procesamiento. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

B.3. Algoritmo de post-procesamiento. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

B.4. Diagrama de estados de pre-procesamiento. . . . . . . . . . . . . . . . . . . . . . . . . . 83

Un nuevo algoritmo en la técnica de velocimetría por imágenes de partículas XI

Índice de figuras

B.5. Diagrama de estados de procesamiento. . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

B.6. Diagrama de estados de post-procesamiento. . . . . . . . . . . . . . . . . . . . . . . . . . 85

B.7. Diagrama de estados del despliegue de resultados. . . . . . . . . . . . . . . . . . . . . . . 85

B.8. Diseño de pantalla de datos generales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

B.9. Diseño de pantalla de datos de pre-procesamiento. . . . . . . . . . . . . . . . . . . . . . 86

B.10.Diseño de pantalla de datos de procesamiento. . . . . . . . . . . . . . . . . . . . . . . . . 87

B.11.Diseño de pantalla de datos de post-procesamiento. . . . . . . . . . . . . . . . . . . . . . 87

E.1. Pantalla de datos generales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126

E.2. Pantalla de datos de pre-procesamiento. . . . . . . . . . . . . . . . . . . . . . . . . . . . 127

E.3. Pantalla de selección de máscara de procesamiento. . . . . . . . . . . . . . . . . . . . . . 128

E.4. Pantalla de datos de procesamiento. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129

E.5. Pantalla de datos de post-procesamiento. . . . . . . . . . . . . . . . . . . . . . . . . . . 130

XII Un nuevo algoritmo en la técnica de velocimetría por imágenes de partículas

ÍNDICE DE TABLAS

3.1. Estimadores de tres puntos para determinar el desplazamiento de los datos de correlación

a nivel subpíxel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

4.1. Resultados aplicación PIV propuesta. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

4.2. Resultados Dantec. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

4.3. Comparación de desviaciones estándar por regiones. . . . . . . . . . . . . . . . . . . . . 62

4.4. Resultados Matpiv. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

A.1. Partículas que pueden usarse en la visualización de flujos líquidos. . . . . . . . . . . . . 71

A.2. Partículas que pueden usarse en la visualización de flujos gaseosos. . . . . . . . . . . . . 72

A.3. Propiedades de un sistema de láser Nd:YAG para PIV. . . . . . . . . . . . . . . . . . . . 74

A.4. Cámaras digitales para PIV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

A.5. Sistemas PIV comerciales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

D.1. Formatos de imágenes soportados por MATLAB. . . . . . . . . . . . . . . . . . . . . . . 122

Resumen.

La velocimetría por imágenes de partículas (PIV) es una técnica basada en el procesamiento digital

de imágenes, la cual permite obtener, en forma indirecta, el campo de velocidades de un fluido en

movimiento.

El presente trabajo está dedicado al diseño e implementación de un nuevo algoritmo en la técnica

de PIV, el cual mejora los resultados obtenidos con el algoritmo clásico al obtener con mayor exactitud

y resolución espacial el campo de velocidades. El algoritmo propuesto está dividido en tres etapas:

pre-procesamiento, procesamiento, y post-procesamiento. En la primera etapa se aplican técnicas de

procesamiento de imágenes para mejorar las imágenes de entrada. En la etapa de procesamiento se

aplica la correlación cruzada de imágenes en forma iterativa para distintas resoluciones de vectores,

con el propósito de aumentar la calidad de resultados. Finalmente, en la tercera etapa, son validados

los resultados de la etapa previa y los vectores incorrectos son sustituidos.

Con el algoritmo propuesto se calcularon los campos de velocidad para diferentes tipos de flujo y

se obtuvieron mejores resultados que con la mayoría de los sistemas existentes de su tipo.

El algoritmo fue implementado en Matlab y se le agregó una interfaz gráfica para que el usuario

pueda elegir las diferentes opciones de procesamiento.

Abstract.

Particle Image Velocimetry (PIV) is a non-intrusive technique to determine the flow velocity field

by image processing.

This work describes the design and implementation of a new PIV algorithm that improves the results

obtained with the classic algorithm by getting a more exact velocity field with greater spatial resolution.

The new algorithm is divided in three phases: pre-processing, processing and post-processing. In the

first phase, the input images are enhanced by applying image processing techniques. In the processing

phase, an iterative cross-correlation with different vector resolution is used to improve the quality of

the results. Finally, in the post-processing phase, the results are validated and the incorrect vectors

are substituted.

Different fluid velocity fields were estimated with the proposed algorithm having obtained better

results than most of the existing similar programs.

The algorithm was implemented in Matlab with a graphic user interface to enable the choice of

different processing options.

CAPÍTULO 1

Introducción a la Técnica de PIV.

1.1. Antecedentes.

El procesamiento digital de imágenes ha mejorado la manera en que los fluidos y sus características

pueden ser estudiados. A través de este procesamiento se pueden obtener de forma automática, rápida

y precisa, características importantes del fluido en cuestión, como densidad, temperatura, velocidad y

estructuras que se forman cuando el fluido está en movimiento.

El determinar la velocidad de un flujo de fluidos es un problema importante en diferentes áreas

de la ingeniería, la industria y la ciencia. Las técnicas clásicas1 emplean diferentes dispositivos que se

introducen al fluido, para medir la velocidad en un punto durante un intervalo de tiempo.

La velocimetría por imágenes de partículas (PIV por sus siglas en inglés) es una técnica no invasiva

para obtener, mediante el procesamiento digital de imágenes, el campo de velocidades en dos dimensio-

nes de un fluido en movimiento. Esta técnica es una herramienta útil en el estudio de flujos de fluidos

en diferentes geometrías y para diferentes regímenes de flujo; desde estacionario hasta turbulento.

El punto de partida de esta técnica consiste en agregar pequeñas partículas al fluido bajo estudio.

Se asume que estas partículas se mueven con la velocidad y dirección del fluido. El fluido y el recipiente

que lo contenga (en caso de estar confinado), deben ser transparentes para que las partículas puedan

ser visibles. Una región del fluido, donde se desea obtener el campo de velocidades, es iluminada con un

plano de luz láser y una cámara digital es colocada en forma perpendicular a este plano. Las partículas

que atraviesan el plano de luz son iluminadas al menos dos veces en un periodo muy corto de tiempo.

El tiempo entre dos pulsos de luz dependerá de la velocidad media del flujo. La cámara y la fuente de

1Algunas de estas técnicas son la Anemometría por Hilo Caliente (Hot Wire Anemometry) y la Anemometría Láser

Doppler (Laser Doppler Anemometry).

1.1. Antecedentes.

luz son sincronizadas por medio de una computadora. Cuando el fluido es iluminado, la luz dispersada

por las partículas que atraviesan el plano de luz, es capturada por la cámara en una imagen o en

una secuencia de imágenes que son almacenadas en la computadora. Las imágenes generadas son la

representación bidimensional del volumen formado por la región del fluido iluminada por el plano láser.

El grosor de este plano es menor a un milímetro. Para poder determinar la velocidad es necesario contar

con las imágenes de al menos dos exposiciones2, el motivo se explica en el siguiente capítulo.

En la Figura 1.1 se muestra un sistema básico de PIV. Se utiliza un láser como fuente de luz, una

cámara digital para capturar las imágenes y una computadora para procesar los datos.

Figura 1.1: Sistema de velocimetría por imágenes de partículas

Una vez que se tienen las imágenes, en caso de ser necesario, éstas pueden ser mejoradas para

optimizar su procesamiento y corregir ciertos defectos experimentales. Después estas imágenes son

procesadas en una computadora mediante un algoritmo de PIV. Para su evaluación todas las imágenes

2La exposición es la cantidad de luz que llega al sensor de la cámara durante un intervalo de tiempo.

2 Un nuevo algoritmo en la técnica de velocimetría por imágenes de partículas

Capítulo 1. Introducción a la Técnica de PIV.

son divididas de igual forma, en pequeñas áreas rectangulares llamadas “zonas de interrogación”. Se

asume que todas las partículas dentro de cada una de estas zonas se mueven homogéneamente entre

las tomas realizadas.

La velocidad de las partículas no se obtiene de manera individual sino en promedio. Para este

propósito se emplean métodos estadísticos para determinar el desplazamiento promedio de las partícu-

las de cada zona de interrogación. Con el desplazamiento y el tiempo transcurrido entre un par de

imágenes, se estima la magnitud y dirección de la velocidad local del fluido en esa zona. Este proceso

es repetido para todas las zonas de interrogación y se construye un mapa con todos los vectores de

velocidad calculados.

A continuación se describen brevemente las principales características de la técnica de PIV.

El proceso es rápido. Determina el campo de velocidades en tiempo real, lo que permite encontrar

estructuras en flujos no estacionarios.

Se pueden determinar dos componentes de velocidad debido la iluminación por medio de un plano

(dos dimensiones).

PIV es una técnica no invasiva, es decir, no se necesita introducir ningún dispositivo en el fluido.

Para medir la velocidad se utilizan pequeñas partículas, pero estas no perturban el flujo de forma

significativa a diferencia de los sensores utilizados en otras técnicas.

La medición de velocidad del flujo del fluido se hace de forma indirecta, pues realmente se calcula

la velocidad de las partículas que se encuentran en el flujo. Es necesario revisar cuidadosamente

las partículas a utilizar en cada experimento, con el fin de que estas tengan una distribución

homogénea y sigan adecuadamente el flujo, sin perturbarlo de forma significativa.

La resolución espacial3 del mapa de vectores de velocidad obtenido con PIV depende principal-

mente de la densidad de partículas y del tamaño de la imagen. Las imágenes de PIV tienen una

densidad de partículas y resolución espacial de resultados media. Debido a que para su procesa-

miento las imágenes son divididas en zonas de interrogación, la resolución del mapa de vectores

es inversamente proporcional al tamaño de las zonas de interrogación. Por lo tanto, para maxi-

mizar la resolución espacial de los resultados, es necesario minimizar el tamaño de las zonas de

interrogación. Sin embargo, si las zonas son muy pequeñas puede que no se encuentre ninguna

partícula y no se podrá determinar la velocidad en esa región.

3Número de vectores por área

Un nuevo algoritmo en la técnica de velocimetría por imágenes de partículas 3

1.1. Antecedentes.

En los inicios de la técnica de PIV, a principios de los ochenta del siglo pasado, las imágenes

de fluidos eran tomadas con cámaras analógicas para después ser procesadas y estudiadas de forma

manual. Este proceso era muy lento; montar el experimento, tomar un par de fotografías y procesarlas,

podía tardar hasta una semana.

El principal avance de la técnica de PIV sucedió cuando la toma analógica de datos y las técnicas de

evaluación manual se reemplazaron por cámaras digitales y la evaluación realizada por computadora.

Este cambio permitió reducir el tiempo de la adquisición de las imágenes y su procesamiento a segundos.

Existe una gran variedad de fuentes donde se describe el desarrollo histórico de los principios básicos

de PIV. En las referencias [1] y [2] se puede encontrar más información de la evolución que ha tenido

esta técnica.

La Figura 1.2 (a), (b) muestra dos imágenes de partículas tomadas en tiempos consecutivos. Al

aplicar la técnica de PIV a estas imágenes se obtiene el mapa de vectores de velocidad de la figura

(c), en la figura (d) se muestra este mismo mapa dibujado sobre la imagen (a). Cada vector representa

el movimiento de las partículas dentro de una zona de interrogación y generalmente es dibujado en el

centro de cada zona.

Figura 1.2: (a) (b): Imágenes de partículas, (c) (d): Mapa de vectores de velocidad



La técnica de PIV tiene una gran variedad de aplicaciones y es utilizada en todo el mundo para

resolver diversos problemas. Algunos ejemplos de su aplicación son:

en medicina, se ha utilizado para estudiar el flujo sanguíneo en el corazón [3], y el flujo del aire

por las fosas nasales durante la respiración [4], [5],

en agricultura, para diseño de sistemas de ventilación [6],

en acústica, para estudiar la propagación del sonido [7],

en ingeniería, para establecer el efecto del aire en las construcciones [8] y en el diseño aerodinámico

de automóviles y aviones [9],

en reología, para elestudio de flujos complejos [10],

y en biología para conocer la manera en que se desplazan aves [11] y peces [12].

El futuro de la técnica depende en buena medida del desarrollo de las áreas que conforman la

técnica: óptica, tecnología láser, tecnología de cómputo y procesamiento digital de imágenes. A casi

veinticinco años de la creación de PIV, se busca producir equipos más modernos, veloces y económicos,

así como desarrollar algoritmos más rápidos y precisos.

1.2. Planteamiento del problema.

La técnica de PIV tiene una gran proyección e importancia a nivel mundial por ser una de las

técnicas de velocimetría que ofrece más ventajas. Esta técnica permite calcular la velocidad de forma

rápida, indirecta y con una resolución espacial grande, a comparación de otras, que solo obtienen la

velocidad en un punto en particular del fluido a analizar.

Sin embargo la técnica aún tiene varias áreas que deben ser mejoradas para aumentar la calidad

de los resultados. PIV puede dar como resultado una gran cantidad de vectores incorrectos, así como

una baja resolución espacial de resultados. Aunque en muchas aplicaciones la resolución obtenida es

aceptable, existen algunos experimentos donde se necesita más detalle de los resultados. Esto no es

posible ya que PIV determina el desplazamiento promedio de las partículas que se encuentran en un

área, y si esta área es muy pequeña, la técnica falla ya que quizá no se encuentre ninguna partícula. Por

otro lado, hay imágenes en las que por una iluminación deficiente durante el experimento, las partículas

prácticamente no se vean. En este caso la técnica también falla y genera resultados erróneos. Para una


1.3. Objetivos.

interpretación correcta del flujo, no se pueden utilizar este tipo de resultados. Por lo tanto, es necesario

utilizar algoritmos que solucionen los problemas mencionados y permitan obtener una mayor exactitud

y mejor resolución espacial del mapa de vectores de velocidad, con el fin de estudiar apropiadamente

el movimiento del fluido.

1.3. Objetivos.

Objetivo General

El objetivo principal de la presente tesis es diseñar y programar un algoritmo que aumente la

resolución espacial y calidad de los resultados obtenidos con velocimetría por imágenes de partículas.

Objetivos Específicos

Comprender los componentes, funcionamiento, aplicaciones e importancia de la técnica de veloci-

metría por imágenes de partículas.

Proponer un algoritmo que sea eficiente, flexible y que utilice la información que se conozca del

fluido y del experimento. El algoritmo debe mejorar las imágenes antes de su procesamiento y validar

los vectores resultantes.



1.4. Solución propuesta.

La solución que se presenta en esta tesis consiste de un nuevo algoritmo de PIV para determinar el

campo de velocidades de un fluido a partir del procesamiento de sus imágenes de partículas. En este

algoritmo se incorporan dos nuevas etapas al algoritmo básico de PIV (ver Capítulo 2); por lo que se

divide en tres módulos: pre-procesamiento, procesamiento y post-procesamiento.

El pre-procesamiento se basa en la teoría de procesamiento de imágenes, se aplican filtros para

mejorar las imágenes de entrada de PIV. En la fase de procesamiento se aplica el algoritmo básico de

PIV, basado en correlación cruzada de imágenes con la opción de realizar un procesamiento iterativo

y utilizar diferentes resoluciones espaciales de vectores. Finalmente en la etapa de post-procesamiento

se hace uso de un algoritmo de validación y sustitución de vectores resultantes.

El algoritmo fue programado totalmente Matlab y cuenta con una interfaz gráfica para facilitar su

uso. Se restringirán los recursos extras de Matlab, como el uso del lenguaje de programación Java o

C/C++ para que la aplicación solo dependa de la versión requerida de Matlab (ver Apéndice D.2).

El sistema recibe como entrada un par de imágenes de partículas, el tiempo que transcurrió entre la

toma de ellas y la configuración de procesamiento deseada. El resultado es el mapa con los vectores de

velocidad del fluido.

1.5. Aportaciones.

Se diseñó un nuevo algoritmo para determinar el campo de velocidades de un fluido a partir de pares

de imágenes consecutivas de partículas, y se implementó en una aplicación de Matlab. A diferencia de

otros sistemas de PIV en está aplicación se creo una interfaz de usuario gráfica para facilitar el uso e

introducción de datos.

En algoritmo propuesto cuenta con dos etapas adicionales al procesamiento básico de la técnica de

PIV, con el fin de mejorar los resultados obtenidos. La etapa pre-procesamiento está diseñada para

mejorar las imágenes de partículas de PIV para optimizar su posterior evaluación. En la siguiente

etapa, se cuenta con la opción de que el procesamiento sea iterativo utilizando distintas resoluciones de

vectores para mejorar la cantidad y precisión de los resultados. Se implementó un post-procesamiento

que permite validar y sustituir los vectores incorrectos, para así minimizar el error en los resultados.

Finalmente se permite graficar estos resultados y guardar el mapa de vectores resultante en un archivo

para un posterior análisis.


1.6. Estructura de la tesis.

1.6. Estructura de la tesis.

A continuación se describe brevemente el contenido de los capítulos que conforman la tesis.

Capítulo 1 Introducción. En este capítulo se da una introducción a la técnica de PIV, se describe

el problema de la técnica y se plantea una solución al mismo, así como los objetivos de este trabajo.

Capítulo 2 Descripción de los componentes de la técnica de PIV. Se describen los componentes

físicos (partículas, láser, computadora) y su función en la técnica de PIV. Además, se explica el pro-

cesamiento para determinar el campo de velocidades de un fluido a partir de sus imágenes.

Capítulo 3 Algoritmo propuesto y su implantación. Aquí se detallan las etapas del algoritmo pro-

puesto, así como su función e importancia. Se muestran distintos ejemplos de la evolución que tiene

una imagen desde que es tomada por la cámara hasta que se genera el mapa de vectores de velocidad.

Capítulo 4 Pruebas y resultados. En este capítulo se muestra una comparación de los resultados de

las pruebas realizadas a distintas imágenes de PIV (tanto imágenes artificiales como experimentales).

Las pruebas se realizaron con el programa realizado para esta tesis así como con los programas Urapiv,

Matpiv y Mpiv.

Capítulo 5 Conclusiones. Se presentan las conclusiones del trabajo realizado en esta tesis, así como

las posibles mejoras a realizar en un trabajo futuro.

Apéndice A Características de los componentes de PIV. Se dan más especificaciones de componentes

como las partículas, el láser y la cámara.

Apéndice B Diagramas. Se muestran los diagramas de flujo y de estados de las diferentes etapas

del sistema.

Apéndice C Código. Se muestra la codificación de la interfaz y rutinas elaboradas para este sistema.

Apéndice D Instalación del sistema. Se detalla toda la información para la instalación del sistema

en cualquier equipo de cómputo.

Apéndice E Manual de Usuario. Se presenta el manual de usuario para la implementación del

sistema propuesto.

Apéndice F Pruebas Realizadas. Se muestran algunas pruebas realizadas con otras imágenes de

PIV, asimismo se presentan comparaciones de resultados con otros sistemas de PIV gratuitos.


CAPÍTULO 2

La Técnica de Velocimetría por Imágenes de Partículas

2.1. Partículas trazadoras.

La técnica de PIV es utilizada para determinar la velocidad de un fluido en forma indirecta.

La velocidad es calculada mediante el procesamiento de imágenes del fluido con partículas, que con

anterioridad se insertan en el fluido. Las partículas son llamadas “partículas trazadoras” y por sus

características, su movimiento será muy parecido al del fluido. Por lo tanto la velocidad de una partícula

será prácticamente la misma que la velocidad del fluido en cada punto.

Las partículas empleadas en PIV deben de satisfacer dos condiciones: ser capaces de seguir el

movimiento del fluido sin retraso excesivo, y reflejar la suficiente luz del láser para ser capturadas

adecuadamente por la cámara. Lo ideal es encontrar partículas que al ser insertadas en el fluido, sigan

el movimiento del fluido de forma exacta, sin modificar el flujo, ni interactuar entre ellas. Sin embargo,

es difícil satisfacer estas condiciones y en la práctica se utilizan aquellas que simplemente sigan el flujo

sin perturbarlo demasiado. Las partículas y el fluido deben tener aproximadamente la misma densidad,

para que al insertar las partículas, éstas permanezcan suspendidas en el fluido. Si las densidades de las

partículas y del fluido son diferentes, las partículas flotarán o se hundirán y no seguirán apropiadamente

el movimiento del fluido.

También es necesario que las partículas dispersen bien la luz del láser, pues de lo contrario éstas no

serían captadas apropiadamente por la cámara. Sin embargo la mejor y más económica opción, es elegir

adecuadamente las partículas que utilizar un láser con mayor potencia. En general la luz dispersada por

las partículas está en función de su índice de refracción y del índice del fluido que las rodea, así como

del tamaño y forma de las partículas. La luz dispersada también depende del ángulo de observación,

el cual generalmente es a 90o de la luz incidente, por ser una de las posiciones con mayor dispersión

de luz. En la referencia [1] se encuentra una descripción más detallada de la teoría de dispersión y los

2.2. Fuente de luz.

ángulos de observación.

Todas las partículas deben ser prácticamente del mismo tamaño, y su elección es un compromiso

entre los beneficios y desventajas, según el tamaño y cantidad utilizados. Las partículas grandes son

las que reflejan mejor la luz. Sin embargo las partículas pequeñas siguen mejor el movimiento del fluido

y distorsionan menos el flujo. La cantidad de partículas utilizada también es muy importante. Tener

muchas partículas o muy pocas afectará los resultados obtenidos. Si se utilizan muchas partículas, se

modifica el flujo y si se utilizan pocas partículas, no se obtiene suficiente información para calcular la

velocidad en todos los puntos del área de estudio. La distribución de las partículas debe ser homogénea,

para poder obtener la velocidad en todas las regiones del fluido. La elección del tipo de partículas en

cada experimento depende del fluido y del tipo de flujo que se desee estudiar. En el Apéndice A se

encuentran unas tablas con las partículas que pueden ser utilizadas para visualizar líquidos y gases.

2.2. Fuente de luz.

La técnica PIV requiere la generación de un plano de luz con características bien definidas de área,

grosor, color e intensidad. Los rayos láser son muy utilizados en la técnica de PIV por ser el único tipo

de luz que, por medio de lentes, puede concentrarse para formar un plano de luz. Este plano se hace

incidir sobre el fluido para iluminar de manera casi uniforme la región de interés. Otra razón por la

que generalmente se utiliza un láser y no otra fuente de luz, es la posibilidad de utilizar los dispositivos

que generan el rayo en pulsos. Este tipo de dispositivos generan los pulsos a diferentes frecuencias y

su duración varía según el modelo. El láser emitido por pulsos es sincronizado, usualmente por medio

de una computadora u otro dispositivo de hardware, con la cámara, de modo que el obturador de la

cámara solo permanezca abierto por un intervalo de tiempo muy corto, mientras se tenga iluminada la

región de interés del fluido. Así, la cámara solo capturará la posición de las partículas en un instante

y el movimiento de las partículas se “congelará”, evitando que se pueda percibir el desenfoque de

movimiento en las imágenes (motionblur).

Las dimensiones del volumen iluminado por el láser son controladas por un conjunto de lentes. Los

sistemas láser comerciales traen integrada la óptica en su interior, y la salida del dispositivo ya es un

plano de luz. La mayoría de estos sistemas permiten ajustar:

El grosor y ángulo del plano de luz.

La orientación del plano en 360o.


Capítulo 2. La Técnica de Velocimetría por Imágenes de Partículas

La distancia de enfoque.

Duración de cada pulso.

Intervalo entre los pulsos

Existen diferentes tipos de láser que cumplen con los requerimientos de la técnica de PIV. La

elección del láser depende de la velocidad del fluido que se quiere estudiar y el tamaño de la región a

iluminar. El tamaño de la región a iluminar en el fluido es importante, pues no todos los dispositivos

tienen la misma potencia para garantizar que el plano generado sea del tamaño del área a estudiar,

así como procurar que la iluminación del área sea uniforme. Además, la fuente de luz debe de ser

suficientemente potente para que, sin importar su tamaño, las partículas sean iluminadas de modo tal,

que reflejen la cantidad de luz necesaria para que sean visibles y puedan ser capturadas por el sensor de

la cámara. El tipo de láser más utilizado en PIV es el de Neodimio:GAI (Granate Alumínico de Itrio)

ó Nd:YAG por sus siglas en inglés. El láser Nd:YAG es de estado sólido y es capaz de seguir flujos muy

rápidos al producir luz muy intensa en pulsos cortos (pulsos de nanosegundos) a una gran frecuencia.

Además, generalmente los equipos láser de este tipo son compactos. El Apéndice A contiene una tabla

con algunas propiedades de los sistemas Nd:YAG PIV láser comerciales, así como una descripción de

otros tipos de láser que son utilizados en PIV.

Actualmente existen varias empresas dedicadas exclusivamente a perfeccionar los sistemas láser

utilizados en PIV (como Oxford Lasers [13] y Big Sky Laser Technologies [14]). El objetivo es que los

nuevos equipos sean cada vez más compactos, fáciles de operar y que se adecuen a las necesidades de

una gran variedad de experimentos.

2.3. Captura y almacenamiento de imágenes.

Es necesario tener un sistema de captura que permita obtener imágenes del fluido con partículas.

Actualmente el medio más utilizado es la cámara digital con CCD (Charge-Coupled Device o Dispo-

sitivo Acoplado por Carga), aunque estas imágenes pueden ser capturadas por otro tipo de cámaras.

El CCD está formado por millones de sensores fotosensibles acomodados en una cuadrícula. Cada sen-

sor del CCD representará un píxel1 en la imagen final. La resolución de una cámara es el número de

sensores que tiene el CCD y por lo tanto el número de píxeles de la cual esta formada la imagen.

En PIV la resolución es muy importante, pues cuando las partículas son muy pequeñas solo pueden1Píxel viene de Picture element, que es la parte más pequeña de la que esta formada una imagen


2.3. Captura y almacenamiento de imágenes.

ser capturadas por cámaras de alta resolución. Sin embargo una imagen con mayor resolución ocupa

mas espacio en la memoria y también involucra más tiempo de procesamiento. Las cámaras digitales

con CCD convencionales necesitan pasar los datos del CCD a memoria y después borrar los valores de

los sensores, antes de poder tomar una nueva imagen. Existen unas cámaras especiales, llamadas de

“correlación cruzada”, que son capaces de tomar dos imágenes en el mismo CCD sin tener que borrar

ni guardar los datos entre una imagen y otra. En PIV, estas cámaras son útiles para analizar fluidos

rápidos, ya que pueden tomar imágenes con una frecuencia muy alta. En el Apéndice A se encuentra

una descripción un poco más detallada del CCD, así como las principales características de las cámaras

utilizadas en experimentos de PIV.

El intervalo de tiempo que debe haber entre cada par de imágenes depende de la velocidad del fluido.

Si el fluido es muy rápido y el intervalo entre las imágenes es muy grande, dos imágenes consecutivas

no tendrán ninguna partícula en común y será imposible determinar correctamente la velocidad del

fluido. Por otro lado, el intervalo de tiempo entre imágenes consecutivas debe ser lo suficientemente

largo para que las partículas tengan un desplazamiento.

Existen diferentes formas de registro de imágenes de PIV y éstas pueden dividirse en dos gran-

des grupos: métodos que generan una sola imagen o los métodos que generan una imagen por cada

exposición. Estos grupos son conocidos como: PIV de una imagen con múltiples exposiciones (single

frame/multi-exposure PIV) y PIV de múltiples imágenes con una exposición (multi-frame/single ex-

posure PIV) respectivamente. En las siguientes figuras se muestran estos dos tipos de registro. En la

Figura 2.1 se muestra el caso de una imagen con una partícula con diferentes exposiciones. En la figura

(a) el obturador estuvo abierto durante un prolongado intervalo de tiempo, por lo tanto la trayectoria

de la partícula se registró formando una línea. En las figuras (b) y (c) se muestran imágenes sobreex-

puestas, donde el obturador estuvo abierto varias veces en cortos periodos de tiempo y solo se registro

la posición de la partícula durante cada una de esas exposiciones. Cada exposición fue realizada en el

tiempo tk = ∆t •k, para k = 0, 1 y 2. Como se puede ver en estas últimas imágenes también es posiblever la trayectoria de la partícula sin embargo el problema es que existe ambigüedad entre la dirección

de la trayectoria.

En la Figura 2.2 se muestran múltiples imágenes con diferentes exposiciones. A diferencia de las

imágenes anteriores, en este tipo de registro no se generan imágenes sobreexpuestas. Por cada exposición

se genera una imagen y al tener múltiples imágenes se tiene de forma inherente el orden de las imágenes

de partículas y su evaluación es más sencilla. En cada imagen se muestra la posición de la partícula

durante una exposición. Como referencia se muestra la posición de la partícula en las tomas anteriores

de cada imagen.



Figura 2.1: (a) Imagen con una exposición, (b) Imagen con doble exposición, (c) Imagen con múltiple

exposición.

Figura 2.2: (a) Dos imágenes con una exposición, (b) Múltiples imágenes con una exposición.


2.4. Evaluación de las imágenes.

En esta tesis solo se considera la evaluación de pares de imágenes con exposición simple, ya que

actualmente éste es el tipo de registro más utilizado.

2.4. Evaluación de las imágenes.

Las imágenes empleadas en la técnica de PIV son capturadas por la cámara digital y son el resultado

de la exposición de las partículas iluminadas. Estas imágenes son la representación bidimensional del

volumen del fluido iluminado por el plano láser (ver Figura 2.3).

Figura 2.3: Representación de las imágenes de PIV.

Una imagen es una función bidimensional de intensidad luminosa y se denota por f(x, y), donde

el valor de f en las coordenadas (x, y) da la brillantez de la imagen en ese punto. En PIV se utilizan

imágenes digitales creadas por una cámara digital y están formadas por un conjunto de valores digitales

llamados píxeles. Los píxeles están espaciados uniformemente y son discretos. En la Figura 2.4 se

muestra una matriz que es la representación de una imagen digital f(x, y) de tamaño N x M. Los

píxeles son almacenados en la memoria de la computadora en un arreglo bidimensional cuyos índices

de fila y columna indican la posición de cada píxel en la imagen. Los valores de este arreglo son

enteros positivos y representan la intensidad o nivel de gris de los píxeles de la imagen. La intensidad

de una imagen f en las coordenadas (x, y) es el nivel de gris f(x, y). El intervalo de valores que

puede tener f(x, y) es [0, L]. Donde f(x, y)= 0 se considera negro y f(x, y) = L es el blanco en la

escala y generalmente es 255. Todos los valores intermedios son tonos de gris que están igualmente



espaciados y varían del negro al blanco. El procesamiento de las imágenes de partículas es realizado

f(x, y) ≈

⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝

f(0,0) f(0,1)

f(1,0) f(1,1)

. . .f(0,M−1)

f(1,M−1)

.... . .

...

f(N−1,0) f(N−1,1) · · · f(N−1,M−1)

⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠Figura 2.4: Representación digital de una imagen

por una computadora de forma digital. Este procesamiento consiste en manipular numéricamente las

imágenes digitales para mejorarlas, transformarlas y obtener información útil de ellas, como el campo

de velocidades del fluido. El algoritmo de la técnica de PIV en su forma más básica consta de una

etapa de procesamiento en la cual se determina el campo de velocidades. Sin embargo, es importante

incluir dos etapas adicionales: pre-procesamiento y post-procesamiento. En las siguientes secciones se

describen brevemente estas etapas.

En la referencia [15] se puede encontrar la teoría del procesamiento digital de imágenes. En los

capítulos 3 y 4, de esta referencia, se encuentran varios filtros y transformaciones útiles para mejorar

las imágenes de PIV en la etapa del pre-procesamiento; así como el uso de la transformada rápida de

Fourier (FFT) para calcular la correlación cruzada de dos imágenes, la cual es utilizada en la etapa de

procesamiento.

2.5. Procesamiento

Una vez adquiridas y mejoradas las imágenes, éstas son procesadas para determinar la velocidad

del fluido. El cálculo de la velocidad de un fluido mediante las técnicas de velocimetría por imágenes

de partículas se basa en la relación:

V elocidad = KDesplazamiento

T iempo.

Para calcular la velocidad es necesario determinar el desplazamiento que tuvo el fluido en un intervalo

de tiempo y conocer el factor de escala K. El intervalo es un dato conocido pues es parte de la

configuración inicial del experimento y el factor de escala Kes la constante escalar que da la relación

entre píxeles en la imagen y unidades de longitud en el experimento. Por lo tanto, solo es necesario

determinar el desplazamiento del fluido en la dirección x e y, para calcular la velocidad del fluido en


2.5. Procesamiento

estas dos direcciones. Estos desplazamientos se calculan mediante las partículas de un par de imágenes

consecutivas. El método de procesamiento depende de la densidad de partículas en la imagen. Se pueden

distinguir tres tipos de densidad de partículas en imágenes: baja, mediana y alta densidad; como se

muestra en la Figura 2.5 la densidad de las partículas en la imagen determina el método de análisis y

procesamiento que debe ser utilizado para determinar la velocidad.

Figura 2.5: Los tres tipos de densidad de partículas en las imágnes: (a) Baja (PTV), (b) Media (PIV),

y (c) Alta (LSV).

La densidad de partículas es pequeña cuando el promedio de separación entre ellas es considera-

blemente mayor a su desplazamiento. Por el contrario, la densidad de partículas es alta, cuando el

promedio de separación entre las partículas es menor a su desplazamiento.

En el caso de la imagen con baja densidad de partículas, Figura 2.5 (a), las partículas pueden ser

detectadas de forma individual. Para la evaluación de este tipo de imágenes, se requieren métodos de

seguimiento de partículas para determinar el desplazamiento individual de cada una. El desplazamiento

de una partícula determina la velocidad del fluido en ese punto. Sin embargo, por la baja densidad de

partículas, la información de la velocidad solo se puede obtener cuando alguna partícula está presente.

Así, pueden existir regiones del flujo donde no se puede obtener una velocidad local. El método para

analizar un fluido en movimiento con estas características se llama “Particle Tracking Velocimetry”

o “PTV” (Velocimetría por seguimiento de partículas). PTV tiene un porcentaje de error alto en

imágenes con ruido. Este problema se debe a que el ruido es confundido con partículas durante el

procesamiento de la imagen.

Las imágenes con densidad media, Figura 2.5 (b), tienen partículas que también pueden ser distin-

guidas de forma individual, sin embargo a simple vista ya no es posible distinguir una misma partícula

en las dos imágenes. Este tipo de imágenes son procesadas con la técnica convencional de PIV, la

cual hace una evaluación estadística y determina el desplazamiento promedio de pequeños grupos de



partículas. Al ser un promedio del desplazamiento la técnica de PIV es más tolerante al ruido que

PTV.

En el caso de imágenes de alta densidad, Figura 2.5 (c), son tantas las partículas, que es imposible

identificar partículas de forma individual porque éstas están encimadas en la mayoría de los casos y se

forman manchas o también llamadas “speckles”. El método utilizado para éste caso es llamada “Láser

Speckle Velocimetry” o “LSV” (Velocimetría de Manchas por Láser). LSV hace una aproximación es-

tadística, que revela el desplazamiento de las manchas formadas por partículas. La principal desventaja

que tiene este método es que, por la alta densidad de partículas, éstas tienden a afectar la dinámica

original del fluido.

La Figura 2.6 (a) y (b) muestra dos imágenes consecutivas con baja densidad de partículas. Este tipo

de imágenes son procesadas para determinar el desplazamiento individual que tuvieron las partículas

entre una imagen y otra (PTV). La Figura 2.6 (c) muestra este desplazamiento. Cada una de estas

partículas puede ser vista como una pareja de partículas (una partícula del tiempo t0 y otra del tiempo

t1).

Figura 2.6: (a) Imagen de partículas tomada en el tiempo t0, (b) Imagen de partículas tomada en el

tiempo t1, (c) Desplazamiento de las partículas en el intervalo t0-t1.

Las imágenes de PIV, Figura 2.7, tienen una densidad de partículas mayor que las mostradas en

la Figura 2.6 y no es posible determinar el desplazamiento individual de las partículas. Por esta razón

se utilizan métodos estadísticos para determinar el desplazamiento promedio de pequeños grupos de

partículas.

El método de evaluación para PIV más usado es la correlación cruzada de las imágenes. En el

procesamiento de imágenes, la correlación cruzada es utilizada para encontrar la posición de mayor

semejanza entre dos imágenes. Para la evaluación con correlación cruzada, las imágenes de PIV son

divididas en pequeñas áreas rectangulares llamadas “zonas de interrogación”. Una zona de interrogación


2.5. Procesamiento

Figura 2.7: Imagen de partículas de PIV con la amplificación de una parte de la imagen.

es una región de una imagen f , es decir es una subimagen de f .

Se asume que todas las partículas que están dentro de una de éstas zonas, se mueven de manera

uniforme. El algoritmo consiste en procesar dos zonas de interrogación del mismo tamaño y con las

mismas coordenadas dentro de la imagen. Una zona es de la primera imagen y la otra de la segunda

imagen como se muestra en la Figura 2.8.

Figura 2.8: Dos imágenes de PIV divididas en zonas de interrogación.

Las zonas de interrogación deben ser procesadas con el fin de encontrar el desplazamiento promedio

de las partículas en las imágenes. Existen varias técnicas para procesar estas zonas de interrogación

[16]. La más utilizada actualmente consiste en calcular directamente la correlación cruzada de las dos



zonas:

f(x, y)¯ g(x, y) =MX

i=−M

NXj=−N

f(i, j)g(i+ x, j + y) (2.1)

para x = 0,±1,±2...±N−1 y = 0,±1,±2...±M−1, y donde f(x, y) es la subimagen correspondientea la primera zona de interrogación, g(x, y) es la segunda zona, M es el ancho y N la altura en píxeles

de la zona de interrogación. El proceso consiste en desplazar una subimagen sobre la otra y sumar

los productos de los valores en donde haya traslape. Cada valor es guardado en una matriz llamada

“matriz o plano de correlación cruzada”, en la posición que represente el desplazamiento de la segunda

subimagen sobre la primera.

En la Figura 2.9 se muestra como se desplazan las zonas de interrogación. La primera zona es el

cuadrado del centro y la segunda zona es el cuadrado gris que se desplaza. Las zonas de interrogación

son de 4 por 4 píxeles (M por N) y el desplazamiento hacia la derecha y arriba es positivo. El cuadrado

negro es la región en la existe un traslape entre las dos zonas, por ejemplo, en la imagen (b) no existe

desplazamiento y el traslape es la imagen completa. La imagen (e) muestra un desplazamiento positivo

de 2 en x y uno negativo de 2 en y. La suma del producto de este traslape es guardada en el plano

de correlación cruzada en la posición (2,−2). Así se puede ver como se construye este plano que es unarreglo bidimensional con los valores de correlación cruzada y cuyo tamaño es 2M−1 por 2N−1. Estetamaño es porque aunque los desplazamientos se podrían extender hasta ± infinito, generalmente en

la práctica solo se consideran los desplazamientos mientras exista traslape.

Hacer la correlación cruzada directa de dos imágenes de esta forma es más complejo computacio-

nalmente y llega a ser prohibitivo en imágenes grandes. El empleo de la transformada de Fourier es una

alternativa para calcular la correlación cruzada con una complejidad inferior. El Teorema de correla-

ción cruzada establece una equivalencia entre la función (2.1) y la multiplicación de las transformadas

de Fourier de las imágenes. Así, la correlación cruzada puede ser calculada de la siguiente forma:

f(x, y)¯ g(x, y) ⇔ F ∗(u, v)G(u, v) (2.2)

donde

F ∗(u, v) = Re(F (u, v))− iIm(F (u, v)) (2.3)

En estas fórmulas f(x, y) es la primera imagen, g(x, y) es la segunda imagen, F ∗(u, v) es el conjugado

de la transformada de Fourier de la primera imagen y G(u, v) es la transformada de Fourier de la

segunda imagen.


2.5. Procesamiento

Figura 2.9: Correlación de dos zonas de interrogación.

En el plano de correlación cruzada se pueden observar las posiciones de los máximos en los cuales

hubo mayor similitud al traslapar una imagen sobre la otra. Generalmente el máximo absoluto repre-

senta el desplazamiento promedio que tuvieron las partículas entre las dos imágenes. Las coordenadas

de este máximo representan el desplazamiento en píxeles que tuvieron las partículas en cada eje. La

Figura 2.10 (c) y (d) muestra de forma gráfica el plano de correlación cruzada de las imágenes (a) y

(b). Las imágenes tienen un tamaño de 64 por 64 píxeles y por lo tanto, el plano de correlación cruzada

es de 127 por 127 píxeles y fue calculado con la fórmula (2.1). En este ejemplo se puede ver que en

el plano de correlación cruzada existen varios máximos. Estos seis máximos locales se deben a que al

traslapar una imagen sobre la otra hay un empalme de un círculo blanco de la primera imagen con

otro círculo de la segunda. El máximo absoluto se encuentra en la posición (36,13) y corresponde al

empalme de los tres círculos de ambas imágenes.

La Transformada Rápida de Fourier asume que la imagen de entrada es periódica y por lo tanto

se copia o repite la información de los extremos contrarios cuando se realizan los cálculos. Al hacer

esto la salida presenta información repetida (máximos con la misma magnitud, pero en diferentes

posiciones). Por esto es necesario aplicar filtros que permitan eliminar esta información, para que

permanezcan solos los datos importantes que representen el desplazamiento real de las partículas.



Figura 2.10: (a) y (b) Imágenes de partículas, (c) y (d) Plano de correlación cruzada.


2.5. Procesamiento

Estos filtros consisten básicamente en aplicar una máscara para tomar solo la información del centro

de cada zona de interrogación. La máscara puede ser generada por una función Gaussiana o totalmente

cuadrada.

El proceso de correlación cruzada mencionado se repite para todas las zonas de interrogación de

las imágenes y como resultado se obtiene una matriz con los desplazamientos en las coordenadas x e

y de cada zona de interrogación. Con la matriz de desplazamientos, el intervalo de tiempo entre las

imágenes capturadas y el factor de amplificación de la imagen se calcula la velocidad local del fluido

para cada zona. Finalmente la velocidad de cada zona se representa con un vector, y se construye un

mapa de vectores de velocidad. Típicamente las zonas de interrogación son cuadradas, por lo tanto el

mapa de velocidades obtenido presenta los vectores acomodados en una matriz como se muestra en la

Figura 2.11.

El tamaño de las zonas de interrogación determina el número de vectores que son calculados por

imagen y por lo tanto determinará la máxima resolución espacial del mapa de vectores de velocidad.

Entre más pequeña sea la zona de interrogación se obtienen más vectores. Por ejemplo en la Figura 2.12

se pueden ver dos imágenes (a) y (b), y dos mapas de vectores con diferentes resoluciones dependiendo

del tamaño de zona de interrogación utilizada. La Figura 2.12 (c) tiene zonas de 256 por 256 píxeles y

genera 4 vectores, mientras que la Figura 2.12 (d) con zonas de 32 por 32 píxeles genera 256 vectores.

La técnica de PIV permite obtener el campo de velocidades a partir de las imágenes de partículas

de un fluido en movimiento. A continuación se mostraran algunos de los mapas que se pueden obtener

al utilizar PIV en un par de imágenes consecutivas. En las secciones previas se mencionaron las partes

y etapas necesarias para obtener las imágenes. Ahora estas imágenes son procesadas para obtener los

mapas de vectores de velocidad y su análisis e interpretación dependerá de lo que se este estudiando

en el fluido. Además, en estos ejemplos se podrá ver cuál es la importancia de elegir adecuadamente el

tamaño de las zonas de interrogación y su repercusión en la cantidad y validez de los resultados.

La Figura 2.13 muestra dos imágenes de partículas generadas de forma artificial [17] de 256 por 256

píxeles. Cada imagen tiene aproximadamente 4000 partículas de 5 píxeles de diámetro en promedio.

El desplazamiento medio de las partículas es de 7.5 píxeles entre el par de imágenes.

A continuación se muestran algunos mapas de vectores que son resultado de aplicar el algoritmo

de PIV a las imágenes de la Figura 2.13. Las imágenes fueron correlacionadas con distintas zonas de

interrogación, aunque para fines prácticos, todas fueron cuadradas y de tamaño potencia de 2. Cabe

mencionar que la técnica no impone ninguna restricción en cuanto al tamaño o la forma de las zonas

de interrogación, por lo cual éstas podrían ser rectangulares.



Figura 2.11: Proceso para obtener el mapa de vectores de desplazamiento.

Figura 2.12: (a) (b) Imágenes de partículas 512 por 512 píxeles, (c) mapa de vectores con zonas de 256

por 256 píxeles (d) mapa de vectores con zonas de 32 por 32 píxeles.


2.5. Procesamiento

Figura 2.13: (a) Imagen de partículas tomada en el tiempo t0 , (b) Imagen de partículas tomada en el

tiempo t1.

Figura 2.14: Mapa de vectores con zonas de (a) 128 x 128 píxeles (b) 64 x 64 píxeles.



Aunque va a depender de la aplicación que se le de a los datos, generalmente es más conveniente

tener una mayor resolución espacial de los resultados. Como se mencionó previamente, la resolución

(cantidad de resultados) depende del tamaño de las zonas de interrogación utilizadas. En la Figura

2.14 (a) con zonas de 128 por 128 solo se obtienen 4 vectores mientras que en la Figura 2.16 (b) se

obtienen poco más de 4000.

Es conveniente utilizar zonas de interrogación pequeñas para tener una mayor resolución y de-

talle del fluido que se esta estudiando, sin embargo si las zonas son muy pequeñas la validez de los

resultados puede verse afectada. En la Figura 2.15 (b) se pueden ver algunos vectores que difieren

considerablemente en magnitud o dirección de sus vecinos. Estos vectores son llamados incorrectos o

espurios.

Tener vectores incorrectos es inevitable, y aun en condiciones ideales un 5% de los vectores será

incorrecto. Los vectores incorrectos se deben a que al correlacionar las zonas de interrogación, el máximo

absoluto no se produce en el punto de mayor similitud entre ellas. Las causas más comunes son porque

que el fluido o el medio que lo contiene producen reflejos de la luz del láser que son confundidos con

partículas; o porque si las zonas son muy pequeñas puede que se encuentren muy pocas partículas o

ninguna, lo cual haría imposible encontrar un mapeo con las partículas de la segunda imagen. En las

orillas de las imágenes es más común tener vectores incorrectos. Debido a que el fluido se desplaza

entre una imagen y la otra, es normal que algunas partículas que se encuentran en una orilla de la

primera imagen ya no se encuentren en la segunda imagen. Para corregir el problema de los vectores

incorrectos muchos de los sistemas comerciales de PIV tienen rutinas para validar los resultados de los

mapas de vectores y sustituirlos en caso de ser necesario. Sin embargo son pocas las opciones gratuitas

que cuentan con esta etapa de post-procesamiento.

Al utilizar la técnica de PIV, se tiene la limitante de la resolución espacial de sus resultados

pues siempre habrá un momento a partir del cual si se sigue reduciendo el tamaño de las zonas

de interrogación, el número de vectores incorrectos superará el número de vectores correctos. Por esta

razón, es muy importante elegir adecuadamente el tamaño de las zonas de interrogación. Como ejemplo

de esta situación tenemos la Figura 2.16 (b) en la cual la mayoría de los vectores son incorrectos. Un

número que se puede tomar como referencia es buscar que haya por lo menos entre 5 y 10 partículas

por zona de interrogación [18]. También se puede tomar como referencia la regla del cuarto, que dice

que el tamaño de la zona de interrogación debe de ser cuatro veces el desplazamiento promedio de las

partículas.

El principal problema de PIV es que la resolución espacial de los resultados no es suficiente en

muchos experimentos. El tamaño de las zonas de interrogación con el que aún se obtienen buenos


2.5. Procesamiento





resultados generalmente es de 64 por 64 o de 32 por 32 píxeles. Si se desea obtener una mayor resolución

es necesario utilizar otro tipo de procesamiento, como los propuestos en Súper-Resolución y Multigrid.

La téncica de Súper-Resolución es un método híbrido de PIV/PTV en el cual la correlación cru-

zada de PIV es empleada antes de la identificación de las trayectorias individuales de las partículas.

Las tareas críticas de este método son identificar correctamente las partículas de forma individual y

determinar con precisión el centro de las partículas identificadas. Esto debido a que las condiciones

de las imágenes de PIV hacen difícil la identificación individual. Los algoritmos que combinan PIV y

PTV han sido estudiados por varios autores, enfocándose por la forma en que se detectan las partícu-

las y en cómo realizar la asociación correcta de las partículas en imágenes consecutivas. Una de las

últimas propuestas, y más utilizada, ha sido emplear reglas de lógica difusa para determinar los pares

de partículas [19].

La técnica de Multigrid es un método que permite obtener mapas de velocidades con diferentes

resoluciones. El procesamiento con correlación cruzada se realiza de forma iterativa y en cada iteración

se reduce en un determinado porcentaje el tamaño de las zonas de interrogación para obtener una

mayor resolución de los resultados. Generalmente se hacen tres o menos iteraciones y en cada una se

utiliza la información de la iteración anterior [20].


CAPÍTULO 3

Diseño e Implementación

En el capítulo anterior se presentó la técnica básica de PIV. Esta técnica presenta buenos resultados

si las imágenes tomadas son de muy buena calidad es decir, con iluminación uniforme, con una distri-

bución homogénea de partículas y con un pequeño desplazamiento de las mismas. Sin embargo, en la

práctica algunos de estos factores pueden fallar, produciendo resultados incorrectos. Además de todo

esto, hay que tomar en cuenta la limitante de PIV en cuanto a la resolución espacial de sus resultados.

En este capítulo se describe un nuevo algoritmo de PIV, con el fin de mejorar los resultados

obtenidos con el algoritmo básico descrito en el capítulo anterior. En este nuevo algoritmo se incluyen

dos etapas adicionales de procesamiento. La primera dedicada al mejoramiento y acondicionamiento de

las imágenes de partículas antes de ser procesadas. En esta etapa se mejora la iluminación y contraste

de las partículas en la imagen, mediante el procesamiento digital de imágenes y la aplicación de filtros.

La segunda etapa realiza la corrección de errores de los resultados obtenidos por el algoritmo. El

principal objetivo es eliminar los vectores incorrectos que difieren mucho en dirección o magnitud en

comparación a sus vecinos adyacentes. Después se sustituyen estos vectores por nuevos, calculados a

partir de los vecinos de cada vector incorrecto.

Adicionalmente, el procesamiento de correlación cruzada se realiza de forma iterativa y con el

máximo absoluto del plano de correlación se realiza una interpolación para determinar el desplaza-

miento de las partículas. Esto mejora la resolución de resultados al obtener de forma más precisa el

desplazamiento de las partículas en cada zona de interrogación.

Por lo tanto, el algoritmo propuesto de PIV esta conformado por tres etapas:

Pre-Procesamiento.

Procesamiento.

3.1. Etapa de pre-procesamiento.

Post-Procesamiento.

En la Figura 3.1 se muestra el diagrama de estados del sistema que a continuación se describen en

este capítulo. En el Apéndice B se muestran los algoritmos en diagramas de flujo y de estado de cada

una de las etapas.

Figura 3.1: Diagrama de estados de datos del sistema.

El sistema completo fue implementado en Matlab. En el Apéndice D se encuentran los requeri-

mientos y especificaciones del programa, además de una lista de las funciones utilizadas que forman

parte de los toolbox incluidos en el paquete de Matlab.


El propósito de esta etapa es mejorar la apariencia de las imágenes capturadas, antes de realizar

el procesamiento que determina la velocidad del fluido. En esta etapa se realiza el mejoramiento

de las imágenes en caso de ser necesario corregir ciertos defectos experimentales. La etapa de pre-

procesamiento a su vez esta dividida en 3 etapas: modificación del tamaño de las imágenes, ajuste

del contraste y eliminación de objetos no útiles para el procesamiento de PIV (objetos fijos o con

movimiento en una secuencia de imágenes).


Capítulo 3. Diseño e Implementación

La selección del área de la imagen a trabajar es importante pues en algunos casos no es posible

ajustar eso previamente desde la cámara. Aunque esta etapa es opcional y en teoría debería dar los

mismos resultados en la parte región donde están presentes las partículas. Sin embargo, los resultados

pueden verse afectados debido a que la división de las zonas de interrogación puede no estar en las

mismas posiciones. Además de que no tiene caso procesar regiones de la imagen en las cuales no se

desea o puede estudiar la velocidad del fluido. En las imágenes del laboratorio de Reología, Figura

3.2, solo menos de una sexta parte de la imagen es importante. Procesar la imagen completa tardaría

varias veces más en comparación de sólo procesar el área de interés. En la Figura 3.3 se muestra la

imagen ya recortada. En adelante, todo el procesamiento se realiza sobre la imagen recortada.

Figura 3.2: Imagen original con áreas a recortar.

Uno de los problemas más comunes en las imágenes es la iluminación no uniforme. En ocasiones las

imágenes de partículas que están del lado de la fuente de luz son ligeramente más brillantes que aquellas



Figura 3.3: Imagen recortada.

que están del lado opuesto (ver Figura 3.4). Aunque esto no es un problema grave, muchos programas

comerciales cuentan con rutinas para mejorar las imágenes y poder garantizar que las imágenes de

PIV están iluminadas homogéneamente antes de realizar su procesamiento. Entre las cosas que se

pueden realizar en el pre-procesamiento están los cambios de contraste, brillo, color, resaltar bordes,

eliminación de ruido, ampliar o reducir la imagen etc. En las imágenes de partículas las modificaciones

más comunes son las de tamaño y contraste para resaltar las partículas del fondo.

Figura 3.4: Imagen con objetos no deseados estáticos.

El contraste es la diferencia de intensidades de iluminación de una imagen, lo que hace que un objeto

sea distinguible del fondo. En la adquisición de las imágenes algunas situaciones como una iluminación

deficiente o una mala configuración en la apertura de los lentes, pueden generar una imagen con bajo

contraste. Así un aumento en el contraste de una imagen de PIV, dará como resultado una imagen con

las partículas más definidas, en el caso las partículas blancas serán más distinguibles del fondo negro,

además de distinguir partículas que antes no eran visibles debido a la iluminación no uniforme.



Existe una mejora de contraste llamado estiramiento de contraste [15]. Esta mejora consiste en

incrementar el rango dinámico de los niveles de gris de la imagen a procesar, es decir, que se usen

más valores de la gamma de la escala de grises. Debido al procesamiento, almacenamiento y a las

consideraciones de muestreo de hardware durante el proceso de digitalización de la imagen, los niveles

de gris típicamente son un entero potencia de 2 (L = 2k, para k = 1, 2. . . ).

Se asume que los niveles discretos están igualmente espaciados y que son enteros en el intervalo [0,

L-1]. El rango de los valores a lo largo de la escala de grises es llamado rango dinámico de la imagen.

Las imágenes que ocupan una porción significativa de la escala de grises, son imágenes con un alto

rango dinámico. Cuando un gran número de píxeles tienen esta propiedad, la imagen tendrá un buen

contraste. Por el contrario si una imagen tiene un rango dinámico bajo, tendrá un aspecto grisáceo.

Para incrementar el contraste, el estiramiento fue calculado con la función imadjust. Se obscurecie-

ron los niveles de gris por debajo de un valor m y abrillantaron los niveles arriba de m de la imagen

original. Usualmente se definen dos valores (m y n) para que los niveles de gris de la imagen se distri-

buyan en un rango más controlado. En este caso los valores n y m fueron calculados automáticamente

con la función strechlim de Matlab. En el Apéndice D se encuentra más información de las funciones

de Matlab mencionadas previamente. En la Figura 3.5 se muestra la imagen de la Figura 3.6 después

de estirar su contraste, como se puede observar algunas partículas se distinguen con más facilidad.

Figura 3.5: Imagen con el contraste ecualizado.

Las imágenes empleadas en la técnica de PIV pueden contener además de partículas, ruido, objetos

etc. Estos objetos pueden estar fijos o en movimiento en la secuencia de imágenes, pero en ambos casos



deben ser eliminados para disminuir los errores en la siguiente etapa de procesamiento.

Los objetos fijos están presentes siempre en la misma posición en una secuencia de imágenes. Los

objetos más comunes son las paredes del recipiente que contiene al fluido (como el acrílico). Por ejemplo,

como se puede observar en la Figura 3.6, existe un objeto (delineado con una línea punteada) el cuál

es una buza por donde se hace pasar el fluido (Imagen del laboratorio de velocimetría de ESIQIE-IPN,

ver Apéndice A para más detalles del experimento). Con este tipo de objetos, el algoritmo de PIV

interpreta que son partículas sin movimiento.

Figura 3.6: Imágenes de PIV con objetos fijos.

La eliminación de los objetos fijos se puede realizar ya que PIV solo se utiliza para estudiar fluidos

en movimiento, por lo tanto, es poco probable que una partícula se encuentre en el mismo lugar en

todas las imágenes. De modo que no se eliminan partículas en esta etapa, solo los objetos no deseados

presentes y estáticos en todas las imágenes de un experimento.

Para poder hacer ésta eliminación se observó que eran necesarias imágenes adicionales al par original

a procesar. Se realizaron diferentes pruebas para determinar el número adicional de imágenes. En

diferentes experimentos se observó que con un par adicional era suficiente para determinar este tipo

de objetos en un experimento dado. Así se toman las cuatro imágenes dadas y se aplica una operación

lógica AND entre todas las imágenes. La matriz resultante es binarizada con un umbral (en este caso

128), esto quiere decir que los píxeles con valor menor a 128 se convierten en 0, mientras que los

mayores o iguales a 128 son sustituidos por un valor de 255.

En la Figura 3.7 se muestra la imagen resultante, la cual es utilizada como máscara. Esta máscara es



restada a la imagen creada a partir del ajuste de contraste de la imagen original, dando como resultado

una imagen sin los objetos no deseados estáticos como se muestra en la Figura 3.8.

Figura 3.7: Máscara de objetos estáticos.

Otro tipo de objetos presentes en algunas imágenes son aquellos que están en movimiento durante

una secuencia de imágenes. Como ejemplo se muestra la Figura 3.9 en la cual se tiene una imagen

de PIV con burbujas de aire suspendidas en un flujo laminar. Dicha imagen fue adquirida por el

laboratorio de velocimetría de ESIME-IPN en el 2001.

Los objetos en movimiento representan una fuente de error al determinar la velocidad del fluido,

pues puede que no sigan fielmente al flujo. Por lo tanto es necesario hacer un filtrado previo al análisis

de la imagen con el fin de utilizar imágenes en las cuales sólo se encuentren partículas. El filtrado se

realizó en dos etapas y se hace por medio de una máscara1 creada para cada imagen. Debido a que

los objetos se encuentran en movimiento, es necesario crear una máscara por cada imagen del par o

secuencia de imágenes.

En la primera etapa se crea la máscara y para construirla se aplica el filtro de mediana para eliminar

las partículas de la imagen. El filtro de mediana es muy utilizado para la eliminación de ruido aleatorio

1En informática, conjunto de datos que, junto con una operación, permiten extraer selectivamente ciertos datos alma-

cenados en otra.



Figura 3.8: Imagen de partículas sin objetos estáticos.

pues preserva la nitidez de la imagen. Este filtro elimina el ruido llamado “sal y pimienta” [21], el cual

consiste en píxeles de un color opuesto al que predomina en la imagen esparcidos de forma aleatoria.

En este caso las partículas son los píxeles a eliminar, las cuales equivalen a ruido sal (puntos blancos)

en la imagen. De este modo se obtiene una nueva imagen con los objetos no deseados (burbujas o

basura en el fluido). En la segunda etapa esta imagen es utilizada a su vez como máscara y es restada

de la imagen original. Como resultado de esta resta se obtiene una imagen cuyos únicos objetos son

las partículas.

Cada píxel de la imagen se reemplaza con la mediana del valor de los píxeles vecinos. La mediana se

calcula mediante el ordenamiento de todos los valores numéricos de los píxeles vecinos y se reemplaza

por el valor del píxel que se está en medio de todos los valores ordenados, es decir, si n es el número

de vecinos del píxel, la mediana será el valor en la posición n/2+1. Si el número n de valores tomados

es par, la mediana será el promedio del valor de los dos píxeles centrales, que serán los píxeles en las

posiciones n/2 y n/2 + 1. La Figura 3.10, se muestra una vecindad de nueve píxeles.

En el centro del vecindario se observa el píxel a sustituir (con un borde punteado) con el valor 218.

Por lo tanto los valores a considerar serán:



Figura 3.9: Imagen de partículas con objetos no deseados (burbujas).

Figura 3.10: Filtro de Mediana.



130, 125, 120, 129, 218, 120, 128, 127, 124.

Al ordenar estos valores se tiene:

Valor 120 120 124 125 127 128 129 130 218

Posición 1 2 3 4 5 6 7 8 9

Como se puede observar, el valor que está en medio (posición n/2+1 =9/2+1 = 5 ) del número de

valores considerados, es el 127 que será el valor por que se reemplazará el píxel con el borde punteado,

así en vez de tener el valor 218 tendrá 127, ya que se considera que el valor 127 es más representativo

de la región de estudio que el 218.

Por ejemplo el tamaño de este vecindario depende del tamaño de las partículas. Se puede tomar

como referencia el hecho de que si el área de los píxeles a eliminar es menor a n2±2 (la mitad del área

del filtro), los píxeles son eliminados con un filtro de mediana de con una vecindad cuadrad de n píxeles

por lado.

Así finalmente cuando se aplique a toda la imagen se obtendrá una imagen suavizada y sin par-

tículas, ya que en el caso de encontrar puntos blancos, estos serán eliminados, pues a comparación

de sus vecinos, estos píxeles no serán representativos de la región de estudio. En la implementación

del programa el vecindario se tomó de 3× 3, pero si las partículas de las imágenes a estudio son másgrandes, es decir que ocupen más píxeles, es importante que este vecindario sea modificado.

Para la implementación de esta etapa, se utilizó la función medfilt2, que aplica el filtro de mediana a

una imagen. Como se mencionó el objetivo es tener como resultado una imagen que no tiene partículas

como se muestra en la Figura 3.11.

La imagen obtenida se utiliza como máscara y es restada de la imagen original, dando como resulta-

do una imagen sin elementos ajenos a las partículas. En la Figura 3.12 se muestra la imagen resultado

de aplicar el pre-procesamiento previamente mencionado (ajuste del contraste, eliminación de objetos

fijos y objetos en movimiento).



Figura 3.11: Imagen con elementos ajenos sin partículas.

Figura 3.12: Imagen de partículas pre-procesada


3.2. Etapa de procesamiento.


En esta etapa se realiza el procesamiento de PIV para calcular el mapa con los vectores de velocidad.

La etapa de procesamiento propuesto en esta tesis se realiza con correlación cruzada por medio de la

Transformada Rápida de Fourier y cuenta con opciones como: traslape entre las zonas de interrogación,

interpolación subpíxel y proceso iterativo-multigrid.

Como se mencionó en el capítulo anterior, para realizar el procesamiento cada imagen es dividida

en pequeñas regiones cuadradas o rectangulares llamadas zonas de interrogación y estas zonas son

procesadas para determinar el desplazamiento del fluido y de las partículas entre las dos imágenes.

Este desplazamiento se determina por medio de la correlación cruzada, por lo que a las zonas de

interrogación de ambas imágenes, se les aplica la función fft2 de Matlab. Esta función se debe aplicar

a las zonas de ambas imágenes por separado. El resultado de esta función es la transformada de

una subimagen al espacio de frecuencia, la cual se utilizara para obtener la correlación cruzada como

se describió en la sección 2.5. A la transformada de la segunda subimagen se le aplica la función

conj, para obtener el conjugado de la matriz, el cual se multiplica por la transformada de la primera

subimagen. Finalmente con la función ifft2 se obtiene la transformada inversa de Fourier de la matriz

resultante para así obtener el plano de correlación cruzada y así determinar el máximo que representa

el desplazamiento del fluido en esa zona de la imagen. Este desplazamiento se representa con un vector

y el proceso se repite para todas las zonas de las imágenes.

En el algoritmo básico de PIV el traslape entre las zonas de interrogación de una imagen es cero;

es decir, la separación entre dos zonas de interrogación es precisamente el tamaño de las zonas. De este

modo, donde termina una zona empieza la siguiente y no existe intersección entre ellas. El problema de

utilizar este traslape de cero es que las partículas que se encuentran en las orillas de la zona de la primera

imagen difícilmente se encontraran en la segunda zona. Esto es debido al movimiento de las partículas

en las imágenes del fluido. Para solucionar este problema y además obtener una mayor resolución

de vectores, el algoritmo permite utilizar un traslape del 0 al 100% horizontal y verticalmente. Este

porcentaje se aplica al tamaño de la zona de interrogación de forma independiente en la dirección x,

y. Típicamente se utiliza el traslape de 50% ó 75% y es el mismo porcentaje en ambas direcciones.

En la Figura 3.13 (a) se muestra una imagen con cuatro zonas de interrogación sin traslape, lo

que da como resultado cuatro vectores, uno en el centro de cada zona de interrogación. En la figura

(b) las zonas de interrogación tienen el mismo tamaño de las zonas de la figura (a), pero éstas están

traslapadas en un 50% de su tamaño en las direcciones x e y. Al tener este traslape se obtienen más

vectores de velocidad pues se utilizan 9 zonas de interrogación en lugar de 4. Como se puede observar



el primer vector se encuentra en la misma posición sin importar si existe o no traslape, y se utiliza el

traslape a partir del segundo vector.

Figura 3.13: (a)Imagen sin traslape, (b) imagen con traslape del 50%

La evaluación de las imágenes de PIV se realiza de forma digital y los datos de entrada son dis-

cretizados. Por lo tanto los valores de correlación cruzada son enteros, pues el desplazamiento de

las partículas se mide en píxeles. Así el desplazamiento podrá tener una imprecisión de ±1/2 píxel.Sin embargo en la literatura se ha reportado que es posible reducir esta imprecisión, al aproximar el

desplazamiento utilizando una escala a nivel subpíxel.

Para obtener una precisión a nivel subpíxel es necesario conocer las coordenadas del máximo y sus

vecinos en el plano de correlación cruzada. Con esto se trata de ajustar una curva para estimar la

posición real del máximo de correlación cruzada y así obtener mayor precisión.

Existen diferentes formas de estimar la posición del pico de correlación cruzada y por lo tanto

aproximar el desplazamiento. En este algoritmo se implementaron las propuestas en la referencia [1]

que son de ajuste gaussiano, centroide y parabólico. Los estimadores utilizados se enlistan en la Tabla

3.1.

El ajuste centroide es usado frecuentemente. Esta definido como el cociente del momento de primer

orden y el de orden cero. Otro método es el ajustar los datos de correlación cruzada a una función, ya

sea Gaussiana o Parabólica. Lo más común es emplear tres valores consecutivos, incluido el máximo,

para estimar el desplazamiento. El ajuste Gaussiano es el que con más frecuencia se implementa.

La interpolación a nivel subpíxel por omisión no se realiza y en caso de hacerlo solo se aplica a una

iteración (la última en caso de que el procesamiento sea iterativo).

El algoritmo de correlación cruzada además se modificó para que contara con la opción de que

el procesamiento se hiciera de forma iterativa. En cada iteración el tamaño de las zonas de interro-

gación se reduce en un porcentaje determinado (en este caso, 50%) en ambas direcciones (x, y). El



Función de ajuste Estimadores

Pico centroide x0 =(i− 1)V(i−1,j) + iV(i,j) + (i+ 1)V(i+1,j)

V(i−1,j) + V(i,j) + V(i+1,j)

f(x) = momento de orden unomomento de orden cero x0 =

(i− 1)V(i−1,j) + iV(i,j) + (i+ 1)V(i+1,j)

V(i−1,j) + V(i,j) + V(i+1,j)

Ajuste de pico parabólico y0 =(j − 1)V(i,j−1) + jV(i,j) + (j + 1)V(i,j+1)

V(i,j−1) + V(i,j) + V(i,j+1)

f(x) = Ax2 +Bx+ C y0 = j +V(i,j−1) − V(i,j+1)

2V(i,j−1) − 4V(i,j) + 2V(i,j+1)

Ajuste de pico Gaussiano x0 = i+lnV(i−1,j) − lnV(i+1,j)

2 lnV(i−1,j) − 4 lnV(i,j) + 2 lnV(i+1,j)

f(x) = C exph−(x0−x)2

k

iy0 = j +

lnV(i,j−1) − lnV(i,j+1)2 lnV(i,j−1) − 4 lnV(i,j) + 2 lnV(i,j+1)

Tabla 3.1: Estimadores de tres puntos para determinar el desplazamiento de los datos de correlación a

nivel subpíxel

procedimiento da como resultado una mejora en la resolución espacial obtenida con cada iteración.

Este algoritmo se conoce como PIV con Multigrid o Multiresolución [20], debido a que se utiliza una

diferente resolución con cada iteración. Después de cada iteración se reduce el tamaño de las zonas y

se utiliza la información del mapa de vectores (dirección y magnitud) obtenido, para poder determinar

la posición de las zonas de interrogación de la siguiente iteración. El algoritmo iterativo se implementó

para que se puedan iterar de 1 a 4 veces, según sea requerido. Adicionalmente, en cada iteración se

validan y se sustituyen los vectores, para que no se introduzca error en iteraciones posteriores, ya que

la información el desplazamiento obtenido, es un dato muy importante en iteraciones posteriores.

En la Figura 3.14 se ejemplifica el proceso iterativo-multigrid. En la primera iteración se determinó

un desplazamiento hacía el extremo inferior derecho (denotado por la flecha), el cuál será tomado como

referencia en la segunda iteración para desplazar la posición de la zona de interrogación de la segunda

imagen. El mapa de vectores resultante es utilizado de la misma forma para desplazar las zonas en la

siguiente iteración.

En la Figura 3.15 se muestra el resultado del procesamiento básico de PIV sin pre-procesamiento a

la imagen de la Figura 3.2 y a su pareja. A comparación se muestra la Figura 3.16 que es el resultado del

procesamiento con 3 iteraciones y con pre-procesamiento. El par de imágenes elegido para esta prueba

es muy difícil de evaluar, pues la velocidad del fluido es muy alta y la iluminación no es uniforme. Por

esta razón los resultados no son tan exactos.



Figura 3.14: Proceso iterativo con multiresolución. Reducción del 50% del tamaño de las zonas de

interrogación con cada iteración.

Figura 3.15: Mapa de vectores sin pre-procesamiento.


3.3. Etapa de post-procesamiento.

Figura 3.16: Mapa de vectores con pre-procesamiento.


Una vez que se tienen los datos de PIV estos pueden ser procesados antes de realizar un análisis.

Esta etapa llamada post-procesamiento y se divide en los siguientes pasos:

Validación de los datos.

Después de la evaluación de las imágenes se puede distinguir a simple vista un cierto número de

vectores incorrectos. Con algoritmos especiales se debe realizar una validación automática para detectar

y eliminar los vectores incorrectos.

Reemplazo de los datos incorrectos.

Para poder hacer un análisis correcto de los resultados de PIV es necesario contar con el mapa de

vectores completo, por lo tanto es necesario reemplazar con nuevos vectores todos los espacios de los

vectores que fueron eliminados durante la validación.

Reducción de datos.

Es complicado analizar los cientos de vectores resultantes para poder describir las características

dinámicas del fluido. Con el fin de simplificar el estudio del flujo, se pueden utilizar ciertas técnicas

para obtener estadísticas o datos importantes como velocidad media.

Presentación y animación de la información.

Esta etapa es de utilidad pues ayuda al entendimiento de las características principales del flujo.

La animación de los datos es muy útil para el entendimiento de series de imágenes y su cambio con el

tiempo.



Prácticamente todos los mapas de vectores calculados con PIV tienen vectores incorrectos, y tra-

bajar con estos vectores pueden generar interpretaciones incorrectas. En esta etapa se realiza una

validación del mapa de vectores que se obtuvo como resultado de la etapa anterior. El algoritmo de

validación consta a su vez de dos etapas: determinación de vectores incorrectos y su sustitución. La

primera etapa se realiza por medio de la comparación de cada uno de los vectores con sus vecinos

adyacentes. Si la diferencia entre un vector en magnitud y dirección con respecto a sus vecinos, es

mayor a la de un umbral determinado, el vector es eliminado. El umbral que se utiliza es un depende

del promedio de la desviación estándar de la magnitud y dirección de cada vector y sus vecinos. Poste-

riormente los vectores eliminados serán remplazados por unos nuevos calculados a partir del valor de

sus vecinos.

En la Figura 3.17 se muestra un mapa típico de vectores antes y después de la validación y susti-

tución. En la Figura 3.17 (a), se observan englobados en círculos, los vectores incorrectos, en la Figura

3.17 (b), estos vectores fueron eliminados, y posteriormente sustituidos (c).

El algoritmo para la validación de vectores debe asegurar que al final de la evaluación, los vectores

que permanezcan serán correctos. Actualmente existen varias técnicas de validación de datos que

buscan detectar valores incorrectos. Una forma de validar el mapa de vectores es revisar cada uno

de los vectores y compararlo con los que se encuentren en su vecindad, este proceso se describe a

continuación.

Sea V una matriz de vectores. Se define la q − vecindad del vector V (i, j), como Nq [V (i, j)] , al

conjunto de los q vectores vecinos más cercanos al vector V con coordenadas (i, j). Es decir para un

vector V con 8 vecinos, su 8− vecindad está determinada por:

N8 [V (i, j)] =

⎧⎪⎪⎨⎪⎪⎩V ( i+ 1, j), V ( i− 1, j ), V ( i, j + 1 ), V ( i, j − 1),

V ( i+ 1, j + 1), V ( i+ 1, j − 1 ), V ( i− 1, j + 1 ), V ( i− 1, j − 1)

⎫⎪⎪⎬⎪⎪⎭ (3.1)

En la Figura 3.18 se muestra la matriz que representa el mapa de vectores con la 8− vecindad. Esta

vecindad es muy importante ya que es la más utilizada durante la validación y sustitución. La distancia

en píxeles entre el vector central V (i, j) y sus vecinos es ∆XZONA,∆YZONA, oq∆X2

ZONA +∆Y2ZONA.

Donde ∆XZONA y ∆YZONA es el tamaño en píxeles en x, y respectivamente, de las zonas de interro-

gación.

Como se mencionó, en la validación de vectores generalmente se utiliza una vecindad N8 [V (i, j)], a

excepción de que V esté en algún borde de la matriz. En este caso, sólo se podrán tener 5 o 3 vectores

vecinos cuando el vector se encuentra en una esquina de la matriz. Una vez que se tienen los vectores



Figura 3.17: (a) Mapa con vectores incorrectos (b) mapa sin vectores incorrectos (c) mapa con vectores

sustituidos.

Figura 3.18: Matriz de vectores.



vecinos, se hace una comparación del vector central y las estadísticas de sus vecinos. Los métodos más

comunes de validación son: Media global, Media local, Mediana local e Histograma global [1] y [22].

La desviación estándar entre el vector central y sus vecinos es comparada con un umbral previamente

definido por el usuario. Si un vector sobrepasa este umbral entonces es considerado como un vector

incorrecto y es eliminado. En la literatura existen diferentes propuestas de cómo definir este umbral,

una opción es definir un porcentaje de diferencia entre un vector y la desviación estándar de sus vecinos.

Los vectores eliminados son reemplazados con nuevos vectores creados a partir de la interpolación de

sus vecinos. El método más utilizado es el de interpolación bicúbica de sus 4-vecinos, aunque esto solo

funciona si estos vecinos son correctos y no fueron eliminados. De otro modo se tiene que utilizar otro

método para sustituir los vectores. Finalmente el mapa de vectores de velocidad completo, resultado

de utilizar la técnica de PIV, es analizado y su interpretación dependerá del experimento que se haya

realizado.

Como se menciono en el capítulo anterior, un vector puede tener entre 3 y 8 vectores vecinos. Sin

embargo, en esta etapa del procesamiento es posible que no se tenga ningún vector vecino, como se

muestra en la Figura 3.19. En esos casos se utiliza la información de los vectores vecinos más cercanos,

aunque no se encuentren en su vecindad de 8, ya que se redefine la vecindad, al hacerla más grande

pero con la restricción de que no se encuentren a una distancia mayor a 3 veces ∆XZONA, ∆YZONA oq∆X2

ZONA +∆Y2ZONA.

También se tiene la restricción de que al menos se deben de tener la mitad de los vectores vecinos

para poder calcular un nuevo vector. En caso de que en la primera iteración no se reemplacen todos

los vectores eliminados previamente, el proceso se repite para estos vectores hasta llenar por completo

el mapa.

En la implementación de esta etapa se utiliza el promedio y desviación estándar de cada vector

y sus vecinos. Después de realizar varias pruebas con diferentes imágenes del laboratorio de reología,

se utilizó el criterio siguiente: si un vector cae fuera del intervalo definido por el promedio ±1,75la desviación estándar entonces éste es eliminado. Sin embargo, es necesario ajustar el tamaño del

intervalo de acuerdo a los resultados de cada experimento.

Es posible hacer la validación sin tener que reemplazar los vectores incorrectos, los cuales son

sustituidos por un vector con velocidad cero en sus dos componentes. Sin embargo por obvias razones

para poder hacer la sustitución es necesario haber hecho previamente la validación de los vectores

incorrectos. Estos procedimientos son opcionales y se realizan al final del procesamiento. En caso de

que el procesamiento sea iterativo, la validación y sustitución se realizará al final de cada iteración.



El resultado de la validación aplicada a los mapas de vectores de las figuras de la sección anterior

(Figura 3.15 y Figura 3.16) se muestra en la Figura 3.20 y Figura 3.21. Posteriormente se muestran

los mapas resultado de sustituir estos vectores incorrectos en la Figura 3.22 y Figura 3.23.

Después de terminar estas etapas, se calculan el resultado final a partir del mapa de vectores

obtenido. Este resultado es la matriz del campo de velocidades(calculada con el tiempo entre imágenes

y la distancia en píxeles). Adicionalmente, se incluyó una opción para permitir guardar estos resultados

en un archivo de texto.



Figura 3.19: Mapa con varios vectores a reemplazar.

Figura 3.20: Mapa de vectores sin pre-procesamiento, con validación.



Figura 3.21: Mapa de vectores con pre-procesamiento, con validación.

Figura 3.22: Mapa de vectores sin pre-procesamiento, con validación y sustitución.

Figura 3.23: Mapa de vectores con pre-procesamiento, con validación y sustitución.


CAPÍTULO 4

Pruebas y Resultados

En este capítulo se muestra un comparativo de los resultados obtenidos con el algoritmo aquí

propuesto y con otros programas. Durante la elaboración de esta tesis se realizaron pruebas a distintos

grupos de imágenes:

Imágenes generadas de forma artificial

Imágenes de experimentos del IPN:

• Laboratorio de Reología, ESFM. En este laboratorio utilizan la técnica de PIV para estudiarflujos de fluidos complejos.

• Laboratorio de Velocimetría, ESIQIE. Aquí estudian el comportamiento del acero líquidocuando pasa a través de distintas buzas y buscan diseñar buzas que distribuyan de forma

uniforme el fluido.

• Laboratorio de ESIME. En este laboratorio utilizan PIV para determinar el movimiento deburbujas dentro de fluidos.

De estos grupos se seleccionaron las imágenes de un experimento del laboratorio de Reología, por

ser las únicas de las cuales se cuenta con los resultados teóricos del movimiento del fluido. Para realizar

la comparación de los mapas de vectores, se tomaron como referencia los datos obtenidos con el modelo.

Además de las pruebas con el algoritmo propuesto, las imágenes fueron evaluadas con otros tres

programas gratuitos. Estas pruebas se realizaron con el fin de comparar los resultados obtenidos con

los distintos programas. Los paquetes gratuitos son Urapiv [23], Matpiv [24],[25] y Mpiv [26].Estos se

distribuyen en Internet en un paquete para Matlab llamado toolbox. Un toolbox en Matlab es una

serie de funciones que trabajan para un mismo fin (como lo es una librería en C/C++ o un paquete

4.1. Descripción de los programas de PIV evaluados.

en Java). Las pruebas fueron realizadas con los parámetros utilizados con el programa de Dantec, sin

embargo, cuando el programa no aceptó esos mismos parámetros se utilizaron unos semejantes.

En la Figura 4.1 y la Figura 4.2 se muestran las dos imágenes de partículas del laboratorio de

Reología, las cuales se muestran a un 45% de su tamaño original que es de 1024× 768.

El sistema de velocimetría para este experimento consistió de una cámara CCD, la fuente de

luz basada en dos láseres, una interfaz de control y análisis Dantec 1100, interconectado por una

computadora personal y validado por el software FlowMap de Dantec.

La cámara empleada fue una Kodak ES 1.0 de 30 Hz con una resolución de 1.1 mega píxeles, a la

cual se le asoció un sistema de lentes AF Micro Nikkor. La cámara se colocó a una distancia aproximada

de 0.085 m del área a capturar. El tiempo entre dos fotografías consecutivas varió desde 25 ms hasta

0.319 ms. El tiempo entre cada par de fotografías de acuerdo al régimen de flujo estudiado, variando

este tiempo desde 100 ms hasta 2000 ms. Las fotografías obtenidas cubrieron un área transversal de

0.00673 m x 0.003 m aproximadamente.

La fuente de luz consistió de dos láseres Nd:YAG de 0.1 Joules de potencia. Estos estaban sin-

cronizados para enviar pulsos de 6 ns de duración con una longitud de onda λ = 525 nm. En este

experimento se utilizó un espesor del plano de luz de menos de 0.0005m.

Las partículas utilizadas y sembradas en la solución fueron esferas huecas de vidrio recubiertas de

plata de 10 μm de diámetro elaboradas por Dantec. Estas partículas fueron añadidas a la solución de

CPyCL/NaSal en dos concentraciones, 0.0025% y 0.00375% en peso. El área de estudio consistió de

dos tubos capilares de vidrio de 0.002942 m de diámetro.


4.1.1. FlowMap 2000 (Dantec Dynamics, 2000).

Dantec Dynamics es una compañía dedicada al desarrollo y venta de sistemas integrados de medición

para diagnóstico e investigación de dinámica de fluidos, micro fluidos, análisis de spray, tecnología de

combustión y materiales.

Entre sus variados productos se encuentra FlowMap, este es un sistema completo de productos que

han sido diseñados y construidos para obtener mediciones instantáneas de los campos de velocidad


Capítulo 4. Pruebas y Resultados

Figura 4.1: Imagen de partículas 1

Figura 4.2: Imagen de partículas 2, tomada 25 ms después de la primera imagen



con la técnica de PIV. El sistema se basa un procesador de PIV, el cuál es un equipo sofisticado de

hardware que permite la interconexión de todo los componentes de PIV, es decir, la cámara, el láser

y en este caso la computadora, con ayuda de un software llamado FlowManager realiza dos tareas:

sincronizar el láser y la cámara y validar los vectores analizados. Una vez que se inicia un experimento

la computadora sincroniza los componentes, las imágenes de partículas son capturadas y se realiza

el procesamiento de PIV básico por medio de hardware especializado. Después estas imágenes son

transferidas a la computadora para almacenarlas y realizar algún post-procesamiento opcional. Los

resultados obtenidos con este sistema son buenos y rápidos ya que la mayoría del procesamiento es

llevado a cabo por medio del hardware.

El mapa de vectores obtenido con las imágenes de prueba de la Figura 4.1 y la Figura 4.2 se muestra

en la Tabla 4.2 de la sección 4.2. Este mapa de vectores base fue calculado por medio de una iteración

de correlación cruzada con zonas de interrogación de 256 píxeles de ancho por 32 píxeles de alto, y con

un 75% de traslape entre las zonas de interrogación en ambas direcciones.

4.1.2. Matpiv.

Para el procesamiento de PIV, se puede definir una máscara de procesamiento y el sistema de coor-

denadas. El procesamiento se puede hacer de forma iterativa, aunque no siempre da buenos resultados.

Además tiene la limitante de que las zonas de interrogación tienen que ser cuadradas. Por esta razón

las pruebas con este programa no se realizaron con los mismos parámetros utilizados con el programa

de Dantec.

Entre las funciones adicionales al procesamiento PIV de las imágenes, se encuentra la opción de va-

lidar resultados por diferentes métodos como por medio de la mediana local (utilizado en este trabajo),

validación por el operador de histograma global y por medio del filtro SNR (Signal to Noise Ratio).

El programa tiene un procesamiento rápido, sin embargo, tener que teclear o hacer un programa para

poner todos los comandos necesarios para procesar y visualizar los resultados puede ser tedioso para

el usuario, ya que no es sencillo hacer la configuración necesaria para que funcione el programa. El

programa falló en muy pocas ocasiones, sin embargo, no fue posible saber cuál era el error pues el

mensaje de error no proporcionó la suficiente información. Las pruebas fueron realizadas con zonas de

interrogación cuadradas de 256 por 256 píxeles.



4.1.3. Urapiv.

A diferencia de los demás, este programa ha sido actualizado ya en varias ocasiones, la última en

el mes de enero del 2006. En versiones anteriores se incluía una interfaz de usuario muy sencilla de

utilizar, en donde se podían introducir los parámetros para el procesamiento; sin embargo en la última

versión ya no se incluye.

A comparación de MatPIV, el uso de UraPIV es más sencillo, ya que en los parámetros de la

función principal se incluyen las opciones para graficar los resultados, seleccionar la región a procesar

y algunos parámetros para la validación de vectores (con el filtro de SNR). Sin embargo, este programa

no incluye procesamiento iterativo.

Al igual que MatPIV, este programa tiene la limitante de que las zonas de interrogación tienen que

ser cuadradas, por lo que las pruebas tampoco se realizaron con los mismos parámetros utilizados con

el programa de Dantec.

Este programa falló con las imágenes que incluían burbujas y con algunas del laboratorio de reología,

en las demás imágenes de prueba obtuvo buenos resultados en la mayoría de los casos. Si se quieren

procesar varias imágenes, este programa recibe como parámetro la ruta del directorio en donde se

encuentran las imágenes a procesar. En conclusión es sencillo de usar y obtiene resultados aceptables,

pero a pesar de los años que ya lleva, no tiene muchas opciones de procesamiento y validación de datos.

4.1.4. Mpiv.

Al igual que las aplicaciones anteriores Mpiv, no presenta una interfaz de usuario. Así, las principales

funciones son dos: la que realiza el procesamiento y la que valida los vectores incorrectos.

En los parámetros de la función principal resalta la opción de utilizar un algoritmo diferente al

de correlación cruzada para determinar el desplazamiento de las partículas. Este algoritmo se basa

en las Diferencias de Mínimos Cuadrados (Minimum Quadratic Differences, MQD). Sin embargo, el

algoritmo es lento y los resultados que presenta no son correctos. También se tiene la opción de

hacer el procesamiento iterativo. Al finalizar el procesamiento verifica un parámetro para graficar o

no, los resultados. Aunque el programa tiene la opción de procesar con el método de MQD, éste no

presenta buenos resultados. Además para poder utilizar este toolbox es necesario otro toolbox llamado

DACE [27], Design and Analysis of Computer Experiments (Diseño y Análisis de Experimentos de

Computadora), del cuál se utilizan las funciones para la correlación cruzada y la interpolación.



La segunda función obtiene los datos del procesamiento, es decir, la información de la matriz de

vectores. Permite seleccionar tres tipos de filtrado (global, por mediana ó estándar), del mismo modo

permite seleccionar el modo de interpolación de los vectores a sustituir (sin interpolación, interpolación

cúbica, interpolación de Kringing). En esta última opción los autores recomiendan la interpolación de

Kringing.

El programa falló casi en todas las pruebas que se realizaron, pues dependía mucho de la calidad

de las imágenes utilizadas. El mensaje de error no proporciona información para conocer la causa del

error. Aunque los resultados no son los que deberían, cabe mencionar que este fue el único programa

gratuito con el que pudieron realizar las pruebas con esos mismos parámetros, al permitir zonas de

interrogación rectangulares.

4.1.5. PIV con el algoritmo propuesto

Este programa es la implementación del algoritmo propuesto en esta tesis. Para el pre-procesamiento

se puede definir una máscara de procesamiento, modificar el contraste y eliminar los objetos no deseados

en las imágenes. El procesamiento se puede hacer de forma iterativa y con interpolación a nivel subpíxel

para obtener una mayor resolución espacial de resultados, además de incrementar su precisión. En este

programa no existe la restricción de que las zonas de interrogación sean cuadradas, por lo que si así

se desea éstas pueden ser rectangulares. Entre las funciones adicionales al procesamiento PIV de las

imágenes, se encuentra la opción de validar y sustituir los resultados incorrectos.

El programa cuenta con una interfaz gráfica sencilla de utilizar, en la cual se introducen los pa-

rámetros de las tres etapas de procesamiento. Además se cuenta con la opción de graficar el mapa de

vectores y guardarlo en un archivo si así se desea. El programa en general dio buenos resultados, aún

con imágenes que incluían burbujas o tenían una geometría compleja.

Las pruebas se realizaron utilizando distintos parámetros, incluyendo la misma configuración que

la utilizada con el programa de Dantec (zonas de interrogación rectangulares de 256 por 32 píxeles, con

un traslape del 75% entre las zonas de interrogación). La principal diferencia fue el pre-procesamiento

y los resultados son considerablemente mejores cuando se aplica esta etapa.

Los resultados de los programas descritos anteriormente dependen demasiado de la calidad de las

imágenes. Por eso es tan importante la etapa de pre-procesamiento, para que se pueda mejorar de

manera significativa la calidad de las imágenes de entrada y a su vez la calidad de los resultados

obtenidos. Además, al contar con una etapa de pre-procesamiento es posible seleccionar la región de



interés a ser evaluada y ahorrar tiempo de procesamiento.

4.2. Resultados.

Como se señaló anteriormente, el experimento estudiado corresponde al uso de un capilar donde se

hace pasar a través de éste un flujo laminar, que corresponden a las Figuras 4.1 y 4.2.

Las comparaciones con los programas se hicieron de acuerdo al modelo que describe el flujo del fluido

en el capilar. Para la rapidez de deformación empleada, el fluido se comporta como un Newtoniano. El

perfil de velocidades está dado por:

u(r) =τaa

2μ

µ1− r2

a2

¶(4.1)

Donde τa es el esfuerzo en la pared del capilar, μ es la viscosidad del fluido, a es el radio total del

capilar y res la distancia del centro del capilar a cualquier punto de la zona de estudio. Así la ecuación

dará la magnitud de la velocidad del vector en función del radio del capilar en estudio. Para el caso

del experimento a estudio τa = 1,6944Pa, μ = 0,9279 Pa s y a= 1.45mm.

Todos los resultados obtenidos con los diferentes programas se compararon con este modelo. Para

cada vector obtenido con los programas se utilizó su localización (x, y) para obtener la magnitud del

vector en cada punto según el modelo descrito anteriormente.

El mapa de vectores obtenido con el modelo 4.1 se muestra en la Figura 4.3:

Como se puede observar, cada columna del mapa de vectores es una parábola, es decir, la velocidad

del fluido es mayor en el centro y en las paredes del capilar es cero. Debido a que es un flujo estacionario,

en esta geometría el perfil de velocidades es el mismo a lo largo del capilar. Por lo tanto se puede obtener

una parábola promedio, que describa la velocidad del flujo, ésta se muestra en la Figura 4.4:

En la Figura 4.4 los vectores parecen tener diferente magnitud al mapa de la Figura 4.3, esto es

porque la función utilizada para graficar el mapa de vectores realiza un reescalamiento de todos los

vectores a desplegar. Se toma el vector con mayor magnitud y le da cierta proporción con respecto a

los píxeles desplegados (200× 400en este caso); los demás los despliega en escala a éste. Se mencionaesto ya que en algunas comparaciones aunque en la gráfica parezca que los vectores no tienen una

magnitud similar, lo importante es la desviación estándar de las dos parábolas promedio obtenidas.

A continuación se presentan tablas con los datos resultantes entre los programas evaluados y el

modelo 4.1. Estas pruebas sirven para comparar la diferencia en resultado conforme al modelo analítico.


4.2. Resultados.

Figura 4.3: Mapa de vectores resultante del modelo.



Figura 4.4: Parábola promedio del modelo.

En el Apéndice F se encuentran algunas pruebas realizadas con diferentes fluidos y geometrías de

diversas aplicaciones. En el caso del programa Mpiv no se obtuvieron resultados comparables, ya que

los vectores obtenidos eran NaN .

En la Tabla 4.4 se muestran los resultados del programa Matpiv y su comparación con el modelo 4.1.

En este caso debido a las limitantes del programa, las zonas de interrogación fueron cuadradas. Como

se puede observar por la desviación estándar de esta comparación, los resultados obtenidos presentan

un error más grande en comparación de los programas aquí presentados.


4.2. Resultados.


4.2. Resultados.

Tabla 4.3: Comparación de desviaciones estándar por regiones.


CAPÍTULO 5

Conclusiones y Trabajo Futuro.

En esta tesis se presentó el diseño e implementación de un algoritmo para determinar el campo

de velocidad de un fluido en movimiento por medio de la técnica de PIV. Este es un primer

intento de realizar un programa de PIV de propósito general para el estudio de diversos flujos

de fluidos. El algoritmo presentado consta de tres etapas: pre-procesamiento, procesamiento y

post-procesamiento.

El sistema implementado ofrece la posibilidad de mejorar las imágenes y acondicionar su tamaño

antes de aplicar el algoritmo para obtener el campo de velocidades. El programa esta diseñado

para aplicar únicamente PIV con el algoritmo más básico si así se desea. En caso de necesitar

una mayor resolución, se encuentra implementado el procesamiento iterativo-multigrid y permite

hacer una interpolación a nivel subpíxel. Adicionalmente la aplicación permite validar los vectores

resultantes. Se tienen las opciones de graficar el mapa de vectores y de guardarlo en un archivo

de texto, para visualizar los resultados y facilitar el posterior uso de ellos. Las etapas agregadas

de pre-procesamiento y post-procesamiento son importantes, ya que gracias a ellas se obtienen

mejores resultados, tanto en precisión como en resolución espacial. Con lo cuál se sugiere que los

sistemas de PIV deben incluir estas etapas.

La interfaz de usuario, facilita el uso de la aplicación ya que en muchos de los programas que se

encuentran en Internet (véase el Capítulo 4) no incluyen ésta y su uso se complica. Además si así

se desea, se puede realizar el procesamiento sin el uso de la interfaz al llamar la función principal

de la aplicación.

El algoritmo empleado para generar el campo de velocidades da mejores resultados comparado con

del algoritmo básico. El algoritmo está implementado en una aplicación para Matlab y a su vez,

esta genera en la mayoría de los casos mejores resultados en comparación con los otros sistemas de

su tipo. Las pruebas del algoritmo aquí propuesto con las imágenes utilizadas (tanto artificiales

como reales) muestran mejores resultados que con los otros programas gratuitos, como se mostró

en el capítulo anterior. En estos programas, con algunas de las imágenes utilizadas, se obtuvieron

resultados con un gran margen de error, algunos de hasta el 100%. El sistema comercial de Dantec

[28], utilizado en los experimentos del Laboratorio de Reología del IPN, genera resultados con

un error del 5% en comparación con los obtenidos analíticamente. Con estos experimentos el

sistema aquí presentado obtiene resultados que difieren en un 7% de los resultados analíticos.

Como se puede observar, para un flujo en capilar el algoritmo genera resultados muy similares

a los de Dantec, aunque ambos tienen problemas en las fronteras del capilar. Por lo que este

nuevo sistema representa una buena opción gratuita para el procesamiento de PIV, aunque no

se cuenta con todas las opciones que tiene el sistema comercial.

En un futuro se tiene que solucionar el problema que se presenta cuando las zonas de interrogación

no están completas debido al tamaño de la imagen y a las fronteras físicas del área de estudio.

También se debe mejorar el algoritmo de validación y sustitución pues no en todos las pruebas

genera los resultados correctos. Una vez que se tengan los datos correctos, éstos pueden ser

utilizados para calcular características importantes del flujo (gastos, esfuerzos, gradientes, etc.).

La resolución de los resultados es otro punto a mejorar, pues sería útil realizar el seguimiento

individual de las partículas con la técnica de Súper-Resolución PIV (PIV con súper resolución o

PIV con PTV). El problema principal de esta técnica es encontrar una forma eficiente y precisa de

detectar todas y cada una de las partículas. Además requiere determinar los pares de partículas

correctamente.

Otra mejora a realizar es aprovechar el acceso a hardware que provee Matlab por medio de los

controladores necesarios para la cámara y el láser. Esto con el fin de integrar el sistema de PIV con

el hardware de captura, de modo que se pueda controlar todos los componentes del experimento

y así obtener un sistema completo de PIV.

Para mejorar el desempeño se puede usar el lenguaje de programación Java para realizar el pro-

cesamiento en diferentes hilos de ejecución y así acelerar el tiempo de obtención de los resultados

y tal vez generar una visualización de los vectores de velocidad en tiempo real.


REFERENCIAS

[1] M.Raffel, C.Willert, and J.Kompehans. PIV: A Practical Guide. Springer, 1st edition edition,

2002.

[2] R.J.Adrian. Twenty years of particle image velocimetry. 12th International Symposium on Appli-

cations of Laser Techniques to Fluid Mechanics, Lisbon, 2004.

[3] P.Hochareon, K.B.Manning, A.A.Fontaine, S.Deutsch, and J.M.Tarbell. Development of high

resolution particle image velocimetry for use in artificial heart research. In Proceedings of the

Second Joint EMBS/BMES Conference, 2002.

[4] J.T.Kelley and A.K.Prasad. Detailed flow patterns in the nasal cavity. Journal of Applied Phy-

siology, 89:323—337, 2000.

[5] S.K.Kim and S.K.Chung. An investigation on airflow in disordered nasal cavity and its corrected

models by tomographic piv. Measurement Science and Technology, 15:1090—1096, 2004.

[6] In-Bok Lee, K.Choi, J.Yun, J.Jeun, and G.Kim. Study of aerodynamics in agriculture - modern

tecnologies. Agricultural Engineering International: the CIGR Journal of Scientific Research and

Development, 5, 2003.

[7] T.Zhang and S.W.Van Sciver. Use of the particle image velocimetry technique to study the

propagation of second sound shock in superfluid helium. Physics of Fluids, 16:99—102, 2004.

[8] H.Okada, Y.Okuda, H.Kikitsu, and M.Ohashi. Strategy on wind engineering researches in bri

and nilim. Technical Report, Building Research Institute and National Institute for Land and

Infrastructure Management, 2002.

Referencias

[9] I.Grant, M.Mo, X.Pan, P.Parkin, J.Powell, and D.Hurst. A piv and lsv study in the wake of an

aircraft model. Journal of Visualization, 2(1):85—92, 1999.

[10] Marín-Santibañez B. Rheo-piv of shear thickening wormlike micells. Unpublished Work, 2005.

[11] G.R.Spedding, Hedenström, and M.Rosén. Quantitative studies of the wakes of freely flying birds

in a low-turbulence wind tunnel. Experiments in Fluids, 34(2):291—303, 2003.

[12] U.K.Muller, E.J.Stamhuis, and J.J Videler. Hydrodynamics of unsteady fish swimming and the ef-

fects of body size: Comparing the flow fields of fish larvae and adults. The Journal of Experimental

Biology, 203:193—206, 2000.

[13] Oxford Lasers. http://www.oxfordlasers.com. Electronic Citation, 2005.

[14] Big Sky Laser Technologies. http://www.bigskylaser.com. Electronic Citation, 2005.

[15] R.González. Digital Image Processing. Prentice-Hall, 2nd edition, 2002.

[16] R.J.Adrian. Particle-imaging techniques for expiremental fluid mechanics. Annual Review of Fluid

Mechanics, 23:261—304, 1991.

[17] K.Okamoto, S.Nishio, T.Saga, and Kobayashi. Standard images for particle-image velocimetry.

Measurement Science and Technology, 11:685—691, 2000.

[18] A.K.Prasad. Particle image velocimetry. Current Science, 79(1):51—60, 2000.

[19] Mark P.Wernet. Fuzzy logic enhanced digital piv processing software. Report, 1999.

[20] F.Scarano and ML.Riethmuller. Iterative multigrid approach in piv image processing with discrete

window offset. Experiments in Fluids, pages 513—523, 1999.

[21] A.Marion. An introduction to image processing. Chapman and Hall, 1991.

[22] J.Westerweel. Digital Particle Image Velocimetry-Theory and Application. Thesis/dissertation,

Delft University, 1993.

[23] URAPIV. www.urapivpress.com. Generic, 2006.

[24] J.Kristian Sveen. An introduction to matpiv v.1.6.1. Report, 2004.

[25] MatPIV. www.math.uio.no/ jks/matpiv. Generic, 2006.

[26] MPIV. sauron.civil.eng.osaka-cu.ac.jp/ mori/softwares/index.php. Generic, 2006.

[27] H.B.Nielsen. Dace. Electronic Citation, 2006.


Referencias

[28] Dantec Dynamics. http://www.dantecdynamics.com. Electronic Citation, 2005.

[29] TSI Incorporated. http:www.tsi.com. Electronic Citation, 2005.

[30] La Vision. www.lavision.de. Electronic Citation, 2006.

[31] PIVTEC GmvH. www.pivtec.com. Generic, 2006.

[32] PIV Challenge. http://www.pivchallenge.org. Electronic Citation, 2005.


APÉNDICE A

Características de los componentes de PIV.

A.1. Partículas

Las Tabla A.1 y la Tabla A.2 son unas listas de las partículas que pueden ser utilizadas para

visualizar líquidos y gases respectivamente [1].

Tipo Material Diámetro promedio en micras

Sólido Poliestireno 5 — 100

Aluminio 2 — 7

Esferas de vidrio 10 — 100

Gránulos con cubiertas sintéticas 10 — 500

Líquido Diferentes aceites 50 — 500

Gaseoso Burbujas de Oxígeno 50 — 1000

Tabla A.1: Partículas que pueden usarse en la visualización de flujos líquidos.

En el laboratorio de Velocimetría de ESIQIE estudian el comportamiento del acero líquido para

diseñar buzas o boquillas por las cuales sale el acero para formar láminas. Para sus experimentos,

utilizan agua y partículas de poliestireno de 5, 20 y 50 micras. En la Figura A.1 se muestra una foto

del experimento de ese laboratorio y la forma de la buza de acrílico utilizada cuyo tamaño es la tercera

parte de las buzas utilizadas con acero.

A.1. Partículas

Figura A.1: (a) Foto tomada en el 2004 de un experimento en el laboratorio de Velocimetría, (b) Buza

de acrílico por donde sale el agua.

Tipo Material Diámetro promedio en micras

Sólido Poliestireno 0.5 — 10

Aluminio 2 — 7

Magnesio 2 — 5

Micro-globos de vidrio 30 — 100

Gránulos con cubiertas sintéticas 10 — 50

Ftalato de dioctilo 1 — 10

Líquido Diferentes aceites 0.5 — 10

Humo < 1

Tabla A.2: Partículas que pueden usarse en la visualización de flujos gaseosos.


Apéndice A. Características de los componentes de PIV.

A.2. Láser

El término láser proviene del acrónimo de Light Amplification by Stimulated Emission of Radiation,

esto es, amplificación de la luz por emisión estimulada de radiación.

Los dispositivos que generan luz láser están formados básicamente por tres componentes:

El material láser es un compuesto químico que puede ser un gas, un semiconductor o un material

sólido,

la fuente de energía que alimenta al material láser por radiaciones electromagnéticas ó energía

química,

y el sistema óptico que es un arreglo de espejos que permiten la amplificación de la luz emitida

y una vez que esta alcanza cierta magnitud, la luz láser sale del dispositivo.

El láser es un dispositivo que emite un haz de luz con características muy especiales. Esta luz

del láser es monocromática, direccional, coherente e intensa; esta luz puede ser generada en pulsos o

de forma continua. Se produce cuando la fuente de energía excita a los electrones de los átomos del

material láser, estos pasan a un estado de mayor energía y espontáneamente pasan a otro estado de

menor energía liberando la energía sobrante en forma de un fotón. Los fotones emitidos, estimulan

a otros electrones para que éstos liberen más fotones con la misma energía y longitud de onda que

los demás ya generados. Esta luz se mantiene reflejada de un lado a otro por los espejos, para así

amplificarla y dejar salir un haz de luz. La luz del láser es altamente monocromática porque está

formada por ondas de una sola longitud o un rango de longitudes muy pequeño. La longitud de onda

de la luz depende de la cantidad de energía liberada, del sistema óptico, del material y la forma en

que son estimulados los átomos. Un láser puede emitir luz infrarroja, visible o ultravioleta. La luz es

coherente si sus ondas están en fase, es decir, que sus máximos, mínimos y puntos de inflexión están

alineados en los mismos puntos. La luz del láser es intensa y puede enfocarse con gran nitidez en

una superficie muy pequeña, y es altamente direccional o concentrada porque la luz se envía en una

dirección y ésta se dispersa muy poco. A diferencia de una luz común en la que la superficie iluminada,

aumenta a medida que se aleja de la fuente, con el láser esta relación se mantiene casi constante.

La forma más simple de crear un plano de luz es utilizar un lente cilíndrico en combinación con un

lente esférico. El lente cilíndrico expande el rayo de luz en una dirección al aumentar la divergencia del

rayo y generar un plano. Después se puede utilizar un lente esférico para controlar el espesor de este

plano. En la referencia [1] se puede encontrar más información de la óptica necesaria para generar el


A.2. Láser

plano de luz láser y en la Figura A.2 se muestra un ejemplo de cómo pueden ser acomodados los lentes

para generar el plano de luz.

Figura A.2: Óptica para generar el plano de luz láser.

La Tabla A.3 muestra algunas propiedades de los sistemas Nd:YAG PIV láser comerciales.

Tasa de Repetición 10-100 Hz

Energía del pulso por cada uno de dos pulsos 150-450 mJ

Duración del pulso 3-12 ns

Temperaturas de funcionamiento 15o - 35oC

Requerimientos de poder 220 — 240 V, 50 Hz

Tabla A.3: Propiedades de un sistema de láser Nd:YAG para PIV.

Existen varias empresas que venden el sistema de láser solo o el equipo completo con el software y

hardware necesario para la técnica de PIV [28], [13] y [29]. El equipo tiene lo último en tecnología, con el

fin de mejorar su rendimiento. Las computadoras cuentan con los procesadores más veloces, tecnología

multi-procesador, así como gran capacidad de disco duro y memoria RAM. Además el software incluido

puede adecuarse a una gran variedad de aplicaciones.



A.3. Cámara digital

El CCD es un sensor electrónico que convierte luz en una carga eléctrica y esta formado por millones

de sensores fotosensibles de tamaño microscópico, acomodados en una cuadrícula. Los sensores del CCD

registran la intensidad del brillo de un color (rojo, verde o azul). Dependiendo de esta intensidad, los

sensores reaccionan creando una carga eléctrica y a partir de su posición en la cuadrícula, se construye

la imagen. Las imágenes digitales están formadas por una cuadrícula de puntos llamados píxeles. La

resolución es medida en mega píxeles (por cada millón de sensores). Desde el punto de vista visual,

una mayor resolución en una imagen es mejor, pues se tiene mejor detalle. Estas imágenes pueden ser

guardadas en la memoria de la cámara o en algunos casos en la memoria de la computadora. En el caso

de PIV es preferible utilizar cámaras que transfieran las imágenes de forma directa a la computadora,

con el fin de tener mayor capacidad de almacenamiento. A continuación se muestra en la Tabla A.4

algunas cámaras digitales para PIV de distintas empresas.

A.4. Sistemas de PIV

A continuación se muestra una tabla con los principales paquetes comerciales de PIV, de distintas

empresas.


A.4. Sistemas de PIV

Marca PCO La Vision Redlake Mo-

tionPro

TSI

Modelo 1200 HS Imagerpro HS HS-3 PIVCAM HS-

500

Máxima Re-

solución

1280 x 1024

píxeles

1280 x 1024

píxeles

1280 x 1024

píxeles

1280 x 1024

píxeles

Captura de

imágenes por

segundo

636 fps < 635 fps 1000 fps 500 fps

Tipo de sen-

sor

CCD CCD CCD CCD

Rango Di-

námico

10 bit 10 bit No Disponible 8 bit

Almacenamiento4 Gb 4 Gb 4 Gb 4 Gb

Tabla A.4: Cámaras digitales para PIV



Dantec Dynamics [28]FlowMap PIV

2DFlowMap PIV 3D

Oxford Lasers [13]

VisiVector 2D

Max VisiVector 2D

Micro VisiVector

VisiVector 3D

Max VisiVector 3D

LaVision [30] FlowMaster HS (2D,3D, Micro, PLIF)

TSI Inc. [29]

PIV System High Frame Rate

PIV System PLIF

PIV System 3D

PIV System MicroPIV

PIV System PIF/PLIF

PIVTEC [31] PIVview (2D,3D)

Tabla A.5: Sistemas PIV comerciales


APÉNDICE B

Diagramas

B.1. Diagramas de flujo

Figura B.1: Algoritmo de pre-procesamiento.

B.1. Diagramas de flujo

Figura B.2: Algoritmo de procesamiento.


Apéndice B. Diagramas

Figura B.3: Algoritmo de post-procesamiento.


B.2. Diagrama de estados

B.2. Diagrama de estados

B.3. Pantallas de la interfaz gráfica



Figura B.4: Diagrama de estados de pre-procesamiento.



Figura B.5: Diagrama de estados de procesamiento.



Figura B.6: Diagrama de estados de post-procesamiento.

Figura B.7: Diagrama de estados del despliegue de resultados.



Figura B.8: Diseño de pantalla de datos generales.

Figura B.9: Diseño de pantalla de datos de pre-procesamiento.



Figura B.10: Diseño de pantalla de datos de procesamiento.

Figura B.11: Diseño de pantalla de datos de post-procesamiento.


APÉNDICE C

Código

C.1. Programa principal.

% PREPROCESAMIENTO

tam=size(mascara);

if(tam(2)==4)

f_original = f_original(mascara(3):mascara(4),mascara(1):mascara(2));

g_original = g_original(mascara(3):mascara(4),mascara(1):mascara(2));

if(ruido)

h_original = h_original(mascara(3):mascara(4),mascara(1):mascara(2));

k_original = k_original(mascara(3):mascara(4),mascara(1):mascara(2));

end

end

f = funcionPreprocesamiento(f_original, ecualiza);

g = funcionPreprocesamiento(g_original, ecualiza);

if(ruido)

mascaraRuido = funcionEliminaRuido(f_original,g_original,h_original,k_original);

else

mascaraRuido = 0;

end

f = f - mascaraRuido;

g = g - mascaraRuido;

Tam = size(f_original);

height = Tam(1);

width = Tam(2);

C.1. Programa principal.

% PROCESAMIENTO

[CX, CY, U, V,window_X,window_Y] =

funcionCorrelacion( f,g,height,width,window_X,window_Y,aumentoX,aumentoY, traslapeX,

traslapeY,subpixel_interpolacion,iteraciones,validacion,sustitucion);

% POSTPROCESAMIENTO

if(validacion)

NH = ceil(height/traslapeX);

NW = ceil(width/traslapeY);

Num_ventanas = NH*NW;

if(NH>=3 && NW>=3)

[U,V] = funcionValidacion(CX, CY, U, V, sustitucion);

end

end

Tam = size(CX);

fil_vect = Tam(1);

col_vect = Tam(2);

tot = fil_vect*col_vect;

MatrizVectores = zeros(tot,4);

cont = 1;

for j=1:col_vect

for i=1:fil_vect

MatrizVectores(cont,1) = CX(i,j);

MatrizVectores(cont,2) = CY(i,j);

MatrizVectores(cont,3) = U(i,j);

MatrizVectores(cont,4) = V(i,j);

cont = cont + 1;

end

end

% RESULTADOS

if(max(size(mascara))==1)

des_x=0;

des_y=0;

else

des_x=mascara(1);

des_y=mascara(3);


Apéndice C. Código

end

if(guardar_datos|grafica)

f = imread(imagen1);

imprime =

funcionImprime(f,g,MatrizVectores,guardar_datos,grafica,window_X,window_Y,des_x,des_y);

end

C.2. Pre-procesamiento

function imagen = funcionPreprocesamiento(imagen_original,ecualiza)

if(ecualiza)

imagen_filt= medfilt2(imagen_original,[3 3]);

elemento = strel(’square’,4);

imagen_filt = imdilate(imagen_filt,elemento);

imagen_original = imagen_original-imagen_filt;

imagen_original = imadjust(imagen_original,stretchlim(imagen_original),[]);

end

imagen = imagen_original;

function mascara = funcionEliminaRuido(f_original,g_original,h_original,k_original)

tam = size(f_original);

tamy = tam(1);

tamx = tam(2);

mascaraAnd = f_original;

mascaraMean4 = f_original;

for i=1:tamx-1

for j=1:tamy-1

w = (bitand(uint8(f_original(j,i)), uint8(g_original(j,i))));

w = (bitand(uint8(h_original(j,i)), uint8(w)));

w = (bitand(uint8(k_original(j,i)), uint8(w)));

if(w >= 128)

w = 255;

else

w = 0;

end


C.3. Procesamiento

mascaraAnd(j,i) = w;

end

end

mascara = uint8(mascaraAnd);

C.3. Procesamiento

function [CX, CY, U, V,window_X,window_Y,archivo] =

funcionCorrelacion(f,g,height,width,window_X,window_Y,aumentoX,aumentoY, traslapeX,

traslapeY,subpixel_interpolacion,iteraciones,validacion,sustitucion, archivo)

mid_Wx = window_X/2;

mid_Wy = window_Y/2;

window_X = 2*mid_Wx;

window_Y = 2*mid_Wy;

tx = 0;

ty = 0;

if(aumentoX>0)

tx = floor((aumentoX/100)*window_X);

end

if(aumentoY>0)

ty = floor((aumentoY/100)*window_Y);

end

NH = ceil(height/(window_Y*traslapeY));

NW = ceil(width/(window_X*traslapeX));

z=NH*NW;

c_x=zeros(z,1);

c_y=zeros(z,1);

u=zeros(z,1);

v=zeros(z,1);

uS=zeros(z,1);

vS=zeros(z,1);

indice=1;

iteraciones=iteraciones-1;

cont1=0;



% Primera Iteraci\U{f3}n para calcular la correlaci\U{f3}n

for i = 1:(traslapeY*window_Y):height

lim_iy = (i+window_Y)-1;

if lim_iy > height

lim_iy=height;

end

cont1 = cont1 + 1;

for j = 1:(traslapeX*window_X):width

lim_jx = (j+window_X)-1;

if lim_jx > width

lim_jx=width;

end

k = i-ty;

lim_ky = lim_iy + ty;

l = j-tx;

lim_lx = lim_jx + tx;

if k < 1

k=1;

end

if l < 1

l=1;

end

if lim_ky > height

lim_ky = height;

end

if lim_lx > width

lim_lx = width;

end

zona1 = f(i:lim_iy,j:lim_jx,1);

zona2 = g(k:lim_ky,l:lim_lx,1);

Z1 = size(zona1);

Z2 = size(zona2);

if(Z1==Z2)

window = max(Z2(1), Z2(2));

T = 2*window;


C.3. Procesamiento

mat_corr = zeros(T,T);

zona1 = zona1 - mean2(zona1);


zona2 = zona2(end:-1:1,end:-1:1);

fft_1=fft2(double(zona1),T,T);


mat_corr = real(ifft2(fft_1.*fft_2));

else

mat_corr = normxcorr2(double(zona1),double(zona2));

end

[valor,cord_x]= max(max(mat_corr,[],1));

[valor,cord_y]= max(max(mat_corr,[],2));

c_w(indice) = cord_x;

c_z(indice) = cord_y;

c_v(indice) = valor;

c_y(indice) = i+mid_Wy;

c_x(indice) = j+mid_Wx;

if(valor~=0)

u(indice) = Z2(2)-cord_x;

v(indice) = Z2(1)-cord_y;

if(subpixel_interpolacion>0)

picos = funcionSubpixel(mat_corr,cord_x,cord_y,subpixel_interpolacion);

uS(indice) = Z2(2)-picos(1,1);

vS(indice) = Z2(1)-picos(1,2);

else

uS(indice) = 0;

vS(indice) = 0;

end

else

u(indice) = 0;

v(indice) = 0;

end

indice=indice+1;

end

end



%%-------Resultados----------%%

if( exist(’c_v’))

H=size(c_x,1);

MatrizResultados = zeros(H(1),6);

MatrizResultados(:,1) = c_x(:,1);

MatrizResultados(:,2) = c_y(:,1);

MatrizResultados(:,3) = u(:,1);

MatrizResultados(:,4) = v(:,1);

fil_vect = cont1;

col_vect = H/cont1;

CX = zeros(fil_vect, col_vect);

CY = zeros(fil_vect, col_vect);

U = zeros(fil_vect, col_vect);

V = zeros(fil_vect, col_vect);


MatrizResultados(:,5) = uS;

MatrizResultados(:,6) = vS;

else

MatrizResultados = MatrizResultados(:,1:4);

End

if (iteraciones>=1)

[MatrizResultados,window_X,window_Y,Num_vectores_Nuevos,NW] =

funcionIterativa(f,g,height,width,window_X,window_Y,MatrizResultados,

subpixel_interpolacion,iteraciones,validacion,sustitucion,fil_vect,

col_vect,traslapeX, traslapeY, archivo);

fil_vect = NW;

col_vect = Num_vectores_Nuevos/NW;

c_x = MatrizResultados(:,1);

c_y = MatrizResultados(:,2);

u = MatrizResultados(:,3);

v = MatrizResultados(:,4);


uS = MatrizResultados(:,5);

vS = MatrizResultados(:,6);

end


C.3. Procesamiento

end

cont=1;

for i=1:fil_vect

for j=1:col_vect

CX(i,j) = c_x(cont);

CY(i,j) = c_y(cont);


U(i,j) = uS(cont);

V(i,j) = vS(cont);

else

U(i,j) = u(cont);

V(i,j) = v(cont);

end

cont = cont + 1;

end

end

end

function picos = funcionSubpixel(mat_corr,cord_x,cord_y,tipo)

picos=zeros(1,2);

T = size(mat_corr);

lim_y = T(1);

lim_x = T(2);

if cord_x == 1

izquierda = 1;

else

izquierda = cord_x-1;

end

if cord_x == lim_x

derecha = lim_x;

else

derecha = cord_x+1;

end

if cord_y == 1

arriba = 1;



else

arriba = cord_y-1;

end

if cord_y == lim_y

abajo = lim_y;

else

abajo = cord_y+1;

end

% INTERPOLACION SUBPIXEL

if (tipo == 1)%peak centroid

x0PC = (izquierda)*mat_corr(cord_y,izquierda) +

(cord_x)*mat_corr(cord_y,cord_x)+(derecha)*mat_corr(cord_y,derecha);

x0PC = x0PC/(mat_corr(cord_y,izquierda) + mat_corr(cord_y,cord_x)

+mat_corr(cord_y,derecha));

y0PC = (arriba)*mat_corr(arriba,cord_x) +

(cord_y)*mat_corr(cord_y,cord_x)+(abajo)*mat_corr(abajo,cord_x);

y0PC = y0PC/(mat_corr(arriba,cord_x) + mat_corr(cord_y,cord_x)+mat_corr(abajo,cord_x));

picos(1,1) = real(x0PC);

picos(1,2) = real(y0PC);

elseif (tipo == 2)%parabolic peak fit

x0PPF = mat_corr(cord_y,izquierda)- mat_corr(cord_y,derecha);

x0PPF = x0PPF/(2*mat_corr(cord_y,izquierda) - 4*mat_corr(cord_y,cord_x)+

2*mat_corr(cord_y,derecha));

x0PPF = x0PPF + cord_x;

y0PPF = mat_corr(arriba,cord_x)-mat_corr(abajo,cord_x);

y0PPF = y0PPF/(2*mat_corr(arriba,cord_x) - 4*mat_corr(cord_y,cord_x)+

2*mat_corr(abajo,cord_x));

y0PPF = y0PPF + cord_y;

picos(1,1) = real(x0PPF);

picos(1,2) = real(y0PPF);

elseif (tipo == 3)%gaussian peak fit

x0GPF = log(mat_corr(cord_y,izquierda))-log(mat_corr(cord_y,derecha));

x0GPF = x0GPF/(2*log(mat_corr(cord_y,izquierda)) - 4*log(mat_corr(cord_y,cord_x))+

2*log(mat_corr(cord_y,derecha)));

x0GPF = x0GPF + cord_x;


C.3. Procesamiento

y0GPF = log(mat_corr(arriba,cord_x))-log(mat_corr(abajo,cord_x));

y0GPF = y0GPF/(2*log(mat_corr(arriba,cord_x)) - 4*log(mat_corr(cord_y,cord_x)) +

2*log(mat_corr(abajo,cord_x)));

y0GPF = y0GPF + cord_y;

picos(1,1) = real(x0GPF);

picos(1,2) = real(y0GPF);

end

function [MatrizResultados,window_X,window_Y,Num_vectores_Nuevos,TW] =

funcionIterativa( f,g,height,width,window_X,window_Y,MatrizResultados,

subpixel_interpolacion, iteraciones,validacion,sustitucion,

fil_vect,col_vect,traslapeX, traslapeY, archivo)

while (iteraciones>0)

if(validacion)

c_x(:,1) = MatrizResultados(:,1);

c_y(:,1) = MatrizResultados(:,2);


u(:,1) = MatrizResultados(:,5);

v(:,1) = MatrizResultados(:,6);

else

u(:,1) = MatrizResultados(:,3);

v(:,1) = MatrizResultados(:,4);

end

CX = zeros(fil_vect, col_vect);

CY = zeros(fil_vect, col_vect);

U = zeros(fil_vect, col_vect);

V = zeros(fil_vect, col_vect);

cont=1;

for j=1:col_vect

for i=1:fil_vect

CX(i,j) = c_x(cont);

CY(i,j) = c_y(cont);

U(i,j) = u(cont);

V(i,j) = v(cont);

cont = cont + 1;



end

end

if(col_vect>=3 && fil_vect>=3)

[U,V] = funcionValidacion(CX, CY, U, V, sustitucion);

cont = 1;

for j=1:col_vect

for i=1:fil_vect

MatrizResultados(cont,3) = U(i,j);

MatrizResultados(cont,4) = V(i,j);

cont = cont + 1;

end

end

end

end

MatrizVectores = MatrizResultados;

T=size(MatrizResultados);

Num_vectores = T(1);

Desplazamiento = MatrizResultados(:,1:4);

c_xOriginal = Desplazamiento(:,1);

c_yOriginal = Desplazamiento(:,2);

NW = floor(width/window_X);

window_X = floor(window_X*.5);

window_Y = floor(window_Y*.5);

TW = ceil(width/(window_X*traslapeX));

TH = ceil(height/(window_Y*traslapeY));

Num_vectores_Nuevos = TW*TH;

c_x = zeros(Num_vectores_Nuevos,1);

c_y = zeros(Num_vectores_Nuevos,1);

u = zeros(Num_vectores_Nuevos,1);

v = zeros(Num_vectores_Nuevos,1);

uS = zeros(Num_vectores_Nuevos,1);

vS = zeros(Num_vectores_Nuevos,1);

c_x(1:Num_vectores,1) = Desplazamiento(:,1);

c_y(1:Num_vectores,1) = Desplazamiento(:,2);

u(1:Num_vectores,1) = Desplazamiento(:,3);


C.3. Procesamiento

v(1:Num_vectores,1) = Desplazamiento(:,4);

indice=1;

indiceNuevo=0;

banderay = 0;

banderax = 0;

incrementex=0;

incrementey=0;

for i=1:ceil(traslapeY*window_Y):height

for j=1:ceil(traslapeX*window_X) :width

indiceNuevo = indiceNuevo+1;

lim_iy = (i+ceil(traslapeY*window_Y))-1;

if lim_iy > height

lim_iy=height;

end

lim_jx = (j+ceil(traslapeX*window_X))-1;

if lim_jx > width

lim_jx=width;

end

k1 = v(indice);

l1 = u(indice);

k = i + ceil(k1);

lim_ky = ceil(k+(traslapeY*window_Y)-1);

l = j + ceil(l1);

lim_lx = ceil(l+(traslapeX*window_X)-1);

if k < 1

k=1;

end

if l < 1

l=1;

end

if lim_ky > height

lim_ky = height;

end

if(lim_ky<1)

lim_ky=1;



end

if lim_lx > width

lim_lx = width;

end

if(lim_lx<1)

lim_lx=1;

end

uOriginal = u(indice);

vOriginal = v(indice);

zona1 = f(i:lim_iy,j:lim_jx,1);

zona2 = g(k:lim_ky,l:lim_lx,1);

Z1 = size(zona1);

Z2 = size(zona2);

if(Z1(1)~=Z2(1)||Z1(2)~=Z2(2))

mH = max(Z1(1),Z2(1));

mW = max(Z1(2),Z2(2));

z1 = f(:,:,1);

z2 = f(:,:,1);

z1(1:Z1(1),1:Z1(2))=zona1;

z1(Z1(1)+1:end,:,1)=0;

z1(:,Z1(2)+1:end,1)=0;

z2(1:Z2(1),1:Z2(2))=zona2;

z2(Z2(1)+1:end,:,1)=0;

z2(:,Z2(2)+1:end,1)=0;

zona1=z1(1:mH,1:mW);

zona2=z2(1:mH,1:mW);

Z1 = size(zona1);

Z2 = size(zona2);

end

window = max(Z2(1), Z2(2));

T = 2*window;

mat_corr = zeros(T,T);



zona2 = zona2(end:-1:1,end:-1:1);


C.3. Procesamiento



mat_corr = real(ifft2(fft_1.*fft_2));

[valor,cord_x]= max(max(mat_corr,[],1));

[valor,cord_y]= max(max(mat_corr,[],2));

c_w(indiceNuevo) = cord_x;

c_z(indiceNuevo) = cord_y;

c_v(indiceNuevo) = valor;

if(valor~=0)

c_y(indiceNuevo) = floor(i+window_Y/2);

c_x(indiceNuevo) = floor(j+window_X/2);

u(indiceNuevo) = Z2(2)-cord_x+uOriginal;

v(indiceNuevo) = Z2(1)-cord_y+vOriginal;


picos = funcionSubpixel(mat_corr,cord_x,cord_y,

subpixel_interpolacion);

uS(indiceNuevo) = Z2(2)-picos(1,1)+uOriginal;

vS(indiceNuevo) = Z2(1)-picos(1,2)+vOriginal;

end

else

c_y(indiceNuevo) = 0;

c_x(indiceNuevo) = 0;

uS(indiceNuevo) = 0;

vS(indiceNuevo) = 0;

end

if(j>width - 2*window_X)

[banderax, indiceR, incrementex] = decideX2(banderax,indice);

else

[banderax, indice, incrementex] = decideX2(banderax,indice);

end %if else de caso especial si cabeX

end %for j (x)

banderax = 0;

[banderay,indice,incrementey] = decideY2(banderay, indice);

if(banderay==0)

indice = indice-(NW-1);



elseif(i>height - 2*window_Y)

if(incrementey)

indice=indice-1;

end

indice = indice-(NW-1);

end%if else caso especial cabeY

end %for i (y)

%%-------Resultados----------%%

if( exist(’c_v’))

Desplazamiento = zeros(Num_vectores_Nuevos,6);

MatrizResultados = zeros(Num_vectores_Nuevos,6);

Desplazamiento(:,1) = c_x;

Desplazamiento(:,2) = c_y;

Desplazamiento(:,3) = u;

Desplazamiento(:,4) = v;


Desplazamiento(:,5) = uS;

Desplazamiento(:,6) = vS;

else

Desplazamiento= Desplazamiento(:,1:4);

end

MatrizResultados = Desplazamiento;

iteraciones=iteraciones-1;

end

end

function[resbx,indiceR,incremente] = decideX2(banderax, indice)

switch(banderax)

case 0

resbx = banderax + 1;

indiceR = indice;

incremente = 0;

case 1

resbx = 0;

indiceR = indice + 1;


C.3. Procesamiento

incremente = 1;

end

function [resby, indiceR, incremente] = decideY2(banderay, indice)

switch(banderay)

case 0

resby = banderay + 1;

indiceR = indice;

incremente = 0;

case 1

resby = 0;

indiceR = indice + 1;

incremente = 1;

end

function vecindad = funcionVecinos(U,V,i,j,tipo)

Tam = size(U);

NH = Tam(1);

NW = Tam(2);

if (tipo == 1)


vecindad = zeros(12,1);

promedioU = 0;

promedioV = 0;

diferenciaU = 0;

diferenciaV = 0;

valorU=0;

valorV=0;

valor1U=0;

valor1V=0;

valor2U=0;

valor2V=0;

nuevoU = 0;

nuevoV = 0;

total_vectores=0;



originalU = U(j,i);

originalV = V(j,i);

if(i<=NW && j==1) %Vectores primera fila

if(mod(i,NW)==0) %Vector primera fila, ultima columna

num_vectores=0;

for x=-1:1:0

for y=0:1:1

vectorU = U(j+y,i+x);

vectorV = V(j+y,i+x);

if(vectorU~=0 || vectorV~=0)

promedioU = promedioU + vectorU;

promedioV = promedioV + vectorV;

valor1U = valor1U + abs(vectorU-originalU);

valor1V = valor1V + abs(vectorV-originalV);

num_vectores = num_vectores+1;

end

end

end

total_vectores = 3;

num_vectores = num_vectores-1;

valor2U = promedioU-originalU;

valor2V = promedioV-originalV;

if(num_vectores>0)

nuevoU = valor2U/num_vectores;

nuevoV = valor2V/num_vectores;

else

nuevoU = 0;

nuevoV = 0;

end

elseif (mod(i,NW)==1) %Vector primera fila, primera columna

num_vectores=0;

for x=0:1:1

for y=0:1:1




C.3. Procesamiento







end

end

end

total_vectores = 3;




if(num_vectores>0)



else

nuevoU = 0;

nuevoV = 0;

end

else %Vectores primera fila

num_vectores=0;

for x=-1:1:1

for y=0:1:1









end

end

end



total_vectores = 5;




if(num_vectores>0)



else

nuevoU = 0;

nuevoV = 0;

end

end

elseif (j==NH) %Vectores ultima fila

if(mod(i,NW)==0) %Vector ultima fila, ultima columna

num_vectores=0;

for x=-1:1:0

for y=-1:1:0









end

end

end

total_vectores = 3;




if(num_vectores>0)




C.3. Procesamiento

else

nuevoU = 0;

nuevoV = 0;

end

elseif(mod(i,NW)==1) %Vector ultima fila, primera columna

num_vectores=0;

for x=0:1:1

for y=-1:1:0









end

end

end

total_vectores = 3;




if(num_vectores>0)



else

nuevoU = 0;

nuevoV = 0;

end

else %Vectores ultima fila

num_vectores=0;

for x=-1:1:1

for y=-1:1:0











end

end

end

total_vectores = 5;




if(num_vectores>0)



else

nuevoU = 0;

nuevoV = 0;

end

end

elseif (mod(i,NW)==0) %Vectores ultima columna

num_vectores=0;

for x=-1:1:0

for y=-1:1:1









end


C.3. Procesamiento

end

end

total_vectores = 5;




if(num_vectores>0)



else

nuevoU = 0;

nuevoV = 0;

end

elseif(mod(i,NW)==1) %Vectores primera columna

num_vectores=0;

for x=0:1:1

for y=-1:1:1









end

end

end

total_vectores = 5;




if(num_vectores>0)





else

nuevoU = 0;

nuevoV = 0;

end

else %Vectores centrales con 8 vecinos

num_vectores=0;

for x=-1:1:1

for y=-1:1:1









end

end

end

total_vectores = 8;




if(num_vectores>0)



else

nuevoU = 0;

nuevoV = 0;

end

end

diferenciaU = abs(originalU-valor2U);

diferenciaV = abs(originalV-valor2V);

valorU = abs(valor1U);

valorV = abs(valor1V);


C.3. Procesamiento

if(valor2U>0)

valorU = abs(valor1U)/abs(valor2U);

end

if(valor2V>0)

valorV = abs(valor1V)/abs(valor2V);

end

vecindad(1,1) = originalU;

vecindad(2,1) = originalV;

vecindad(3,1) = total_vectores;

vecindad(4,1) = num_vectores;

vecindad(5,1) = nuevoU;

vecindad(6,1) = nuevoV;

vecindad(7,1) = promedioU;

vecindad(8,1) = promedioV;

vecindad(9,1) = diferenciaU;

vecindad(10,1) = diferenciaV;

vecindad(11,1) = valorU;

vecindad(12,1) = valorV;

elseif (tipo == 2) % Puntos = funcionVecinos(U,V,i,j,2);

contador = 1;

vecindad = zeros(6,2);

for k =-1:1

for l=-1:1

ik = i + k;

jl = j + l;

if(ik>=1 && ik<=NW && jl>=1 && jl<=NH)

if(U(jl,ik)~=0 && V(jl,ik)~=0 && contador<=6)

vecindad(contador,1) = jl;

vecindad(contador,2) = ik;

contador = contador +1;

end

end

end

end

if(contador<3)



for k =-2:2

for l=-2:2

ik = i + k;

jl = j + l;

if(ik>=1 && ik<=NW && jl>=1 && jl<=NH)





end

end

end

end

end

if(contador<3)

for k =-3:3

for l=-3:3

ik = i + k;

jl = j + l;

if(ik>=1 && ik<=NW && jl>=1 && jl<=NH)





end

end

end

end

end

vecindad = vecindad(1:contador-1,:);

end


C.4. Post-procesamiento


function [U,V] = funcionValidacion(CX, CY, U, V, sustitucion)

Tam = size(U);

NH = Tam(1);

NW = Tam(2);


Errores = zeros(Num_ventanas,2);

indiceErrores = 0;

Diferencias_U = zeros(Num_ventanas,1);

Diferencias_V = zeros(Num_ventanas,1);

%VALIDACION

indice=1;

for i=1:NW

for j=1:NH

vecindad = funcionVecinos(U,V,i,j,1);

Diferencias_U(indice,1) = vecindad(9,1);

Diferencias_V(indice,1) = vecindad(10,1);

Valores_U(indice,1) = vecindad(1,1);

Valores_V(indice,1) = vecindad(2,1);

indice = indice + 1;

end

end

desDU = std(Diferencias_U);

desDV = std(Diferencias_V);

meanU = mean(Valores_U);

meanV = mean(Valores_V);

desU = std(Valores_U);

desV = std(Valores_V);

for j=1:NH

for i=1:NW


diferenciaU = vecindad(9,1);

diferenciaV = vecindad(10,1);

if( diferenciaU > desDU*1.65 || diferenciaV > desDV*1.65



|| abs(U(j,i)-meanU)>=desU

|| abs(V(j,i)-meanV)>=desV)

indiceErrores = indiceErrores+1;

U(j,i) = 0;

V(j,i) = 0;

Errores(indiceErrores,1) = j;

Errores(indiceErrores,2) = i;

end

end

end

Errores = Errores(1:indiceErrores,:);

numero_errores = indiceErrores;

% SUSTITUCION

numero_anterior = 0;

bandera=-1;

if(sustitucion==1)

while(numero_errores>0)

if(numero_anterior==numero_errores)

if(bandera==5)

break;

else

bandera=bandera+1;

end

end

for k = 1:indiceErrores

j = Errores(k,1);

i = Errores(k,2);

if(j>0 && i>0)


total_vectores = vecindad(3,1);

num_vectores = vecindad(4,1);

if(num_vectores/total_vectores >=.35)

nuevoU = vecindad(5,1);

nuevoV = vecindad(6,1);

U(j,i) = nuevoU;



V(j,i) = nuevoV;

numero_errores = numero_errores-1;

bandera=0;

Errores(k,:) = 0;

end

end

end

numero_anterior = numero_errores;

end

elseif(sustitucion==2)

while(numero_errores>0)

if(numero_anterior==numero_errores)

if(bandera==5)

break;

else

bandera=bandera+1;

end

end

for k = 1:indiceErrores

j = Errores(k,1);

i = Errores(k,2);

if(j>0 && i>0)

Puntos = funcionVecinos(U,V,i,j,2);

tam = size(Puntos);

if(tam(1)>=3)

[xU, xV] = funcionCoeficientes(CX, CY, U, V,Puntos);

if( max(isnan(xU))==0 && max(isnan(xV)==0))

if(tam(1)>=6)

u_aprox = xU(1,1)+ xU(2,1)*CX(j,i) + xU(3,1)*CY(j,i) +

xU(4,1)*CX(j,i)*CX(j,i) +

xU(5,1)*CY(j,i)*CY(j,i) + xU(6,1)*CX(j,i)*CY(j,i);

v_aprox = xV(1,1)+ xV(2,1)*CX(j,i) + xV(3,1)*CY(j,i) +

xV(4,1)*CX(j,i)*CX(j,i) +

xV(5,1)*CY(j,i)*CY(j,i) + xV(6,1)*CX(j,i)*CY(j,i);

else



u_aprox = xU(1,1)+ xU(2,1)*CX(j,i)+ xU(3,1)*CY(j,i);

v_aprox = xV(1,1)+ xV(2,1)*CX(j,i)+ xV(3,1)*CY(j,i);

end

U(j,i) = u_aprox;

V(j,i) = v_aprox;

numero_errores = numero_errores-1;

bandera=0;

Errores(k,:)=0;

end

end

end

end

numero_anterior = numero_errores;

end

end

C.5. Resultados

function imprime = funcionImprime(f,g,MatrizResultados,guardar_datos,grafica,window_X,

window_Y,des_x,des_y)

imprime=1;

Tam = size(f);

height = Tam(1);

width = Tam(2);

% GRAFICAR los vectores

if(grafica>0)

c_x = MatrizResultados(:,1);

c_y = MatrizResultados(:,2);

u = MatrizResultados(:,3);

v = MatrizResultados(:,4);

c_x = c_x + des_x;

c_y = c_y + des_y;

figure;

quiver(c_x’,c_y’,u’,v’,’b’), title(’Resultados’,’FontWeight’,’bold’);


C.5. Resultados

axis([1 width 1 height])

xlabel(’EJE x’,’FontWeight’,’bold’)

ylabel(’EJE y’,’FontWeight’,’bold’)

set(gca,’YDir’,’reverse’)

figure;

im1 = double( f’ );

im2 = double( g’ );

x = 1:width;

y = 1:height;

image_max = max(max(im1));

im = 75*im1/image_max;

image(x, y, im’);

colormap(gray)

hold on

quiver( c_x’, c_y’, u’, v’, ’g’ ), title(’Resultados’,’FontWeight’,’bold’);

hold off

xlabel(’EJE x’,’FontWeight’,’bold’)

ylabel(’EJE y’,’FontWeight’,’bold’)

end% GUARDAR los vectores


APÉNDICE D

Instalación del sistema

Para instalar el sistema, basta descomprimir el archivo .ZIP incluido en alguna carpeta.

D.1. Matlab

Para utilizar el sistema aquí propuesto es necesario tener la versión 6.0 (Release 12 ) o alguna

posterior.

D.2. Requerimientos

Los requerimientos para ejecutar la versión 6.0 de Matlab incluyen:

D.2.1. General

Lector de CD-ROM drive (para la instalación)

Netscape Navigator 4.0 o superior o Microsoft Internet Explorer 4.0 o superior.

Adobe Acrobat Reader 3.0 o superior para ver o imprimir la documentación en línea de MATLAB

en formato PDF.

Algunos tipos de licencia requieren una licencia en un servidor corriendo FLEXlm 6.1, que lo

provee el instalador de MathWorks.

Se requiere TCP/IP en todas las plataformas si se usa una licencia de servidor.

D.2. Requerimientos

D.2.2. Específicos de Windows.

Microsoft Windows 95, Windows 98 (original and Second Edition), Windows Millennium Edition,

Windows NT 4.0 (with Service Pack 5 for Y2K compliancy or Service Pack 6a) or Windows 2000

Pentium, Pentium Pro, Pentium II, Pentium III, Pentium IV* or AMD Athlon basados en ar-

quitectura PC.

64 MB RAM mínimo, se recomienda 128 MB RAM.

El espacio en disco duro varía dependiendo en el tamaño de la partición y la instalación de los

archivos de ayuda. El instalador de MATLAB informará del espacio requerido para una partición

en particular.

Adaptador gráfico de 8-bit y un monitor (que despliegue 256 colores simultáneamente).

Se recomiendan otros periféricos que incluyen:

Una tarjeta aceleradora de gráficos compatible con MS Windows.

Una impresora compatible con MS Windows.

Una tarjeta de sonido compatible con MS Windows.

Microsoft Word 7.0 (Office 95), 8.0 (Office 97) o Office 2000 se requiere para correr el Matlab

Notebook.

Se requiere uno de los siguientes compiladores para construir los MEX-files propios:

• Compaq Visual Fortran 5.0 ó 6.1

• Microsoft Visual C/C++ versión 5.0 ó 6.0

• Borland C/C++ versión 5.0, 5.02

• Borland C++Builder versión 3.0, 4.0, ó 5.0

• WATCOM versión 10.6 ó 11

• LCC 2.4 (incluido con MATLAB)

D.2.3. Específicos de UNIX/Linux

90 MB libres en disco duro para MATLAB solamente (215 MB para incluir los archivos de ayuda

de MATLAB)


Apéndice D. Instalación del sistema

Vendedor Hardware S/O Versión

Sun SPARC, ULTRA Solaris 2.6, 2.7, 2.8

HP PA-RISC 1.1 HP-UX 10.20

HP PA-RISC 2.0 HP-UX 11.0

SGI R5000, R8000, R10000, R12000 IRIX, IRIX 64 6.5.8 - 6.5.12

IBM RS/6000 AIX 4.3.3, 5.1

Compaq Alpha Tru64 UNIX 4.0 (mínimo), 5.0, 5.1

Linux Pentium, Pentium Pro, II, III, IV*,

AMD Athlon

ix86 - Linux 2.2.x, 2.4.x** kernel

64 MB RAM mínimo, se recomiendan 128 MB RAM.

64 MB de espacio en swap (recomendado)

Para Linux, MATLAB 6.0 solo correrá bajo un ambiente glibc (libc6). El ambiente libc5 ya no es

soportado.

D.3. Formatos para imágenes de entrada

El formato de las imágenes debe ser alguno soportado por Matlab, que incluyen:

: Formatos de imágenes soportados por MATLAB.

D.4. Funciones utilizadas incluidas en Matlab.

A continuación se muestran las funciones de Matlab implementadas por Mathworks. Para mayor

información de las funciones se recomienda consultar el manual o la ayuda de MATLAB.

Filtro de Mediana (2 Dimensiones). medfilt2(imagen,[n m])


D.4. Funciones utilizadas incluidas en Matlab.

Formato Nombre Completo

’bmp’ Windows Bitmap (BMP)

’cur’ Windows Cursor resources (CUR)

’gif’ Graphics Interchange Format (GIF)

’hdf’ Hierarchical Data Format (HDF)

’ico’ Windows Icon resources (ICO)

’jpg’ or ’jpeg’ Joint Photographic Experts Group (JPEG)

’pbm’ Portable Bitmap (PBM)

’pcx’ Windows Paintbrush (PCX)

’pgm’ Portable Graymap (PGM)

’png’ Portable Network Graphics (PNG)

’pnm’ Portable Anymap (PNM)

’ppm’ Portable Pixmap (PPM)

’ras’ Sun Raster (RAS)

’tif’ or ’tiff’ Tagged Image File Format (TIFF)

’xwd’ X Windows Dump (XWD)

Tabla D.1: Formatos de imágenes soportados por MATLAB.


Apéndice D. Instalación del sistema

El filtro de mediana es una operación no-lineal utilizada en procesamiento de imágenes usualmente

para reducir el ruido de “sal y pimienta”. El filtro de mediana es más efectivo que la convolución

cuando la meta es reducir el ruido preservando los contornos de la imagen. salida = medfilt2(imagen,[n

m]) aplica el filtro de mediana a la matriz de dos dimensiones imagen. Cada píxel de salida contiene

el valor de la mediana del vecindarionpormalrededor del píxel. En caso de que el píxel a evaluar se

encuentre en las orillas de la imagen, el vecindario se rellena con ceros. La imagen resultante salida es

del mismo tipo que la imagen de entrada imagen.

Ajuste de Contraste. imadjust(imagen)

La función imadjust mapea los valores de intensidad de imagen, a unos nuevos valores tales que

el 1% de la información se encuentre en las intensidades más bajas y altas de la imagen. Con esto se

incrementa el contraste de la imagen resultante. Si la función se llama sin más argumentos, los límites

de ajuste de contraste se encuentran con la función stretchlim, y es equivalente a llamar la función con

imadjust(imagen,stretchlim(imagen)).

Transformada Rápida de Fourier (2 Dimensiones). fft2(matriz)

La función ifft calcula la transformada discreta de Fourier en dos dimensiones de matriz, calculada

con el algoritmo de la trasformada rápida de Fourier (FFT). El resultado es del mismo tipo que matriz.

El tiempo de ejecución de fft2 depende del tamaño de matriz, y es más rápido cuando es potencia de

dos.

Transformada Inversa de Fourier (2 Dimensiones). ifft2(matriz)

La función ifft2 calcula la transformada discreta de Fourier inversa en dos dimensiones de matriz,

la cual fue calculada con el algoritmo de la trasformada rápida de Fourier (FFT). El resultado es del

mismo tipo que matriz. El tiempo de ejecución de ifft2 depende del tamaño de matriz, y es más rápido

cuando es potencia de dos.

Conjugado complejo. conj(a)

La función conj regresa el conjugado complejo dea, conj(a) = real(a) - i*imaginaria(a).


APÉNDICE E

Manual de usuario

Como ya se menciono antes, la interfaz gráfica de usuario fue creada con el fin de facilitar al

usuario la introducción de todos los datos necesarios y algunos opcionales. Esta interfaz esta realizada

en un tipo wizard, que es un asistente que guía al usuario paso a paso, en todos los elementos o datos

necesarios para llevar a cabo un proceso, en este caso, el de PIV.

Para ejecutar la aplicación con la interfaz de usuario, se necesita ejecutar MATALB, cambiar al

directorio donde se descomprimieron las funciones necesarias y a continuación ejecutar el comando:

>> pivgui

Si se requiere corregir alguna opción o dato de un paso anterior, se deberá presionar el botón Atrás

y regresar a la pantalla de los datos que se requieren modificar.

A continuación se describen todos los pasos que requiere el asistente.

E.1. Paso 1: Datos generales.

En este paso (véase Figura E.1) se requieren tres datos, que serán generales para todo el procesa-

miento:

Par de imágenes a procesar. Para introducir las imágenes, se proporciona un botón para

cada imagen con la etiqueta “Examinar”, el cuál al presionarlo, permitirá al usuario examinar

en cualquier ruta las imágenes que desee procesar.

El tiempo transcurrido entre las imágenes. Este dato es requerido para calcular la velocidad

de las regiones del fluido.

E.2. Paso 2: Datos de pre-procesamiento.

Para continuar al siguiente paso se deberá presionar el botón de “Siguiente”.

Figura E.1: Pantalla de datos generales.

E.2. Paso 2: Datos de pre-procesamiento.

En esta pantalla, como se puede observar en la Figura E.2, se pedirán las opciones que se tienen

disponibles en caso de que se necesite pre-procesar las imágenes.

Estiramiento de Contraste. Esta opción realiza la ecualización de ambas imágenes como lo

mencionado en la sección 3.1.

Eliminar Ruido. Como la opción lo dice, eliminará el ruido de las imágenes por medio del

algoritmo descrito en la sección 3.1. Como se menciona en dicha sección se requieren dos imágenes

más (de la misma sesión del experimento capturado), las cuáles se introducirán de la misma

manera como en el paso anterior.

Máscara de Procesamiento. Se permite introducir las coordenadas, que describen un rec-

tángulo, las cuales indicarán al programa que solo esa región se deberá procesar. Asimismo, se

permite introducir las coordenadas de forma interactiva. Al presionar el botón con la etiqueta

“Seleccionar”. Aparecerá la pantalla como se ve en la Figura E.3, la cuál permite al usuario


Apéndice E. Manual de usuario

seleccionar la región a procesar. Con la ayuda del ratón se debe de seleccionar un rectángulo.

Esto se hace al dejar presionado el botón izquierdo del ratón, arrastrar el puntero hasta donde

se desee y a continuación se libera el botón izquierdo. Esto se facilita aún más, ya que usuario

selecciona la región encima de la imagen a procesar.

Figura E.2: Pantalla de datos de pre-procesamiento.

E.3. Paso 3: Datos de procesamiento.

Aquí se presentan diferentes opciones a utilizar durante la etapa de procesamiento como lo descrito

en la sección 3.2.

Ventana. Indica el tamaño de la ventana o zonas de interrogación a tomar, en ancho (X) y alto

(Y ).

Iteraciones. Es el número de iteraciones que se procesarán.

Interpolación Subpíxel. Es el tipo de interpolación del máximo de correlación cruzada a nivel

subpíxel. Puede ser calculado de forma: Gaussiana, Parabólica o Centroide.


E.3. Paso 3: Datos de procesamiento.

Figura E.3: Pantalla de selección de máscara de procesamiento.


Apéndice E. Manual de usuario

Figura E.4: Pantalla de datos de procesamiento.

E.4. Paso 4: Datos de post-Procesamiento.

Finalmente, se introducen los datos que se tomarán en cuenta en la etapa de post-procesamiento.

Procesar resultados PIV. Al seleccionar esta opción se realizará la validación y sustitución

del mapa de vectores resultante como se describió en el capítulo 3.

Graficar Resultados. Esta opción permite al usuario visualizar el mapa de vectores, una gráfica

encima de la primera imagen a procesar y otra sin ésta.

Guardar Resultados. Al finalizar el procesamiento, si está opción está seleccionada, se guardará

la matriz de vectores en un archivo de textos. En el archivo las columnas de la matriz estará

separadas por tabulaciones, esto es útil si se requiere abrir el archivo en una hoja de cálculo como

lo es Excel de Microsoft.

Para llevar a cabo el procesamiento bastará presionar el botón Procesar y esperar por los resultados.


E.4. Paso 4: Datos de post-Procesamiento.

Figura E.5: Pantalla de datos de post-procesamiento.


APÉNDICE F

Pruebas Realizadas.

A continuación se presentan una serie de pruebas realizadas, los resultados solo son gráficos ya que

no se tienen los datos analíticos para poder comparar con un modelo teórico. Algunas de las imágenes

de prueba fueron obtenidas del reto PIV presentado en el congreso de PIV Challenge [32] anuales.

Apéndice F. Pruebas Realizadas.


UN NUEVO ALGORITMO EN LA TÉCNICA DE VELOCIMETRÍA POR ...

Documents

Transcript of UN NUEVO ALGORITMO EN LA TÉCNICA DE VELOCIMETRÍA POR ...