# HPC - Map gray levels on image

Last updated: 2022-12-16

Let us consider a grayscale bitmap with $$M$$ lines and $$N$$ columns, where the color of each pixel is encoded by an integer from 0 (black) to 255 (white). Given two values low, high, $$0 \leq \mathit{low} < \mathit{high} \leq 255$$, the function map_levels(img, low, high) modifies img so that pixels whose gray level is less than low become black, those whose gray level is greater than high become white, and those whose gray level is between low and high (inclusive) are linearly mapped to the range $$[0, 255]$$.

Specifically, if $$p$$ is the gray level of a pixel, the function computes the new gray level $$p'$$ as:

$p' = \begin{cases} 0 & \text{if}\ p < \mathit{low}\\ \displaystyle\frac{255 \times (p - \mathit{low})}{\mathit{high} - \mathit{low}} & \text{if}\ \mathit{low} \leq p \leq \mathit{high}\\ 255 & \text{is}\ p > \mathit{high} \end{cases}$

Figure 1 shows the image produced by the command

    ./simd-map-levels 100 255 < Yellow_palace_Winter.pgm > out.pgm

As an interesting practical application we provide the image C1648109 taken by the Voyager 1 probe on March 8, 1979. The image shows Io, one of the four Galilean moons of planet Jupiter. The Flight Engineer Linda Morabito was using this image to uncover the background stars to determine the precise location of the probe. To this aim, she needed to remap the gray levels so that the faint stars would be visible. This lead to one of the most important discoveries of modern planetary sciences: see by yourself by running the program simd-map-levels.c on the image C1648109 with low = 10 and high = 30, and look at what appears next to the disc of Io at ten o’clock…

The file simd-map-levels.c contains a serial implementation of the function map_levels() above. The goal of this exercise is to develop a SIMD version using GCC vector datatypes. We start by defining a vector datatype v4i that represents four integers:

typedef int v4i __attribute__((vector_size(16)));
#define VLEN (sizeof(v4i)/sizeof(int))

Then, the idea is to process the image four pixels at a time. However, the serial code

int *pixel = bmap + i*width + j;
if (*pixel < low)
*pixel = BLACK;
else if (*pixel > high)
*pixel = WHITE;
else
*pixel = (255 * (*pixel - low)) / (high - low);

is problematic because it contains conditional statements that can not be directly vectorized. To address this issue we use the selection and masking technique described in class.

Let pixels be a pointer to a v4i SIMD array. Then, the expression mask_black = (*pixels < low) produces a SIMD array of integers whose elements are -1 for those pixels whose gray level is less than low, 0 otherwise. mask_black can therefore be used as a bit mask to assign the correct values to these pixels.

Using the idea above, we can rewrite the code as follows:

v4i *pixels = (v4i*)(bmap + i*width + j);
const v4i mask_black = (*pixels < low);
const v4i mask_white = (*pixels > high);
const v4i mask_map = ??? ;
*pixels = ( (mask_black & BLACK) |
( ??? ) );

Note that the compiler automatically promotes BLACK and WHITE to SIMD vectors whose elements are all BLACK or WHITE, respectively. (The code above can be further simplified since (mask_black & BLACK) always produces a SIMD array whose elements are all zeros: why?).

The SIMD version requires that

1. The bitmap is aligned in memory at a memory address that is multiple of 16;

2. The image width is an integer multiple of the v4i SIMD vector width, i.e., an integer multiple of 4.

Both conditions above are satisfied by the sequential code and the input images.

To compile:

    gcc -std=c99 -Wall -Wpedantic -O2 -march=native simd-map-levels.c -o simd-map-levels

To execute:

    ./simd-map-levels low high < input_file > output_file

dove $$0 \leq \mathit{low} < \mathit{high} \leq 255$$.

Example:

    ./simd-map-levels 10 30 < C1648109.pgm > C1648109-map.pgm