Skip to content
GitLab
Commit 18fe541e authored by Pablo Quesada Barriuso's avatar Pablo Quesada Barriuso
Browse files

Update Readme and minor cleanup

parent ed96d9f2
Hide whitespace changes
Inline Side-by-side
DWT_1D.c
View file @ 18fe541e
......@@ -24,33 +24,12 @@
* Fecha: 03-10-2013
* Notas: Esta funcion usa BaseType (double) como tipo de datos por defecto. Ver README
*
* OpenMP: Si. Cada thread calcula la 1D_DWT para un bloque consecutivo de pixel vectors
*
* @warning: Esta funcion modifica los datos de entrada!
*
* Solo hacemos convolucion con filtro bajo ya que sólo nos interesa la aproximación
* La salida tiene los resultados de aproximacion en la primera mitad del vector (decimated)
* La convolución se hace de forma circular respecto a los últimos elementos del vector
*
* @param h_idata Imagen hiperespectral con un pixel vector por fila
* @param LEVELS Número de veces que se aplica la 1D-DWT
* @param NCOLS Número de bandas espectrales (dimensión del pixel vector)
* @param ROWS Número total de pixels vector de la imagen (width x height)
*
* Return:
*
* @param h_odata 1D-DWT por filas
*/
* OpenMP: Each thread do a convolution in a one dimension for each pixel-vetor in a chunk
* */
void computeDWT_1D_OMP(baseType* h_idata, baseType* h_odata, int NCOLS, int NROWS)
{
/* OMP SETUP */
// enable nested parallelism.
omp_set_nested(1);
// disables dynamic adjustment (by the run time system) of the number
// of threads available for execution of parallel regions.
omp_set_dynamic(0);
int chunk = NROWS / omp_get_max_threads();
printf("(CONFIG) OMP_NUM_THREADS : %i\n", omp_get_max_threads());
......@@ -85,16 +64,9 @@ void computeDWT_1D_OMP(baseType* h_idata, baseType* h_odata, int NCOLS, int NROW
// do the convolution
for (int j = 0; j < M1; ++j)
{
// NOTA: Hacemos un padding circu, lar en la misma fila! Por eso
// tenemos que sumar NCOLS*k para volver a tomar el primer
// elemento de la fila que estamos procesando
// NOTE: Circular padding for each pixel vector
lp += lpf[ j ] * i_ptr[ ( (idx_r + j) % S) + (NCOLS*k) ];
}
// El resultado de la convolución en la iteración actual se almacena
// en la posición final de salida (respecto a cada pixel vector) calculando la decimation
// Esa posición depende del tamaño actual del vector, por lo que hacemos un
// modulo S, y también depende de la posición del valor actual que es 'i'.
// Ver variable idx_w
o_ptr[ idx_w % S + (NCOLS*k) ] = lp;
}
} // implicit barrier
......@@ -106,9 +78,7 @@ void computeDWT_1D_OMP(baseType* h_idata, baseType* h_odata, int NCOLS, int NROW
S = iDivUp(S,2);
}
// ey man!! copy the final result to h_odata!
// Si los punteros son iguales, se ha hecho un numero par de iteraciones
// y por lo tanto el resultado final esta en h_idata
// keep output in the h_odata buffer in case of even iterations
chunk = (NROWS*NCOLS) / omp_get_max_threads();
if ( o_ptr == h_odata )
#pragma omp parallel for schedule(static, chunk) default(shared)
......
EMP.c
View file @ 18fe541e
......@@ -149,8 +149,7 @@ void MM_imerode(Npp8u *p_in, MM_SE se, Npp8u *p_out, int w, int h)
// erosion within the kernel
if ( se.box[ kdx++ ] == 1 ) p = min( p, v );
}
}
}
// save erosion in the output buffer
p_out[ idx ] = p;
}
......@@ -208,22 +207,15 @@ void MM_imcomplement(Npp8u *i_ptr, Npp8u *o_ptr, int w, int h)
}
// morphological grayscale reconstruction
// References : [1] Vincent, Luc., "Morphological grayscale reconstruction in
// References : Vincent, Luc., "Morphological grayscale reconstruction in
// image analysis: Applications and efficient algorithms." Image
// Processing, IEEE Transactions on 2.2 (1993): 176-201.
// Section 4.2 Sequential reconstruction algorithm
// I: mask image
// J: marker image
// Reconstruction is determined directly in J
void MM_imreconstruct4_2(Npp8u *J, Npp8u *I, int w, int h)
{
// stability flag
int stability;
#ifdef DEBUG
int iterations = 0;
#endif
// repeat until stability
do {
// reset flag
......@@ -301,45 +293,15 @@ void MM_imreconstruct4_2(Npp8u *J, Npp8u *I, int w, int h)
}
/*
* Aplica un Morphological Profile (MP) a cada banda de la imagen de entrada, usando
* 4 openings y 4 closings por reconstrucción con un elemento structural (SE) de tamaño
* creciente (ver morphtools.h).
*
* Al aplicar el MP a cada banda, se crea un Extended Morphologial Profile EMP
*
* @param h_idata Imagen hiperespectral con una banda de w x h pixels detrás de otra
* @param width Ancho de la imagen
* @param height Alto de la imagen
* @param bands Número de bandas espectrales
* @param mpSize Tamano del perfil morfologico para una banda, es decir, numero
* de opening por reconstruccion + numero de closing por reconstruccion + 1
*
* Actualmente mpSize solo se usa para reservar memoria, ya que en esta funcion
* se ha preparado el blucle para hacer siempre cuatro opening y cuatro closing
* por reconstruccion. Por tanto, mpSize = 9 SIEMPRE
*
* Return:
* Do a Morphological Profile for each spectral band. Thus, the result is an Extended Morphological Profile
*
* @param h_odata Extended Morphological Profile
*
* Notas: Antes de llamar a esta función hay que normalizar la imagen entre [0..255]
*
* Esta función usa la misma implementación para opening y closing by reconstruction, así que en el caso
* de "closing by reconstruction" primero hay que complementar la entrada y el resultado de la erosion, y
* al final complementar el resultado de la reconstrucción. Otra opción sería implementar funciones separadas
* para opening y closing by reconstruction (no es el caso en este esquema).
* OpenMP: Each thread get a set of spectral bands and do all the opening and closing by reconstruction
* with an structuring element of increasing size
*/
void computeDWT_EMP_CPU_4_2(Npp8u* h_idata, Npp8u* h_odata, int width, int height, int bands, int mpSize)
{
/* OMP SETUP */
// enable nested parallelism.
omp_set_nested(1);
// disables dynamic adjustment (by the run time system) of the number
// of threads available for execution of parallel regions.
omp_set_dynamic(0);
// setup the chunk size based on the number of spectral bands and max number of threads
int chunk = iDivUp(bands, omp_get_max_threads());
printf("(CONFIG) OMP_NUM_THREADS : %i\n", omp_get_max_threads());
......@@ -375,13 +337,11 @@ void computeDWT_EMP_CPU_4_2(Npp8u* h_idata, Npp8u* h_odata, int width, int heigh
// filter box size for a SE of radius 'r'
SE.r = rSE[ r ];
// morphological grayscale erosion
MM_imdilate(i_ptr, SE, o_ptr, width, height);
// complement the erosion, that is, the marker
MM_imcomplement(o_ptr, o_ptr, width, height);
// morphological grayscale reconstruction
MM_imreconstruct4_2(o_ptr, c_ptr, width, height);
......
README.md
View file @ 18fe541e
......@@ -9,7 +9,9 @@ Use Case: Stencil computation applied to remote sensing hyperspectral preprocess
In parallel programming we frequently find computing patterns that are common in parallel algorithms. An example is the stencil computation where each output element is computed with data from its neighborhood defined by a mask. For this computation the number of operations per pixel is very large. We will see in this use case how to speed up two common pre-processing steps in hyperspectral analysis using OpenMP, from spectral feature reduction using wavelets to morphological profiles by means of the same parallel pattern.”
## Getting started
A good point to start is the youtube recorded session available online at []. The bibliography has some papers related directly to the use case. It is a good reading to enjoy an evening taking a coffee.
A good point to start the HDCRS Summer School - Day 2 - May 31, 2022 recorded session available online at https://www.youtube.com/watch?v=55E9DyoMs1k. Move forward until time 1:01:23.
The bibliography has some papers related directly to the use case. It is a good reading to enjoy an evening taking a coffee.
## Prerequisites
......@@ -25,11 +27,11 @@ g++ WT-EMP-Scheme.c -O2 -fopenmp -lm -o WT-EMP-Scheme
./WT-EMP-Scheme data/PaviaUniversityPV.raw 5
```
You need a hyperspectral image in raw format and stored as pixel vectors. You can find the well-know image of Pavia University in that format in the folder named 'data'.
You need a hyperspectral image in raw format and stored as pixel vectors. You can find the well-know image of Pavia University in that format in the `data/` folder.
You must specify the number of times you want to reduce the spectral dimensionality to a half.
You must specify the number of times you want to reduce the spectral dimensionality to a half as the second parameter.
This version create a fixed Morphological Profile with 8 opening and 8 closing for each band.
This version create a fixed Morphological Profile with 8 opening and 8 closing for each spectral band.
## Output
......@@ -54,4 +56,24 @@ The following is just an example:
(TIME) NORMALIZE EMP 0.095508 sg
[info] Saved 'h_emp_opening_u8.pgm', 340 x 610 (207400 pixels) (202 Kb)
[info] Saved 'h_emp_closing_u8.pgm', 340 x 610 (207400 pixels) (202 Kb)
```
\ No newline at end of file
```
## Bibliography
[1] Quesada-Barriuso, P. Argüello, F., Heras D. B., Benediktsson, J. A., "Wavelet-Based Classification of Hyperspectral Images Using Extended Morphological Profiles on Graphics Processing Units," in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 8, no. 6, pp. 2962-2970, June 2015, doi: 10.1109/JSTARS.2015.2394778.
[2] T.G. Mattson, B.A. Sanders and B.L. Massingill, Patterns for Parallel Programming, Addison-Wesley, 2005.
[3] Kirk, D. B., & Wen-mei, W. H. (2013). Programming massively parallel processors: a hands-on approach. Morgan Kaufmann Publishers.
[4] Benediktsson, J. A., Palmason, J. A., & Sveinsson, J. R. (2005). Classification of hyperspectral data from urban areas based on extended morphological profiles. IEEE
Transactions on Geoscience and Remote Sensing, 43(3), 480-491.
[5] Vincent, Luc., "Morphological grayscale reconstruction in image analysis: Applications and efficient algorithms." Image Processing, IEEE Transactions on 2.2 (1993): 176-201.
## Pavia University Hyperspectral Image
The Pavia University hyperspectral image was provided by Prof. Paolo Gamba from the Telecommunications and Remote Sensing Laboratory, Pavia university (Italy).
This is image is hosted in the Hyperspectral Remote Sensing Scenes collected by: M Graña, MA Veganzons, B Ayerdi.
The version included in this repository is in stored as pixel vector using 4 bytes per value in a raw format.
\ No newline at end of file
WT-EMP-Scheme.c
View file @ 18fe541e
......@@ -48,7 +48,7 @@ int main(int argc, char *argv[])
// host output data
baseType *h_odata = (baseType*)malloc(data.size*sizeof(baseType));
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////
// 1D-DWT Processing :: NOTE: This functions updated the input date
computeDWT_1D_OMP(data.x, h_odata, data.bands, data.rows*data.cols);
......@@ -75,7 +75,8 @@ int main(int argc, char *argv[])
// Export first component extracted by wavelet as pgm for a visual check
exportPGM("h_dwt_u8.pgm", data.xu8, data.cols, data.rows);
// EXTENDED MORPHOLOGOCAL PRPOFILE //
///////////////////////////////////////////////////////////////////
// EXTENDED MORPHOLOGOCAL PRPOFILE
// We will do 4 opening and 4 closing by reconstruction by default
int opn = 4;
......
globals.h
View file @ 18fe541e
......@@ -53,23 +53,6 @@ int LEVELS = 1;
// 1D-DWT filters for Feature Reduction :: low pass filter for CPU (the kernel for the convolution)
const baseType lpf[ M1 ] = {0.0378284555, -0.0238494650, -0.1106244044, 0.3774028556, 0.8526986790, 0.3774028556, -0.1106244044,-0.0238494650, 0.0378284555};
// 1D-DWT REMAINING COEFFICIENTS
int NCOEFFS = 0;
/**
* Adding message to assert
* http://stackoverflow.com/questions/3767869/adding-message-to-assert
*/
#define ASSERT(condition, message) \
do { \
if (! (condition)) { \
std::cerr << "Assertion `" #condition "` failed in " << __FILE__ \
<< " line " << __LINE__ << ": " << message << std::endl; \
std::exit(EXIT_FAILURE); \
} \
} while (false)
// define some color
#define KNRM "\x1B[0m"
......@@ -95,7 +78,6 @@ typedef double timenfo;
resnfo start, end, start_EAP, end_EAP;
timenfo elapsedTime, formatTime = 0.0, totalTime = 0.0;
timenfo tiempo, t_transfer = 0.0, t_malloc = 0.0, t_free = 0.0, t_compute = 0.0, t_cpu = 0.0;
void myElapsedtime(const resnfo start, const resnfo end, timenfo *const t)
{
......@@ -164,16 +146,9 @@ struct HSI {
};
/**
* @brief Load hyperspectral dataset from a RAW file.
*
* Load hyperspectral dataset from a RAW file.
*
* This function loads a hyperspectral dataset (stored as INT) into a DOUBLE variable in pixel-vector format.
*
* @note Modified for using HSI structures
* @see struct HSI {...} in imageHeader.h
*
* @param[in] filename The name of the file for reading the data
*
* @param[out] data A structure for managing the hyperspectral data
*/
void read_dataset_raw(HSI &data, const char* filename)
{
......@@ -245,14 +220,13 @@ void read_dataset_raw(HSI &data, const char* filename)
}
/*
* change HSI 1 pixel vector per row to CUBE data
* Change HSI format from pixel vector to consecutive bands
*/
baseType* Rows2HSI(baseType* iPtr, int w, int h, int b, int bb)
{
baseType* oPtr = NULL;
oPtr = (baseType *) malloc( sizeof(baseType) * w * h * b);
/* change HSI 1 pixel vector per row to CUBE data */
for (int row = 0; row < h*w; ++row)
for (int band = 0; band < b; ++band)
oPtr[ (band*h*w) + row ] = iPtr[ row*bb + band ];
......@@ -265,37 +239,10 @@ baseType* Rows2HSI(baseType* iPtr, int w, int h, int b, int bb)
/*
* Auxiliary function: normalize
*
* Normalizar valores entre [a..b] (0..255 para imagenes en escala de grises)
*
* @param i_ptr Imagen hiperespectral con una banda de w x h pixels detrás de otra
* @param a Valor mínimo
* @param b Valor máximo
* @param w Ancho de la imagen
* @param h Alto de la imagen
* @param d Número de bandas espectrales (por defecto 1)
*
* Return: baseType* o_ptr Una nueva imagen de tipo UNSIGNED CHAR con valores en el rango [a..b]
* Normalize data between a..b] (0..255 for grayscale images)
*/
baseType* normalize32f(unsigned char* i_ptr, baseType a, baseType b, int w, int h, int d)
{
/*
Suppose you have a range or scale from A to B and you want to convert
it to a scale of 1 to 10, where A maps to 1 and B maps to 10.
Furthermore, we want to do this with a linear function, so that for
example the point midway between A and B maps to halfway between 1 and
10, or 5.5.
Then the following (linear) equation can be applied to any number
x on the A-B scale:
y = 1 + (x-A)*(10-1)/(B-A)
Note that if x = A, this gives y = 1+0 = 1 as required, and if x = B,
y = 1+(B-A)*(10-1)/(B-A) = 1+10-1=10, as required.
You can use this equation even if A > B.
*/
baseType* o_ptr = ( baseType*)calloc(w*h*d, sizeof(baseType) );
baseType A = INT_MAX, B = 0.0;
......@@ -313,6 +260,10 @@ You can use this equation even if A > B.
return o_ptr;
}
/*
* Export the array v into a file in a raw format
*/
int exportRAW(baseType* v, const char* o_file, int w, int h, int b)
{
FILE* fp;
......@@ -345,7 +296,7 @@ int exportRAW(baseType* v, const char* o_file, int w, int h, int b)
return 0;
}
/*
* export PGM data to file
* Export the arra v into a file a PGM format
*/
int exportPGM(const char *filename, unsigned char* v, int width, int height)
{
......@@ -376,8 +327,6 @@ int iDivUp(int a, int b) { return (a % b != 0) ? (a / b + 1) : (a / b); }
/**
* Configure OpenMP and display basic information
*
* @param[out] The number of threads used in the execution
*/
int openMPInfo()
{
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment