Acerca de

Bueno, este es un pequeño proyecto que hice en la universidad, acerca de computación multinúcleo y en clúster. Básicamente nos pidieron que escogiéramos una transformada de wavelet, y que la paralelizásemos lo máximo posible. Como éramos libres de escoger cualquier algoritmo, nos decantamos por usar Piecewise Linear Haar, una transformada discreta orientada a imágenes.

Técnicamente era una clase sobre cómo usar OpenMP y Message Passing Interface, pero decidí por mi parte ir un poco más allá y juguetear por mera diversión con las extensiones SIMD disponibles en los procesadores Intel.

A continuación se encuentra la memoria final del proyecto:

Paralelización de PLHaar

Aplicación de una transformada discreta de wavelet en imágenes

Objetivo y utilidad del algoritmo

El objetivo de nuestro proyecto es desarrollar una versión paralelizada del algoritmo PLHaar en su versión discreta. El algoritmo PLHaar es una DWT (Discrete Wavelet Transform) y se emplea en tareas de compresión además de otros campos.

Las transformadas Wavelet descomponen la señal una serie de subseñales organizadas de forma jerárquica por relevancia de los datos. Cada transformación separa la señal en dos subseñales, una aproximación de la señal original y una que recoge los detalles extraídos de la aproximación. Estas subseñales pueden ser utilizadas para reconstruir la señal original sin ninguna pérdida.

En este proyecto vamos a utilizar esta transformada para el tratamiento de imágenes. Como reza el enunciado del trabajo, se utilizará sobre una imagen de un mínimo de 32x2 píxeles, en escala de grises o en RGB, con 4 coeficientes y un nivel arbitrario de profundidad.

Esto nos permitirá obtener una serie de imágenes que se comportarán como las señales propuestas en el párrafo anterior, de forma que tendremos una imagen reducida de la original y tres más adicionales donde estarán representados los detalles verticales, horizontales y diagonales de la primera imagen procesada, obteniendo así los 4 coeficientes.

En cada iteración se reducirá el tamaño de la imagen original, tantas veces como niveles de profundidad se quieran obtener.

Finalizado el proceso, el algoritmo es capaz de recomponer la imagen primigenia sin ninguna pérdida.

Ejemplos

Imagen original

Imagen original

Imagen tras una iteración

Imagen tras una iteración

Imagen tras cuatro iteraciones

Imagen tras cuatro iteraciones

Dibujo original (©TysonTan; CC-BY-SA)

Dibujo original (©TysonTan; CC-BY-SA)

Dibujo tras tres iteraciones. El contraste ha sido aumentado para mejorar la visibilidad.

Dibujo tras tres iteraciones. El contraste ha sido aumentado para mejorar la visibilidad.

Descripción del algoritmo

PLHaar dispone de dos entradas: A y B, que son bytes que contienen dos valores de un mismo canal de dos píxeles contiguos. La salida es L (lowpass), que representa en la nueva imagen una aproximación de la original; y H (highpass), que representa la diferencia entre los píxeles.

La función PLHaar tiene hay dos entradas y dos salidas, y se encuentra definida a trozos:

  1. Si ab0 y a>b:
    • h=(a-b)mod256
    • l=a
  2. Si ab0 y ab:
    • h=(a-b)mod256
    • l=b
  3. Si ab<0 y a>b:
    • h=(-b)mod256
    • l=(a+b-128)mod256
  4. Si ab<0 y ab:
    • h=a
    • l=(a+b-128)mod256

Este proceso se realiza en bloques de 2x2 píxeles para obtener cuatro señales: uno es la aproximación real y las otras 3 contienen información sobre las diferencias horizontales, verticales y diagonales de los píxeles que forman la original.

A continuación una muestra de la implementación del código, que es una versión ligeramente retocada del código de ejemplo original:

void plhaar_int(uint8_t a, uint8_t b, uint8_t * l, uint8_t * h) {
    const uint8_t c = 128;
    const uint8_t s = (a < c), t = (b < c);

    a += s; b += t;         // asymmetry: nudge origin to (+0,+0)
    if (s == t) {           // A * B > 0?
        a -= b - c;         // H = A - B
        if ((a < c) == s) { // |A| > |B|?
            b += a - c;     // L = A (replaces L = B)
        }
    } else {                // A * B < 0
        b += a - c;         // L = A + B
        if ((b < c) == t) { // |B| > |A|?
            a -= b - c;     // H = -B (replaces H = A)
        }
    }
    a -= s; b -= t;         // asymmetry: restore origin
    *l = b; *h = a;         // store result
}
Diagrama explicando cómo se emplea la función

Diagrama explicando cómo se emplea la función

Propuestas de paralelización

En nuestro proyecto vamos a proponer la utilización tanto de paralelización de datos como de las tareas que realiza el algoritmo.

Para la paralelización de los datos vamos a utilizar las extensiones de SIMD de Intel® para que cada núcleo trabaje con múltiples píxeles a la vez, aumentando la eficiencia de las operaciones a realizar.

Para la paralelización de las tareas utilizaremos la librería OpenMP. Dividiremos la imagen a procesar en secciones de píxeles consecutivos de tamaño inversamente proporcional al número de núcleos disponibles para el cómputo. De esta manera cada núcleo se encargará de una zona distinta de la imagen.

También se probará a realizar una implementación empleando MPI, de modo que haya varios procesadores operando con una misma imagen.

Paralelización con SSE

Empleando SSE v4.1, es posible emplear registros especiales de 128 bits, que permiten operar en paralelo. Por lo tanto, empleando cuatro registros de 128 bits, que pueden almacenar 16 bytes, es posible operar con 64 píxeles en una sola ejecución del algoritmo.

Hay un problema con SSE, y es que debido al paralelismo no es posible el uso de operaciones condicionales clásicas, los if no existen por así decirlo.

Sin embargo, existe una operación determinada blend desde SSE v4.1, que dados registros SSE A, B y X, de 16 bytes cada uno, ejecuta el equivalente al siguiente código en C:

for (int i = 0; i < sizeof(x); i++) {
    if (x[i] == 0)
        salida[i] = a[i];
    else
        salida[i] = b[i];
}

Como se puede observar en el código original adjunto en el punto "Descripción del algoritmo", hay cuatro posibles funciones que emplear en base a las entradas.

Lowpass tree

Lowpass tree

Highpass tree

Highpass tree

Si se define la salida como un árbol dependiendo de las entradas, para adaptar el código secuencial a SSE, se ha optado por calcular todos los valores posibles para las salidas, y después podar el árbol mediante blends.

#include <smmintrin.h>

void plhaar_sse(
        __m128i valuesA,
        __m128i valuesB,
        __m128i lowpass,
        __m128i highpass) {
    __m128i s = _mm_cmplt_epi8(valuesA, _mm_setzero_si128());
    __m128i t = _mm_cmplt_epi8(valuesB, _mm_setzero_si128());

    valuesA = _mm_sub_epi8(valuesA, s);
    valuesB = _mm_sub_epi8(valuesB, t);

    __m128i ifA = _mm_sub_epi8(valuesA, valuesB);
    __m128i aLessZero = _mm_xor_si128(s,
            _mm_cmplt_epi8(ifA, _mm_setzero_si128()));
    __m128i ifB = _mm_blendv_epi8(valuesA, valuesB, aLessZero);

    __m128i ifNotB = _mm_add_epi8(valuesB, valuesA);
    __m128i minusB = _mm_sub_epi8(_mm_setzero_si128(), valuesB);
    __m128i bLessZero = _mm_xor_si128(t,
            _mm_cmplt_epi8(ifNotB, _mm_setzero_si128()));
    __m128i ifNotA = _mm_blendv_epi8(minusB, valuesA, bLessZero);

    __m128i positive = _mm_xor_si128(s, t);
    highpass = _mm_add_epi8(_mm_blendv_epi8(ifA, ifNotA, positive), s);
    lowpass = _mm_add_epi8(_mm_blendv_epi8(ifB, ifNotB, positive), t);
}

Para realizar la implementación con SSE, en lugar de operar directamente con ensamblador, se ha optado por emplear los llamados intrinsics: macros que hacen de puente entre C y ensamblador, pudiendo especificar exactamente las instrucciones de ensamblador a ejecutar, pero permitiendo a GCC reordenar código y adaptarlo al número de registros SSE de que el sistema dispone (ocho cuando el procesador está en modo protegido -32 bits-, y dieciséis en modo largo -64 bits-).

En estos intrinsics, se puede notar un patrón. Tomando por ejemplo _mm_add_epi8:

Por lo tanto, _mm_add_epi8 ejecuta una suma de dos registros que contienen enteros de 8 bits, y devuelve el resultado.

Cada pseudo-variable __m128i, que representa un registro de 128 bits que GCC escogerá, puede almacenar 16 bytes. Por tanto, cada hilo ejecutándose en un procesador compatible es capaz de ejecutar 16 operaciones con enteros de 8 bits en paralelo, lo que resulta en un speedup enorme, y por ende una mayor eficiencia, manteniendo el mismo número de procesadores.

OpenMP implementation

Refinando la propuesta anterior, vamos a discutir cómo hemos optado por paralelizar nuestro código.

Tal y como se comenta en el apartado anterior, la idea era dividir la imagen en pequeñas imágenes independientes para su procesado independiente. Esta habría sido una buena aproximación, pero se puede hacer de una forma más simple.

Nuestro código son en realidad dos bucles for anidados que computan la imagen primero en horizontal y luego en vertical. Por la escasa dependencia de datos que comentamos en la sección correspondiente del documento, es posible paralelizar ambos cálculos de forma que resulta innecesario establecer un orden concreto. Por ello hemos optado por la siguiente aproximación:

void plhaar2d_transform_naive_mp(
    const uint8_t * values,
    size_t width,
    size_t height,
    uint8_t * aver,
    uint8_t * hori,
    uint8_t * vert,
    uint8_t * diag
) {
    #pragma omp parallel for \
    collapse(2) \
    schedule(static) \
    shared(values, width, height, aver, hori, vert, diag)
    for (size_t y = 0; y < height; y += 2) {
        for (size_t x = 0; x < width; x += 2) {
            /*
             * Load four pixels into local thread-private variables
             */
            uint8_t topLeft     = values[(x + 0) + (y + 0) * width];
            uint8_t topRight    = values[(x + 1) + (y + 0) * width];
            uint8_t bottomLeft  = values[(x + 0) + (y + 1) * width];
            uint8_t bottomRight = values[(x + 1) + (y + 1) * width];

            /*
             * Compute PLHaar horizontally
             * Top left with top right, and bottom left with bottom right
             */
            uint8_t topLow, topHigh, bottomLow, bottomHigh;
            plhaar_int(topLeft,    topRight,    &topLow,    &topHigh);
            plhaar_int(bottomLeft, bottomRight, &bottomLow, &bottomHigh);

            /*
             * Now compute PLHaar vertically
             * Lowpass from top row with lowpass from bottom row, and
             * highpass from top row with highpass from bottom row.
             * The output is written to aver (average), hori (horizontal),
             * vert (vertical) and diag (diagonal) vectors.
             */
            size_t offset = y / 2 * width / 2 + x / 2;
            plhaar_int(topLow,  bottomLow,  aver + offset, hori + offset);
            plhaar_int(topHigh, bottomHigh, vert + offset, diag + offset);
        }
    }
}

Como podemos ver, el pragma empleado finalmente es:

#pragma omp parallel for collapse(2) schedule(static)

Hemos omitido las directivas private() ya que el compilador automáticamente asigna como private los índices de los bucles for.

Para nuestra función hemos considerado que el valor óptimo para el planificador sería *static, debido a que siempre será una carga fija. Es decir, el coste de cada iteración del bucle es siempre el mismo, con lo que no tiene sentido usar asignación dinámica.

Por último empleamos la cláusula collapse(2) para que el planificador distribuya la carga entre el total de iteraciones a computar agregando las de ambos bucles.

El código fuente está disponible en la rama "master" de GitHub.

Benchmark de SSE+MP

A continuación se muestran gráficas comparando el rendimiento de diferentes configuraciones: variando el número de núcleos, y variando la implementación.

Para dibujar las siguientes gráficas, se ha ejecutado cada configuración cuatro veces, ejecutadas secuencialmente sin pausas, y se ha cogido el mejor tiempo de cada configuración.

Para cada ejecución, se ha generado una imagen nueva, generando así gran cantidad de I/O, evitando que el sistema cachease trozos del programa y que pudiera beneficiar a ejecuciones posteriores.

Variando el número de hilos

Ejecuta secuencial, OpenMP en 2, 4, 6, 8, 12, 16, 24, 32 y 128 hilos, con imágenes entre 512x512 y 32768x32768 píxeles.

Esta prueba se ha ejecutado en un servidor con dos Sandy Bridge Xeon E5-2620 processors a 2.1GHz, con 6 núcleos físicos y 12 virtuales cada uno, para tener un total de 24 núcleos. El sistema dispone de 64GB de memoria DDR3 a 1333MHz.

Para los siguientes gráficos hemos decidido:

Tiempo absoluto de ejecución

Tiempo absoluto de ejecución

Aceleración con OpenMP

Aceleración con OpenMP

Eficiencia de OpenMP

Eficiencia de OpenMP

Como se puede observar en los gráficos:

Variando la implementación

Ejecución secuencial, OpenMP en 24 núcleos, SSE, y OpenMP más SSE en 24 núcleos; con imágenes entre 512x512 y 32768x32768 píxeles.