Appendix A

// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

// Kernel Code that is executed on GPU to convert Raw to SC/SFI

// Two kernels (rows and columns) execute sequentially

// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

// 24-bit multiplication is faster on G80,

// but we must be sure to multiply integers

// only within [-8M, 8M - 1] range

#define IMUL(a, b) __mul24(a, b)

#ifndef _SPECKLECONTRAST_H_

#define _SPECKLECONTRAST_H_

////////////////////////////////////////////////////////////////////////////////////////

// Kernel configuration (constants for GPU calculations)

////////////////////////////////////////////////////////////////////////////////////////

// Assuming ROW_TILE_W, KERNEL_RADIUS_ALIGNED and dataW

// are multiples of maximum coalescable read/write size,

// all global memory operations are coalesced in convolutionRowGPU5()

#defineROW_TILE_W 128

#define WINDOW_RADIUS_ALIGNED 16

// Assuming COLUMN_TILE_W and dataW are multiples

// of maximum coalescable read/write size, all global memory operations

// are coalesced in convolutionColumnGPU5()

#define COLUMN_TILE_W 16

#define COLUMN_TILE_H 48

////////////////////////////////////////////////////////////////////////////////////////

// Loop unrolling templates, needed for best performance (5×5 window)

////////////////////////////////////////////////////////////////////////////////////////

template<int i> __device__ int convolutionRowMeans5(int *data){

return data[2 - i] + convolutionRowMeans5<i - 1>(data);

}

template<> __device__ int convolutionRowMeans5<-1>(int *data){

return 0;

}

template<int i> __device__ int convolutionRowSquares5(int *data){

return data[2 - i] * data[2 - i] + convolutionRowSquares5<i - 1>(data);

}

template<> __device__ int convolutionRowSquares5<-1>(int *data){

return 0;

}

template<int i> __device__ int convolutionColumnFull5(int *data){

return data[(2 - i) * COLUMN_TILE_W] + convolutionColumnFull5<i - 1>(data);

}

template<> __device__ int convolutionColumnFull5<-1>(int *data){

return 0;

}

////////////////////////////////////////////////////////////////////////////////////////

// Row convolution filter (5×5)

////////////////////////////////////////////////////////////////////////////////////////

// dataW and dataH are the width and height of raw data, respectively

// d_Data, d_Means, d_Squares are preallocated device memory for the raw input,

//strip sums of mean buffer, and strip sums of squares buffer, respectively

// WINDOW_RADIUS is NOT current implemented, 5×5 window filter is hardcoded

__global__ void convolutionRowGPU5(

int *d_Data,

int *d_Means,

 int *d_Squares,

int dataW,

int dataH,

 int WINDOW_RADIUS

){

// Shared memory allocation

__shared__ int data[2 + ROW_TILE_W + 2];

// Current tile and apron limits, relative to row start

const int  tileStart = IMUL(blockIdx.x, ROW_TILE_W);

const int  tileEnd = tileStart + ROW_TILE_W - 1;

const int  apronStart = tileStart - 2;

const int  apronEnd = tileEnd + 2;

// Clamp tile and apron limits by image borders

const int tileEndClamped = min(tileEnd, dataW - 1);

const int apronStartClamped = max(apronStart, 0);

const int apronEndClamped = min(apronEnd, dataW - 1);

// Row start index in d_Data[]

const int rowStart = IMUL(blockIdx.y, dataW);

// Aligned apron start. Assuming dataW and ROW_TILE_W are multiples

// of half-warp size, rowStart + apronStartAligned is also a

// multiple of half-warp size, thus having proper alignment

// for coalesced d_Data[] read.

const int apronStartAligned = tileStart - WINDOW_RADIUS_ALIGNED;

// Current load postion (of this thread)

const int loadPos = apronStartAligned + threadIdx.x;

// Set the entire data cache contents

// Load global memory values, if indices are within the image borders,

// or initialize with zeroes otherwise

if(loadPos > = apronStart){

 //shared memory postion

 const int smemPos = loadPos - apronStart;

 data[smemPos] =

  ((loadPos > = apronStartClamped) && (loadPos < = apronEndClamped)) ?

  d_Data[rowStart + loadPos] : 0;

 }

// Ensure the completeness of the loading stage

// because results, emitted by each thread depend on the data,

// loaded by another threads

__syncthreads();

// Current write position

const int writePos = tileStart + threadIdx.x;

// Assuming dataW and ROW_TILE_W are multiples of half-warp size,

// rowStart + tileStart is also a multiple of half-warp size,

// thus having proper alignment for coalesced d_Result[] write.

if(writePos < = tileEndClamped){

 const int smemPos = writePos - apronStart;

 int sumMeans = 0;

 int sumSquares = 0;

// Loop unrolling section - computes rolling sums and square sums

// only uses unrolling templates (defined at begining) if defined

// if not defined, uses for loop

#ifdef UNROLL_INNER

 sumMeans = convolutionRowMeans5<2 * 2>(data + smemPos);

 sumSquares = convolutionRowSquares5<2 * 2>(data + smemPos);

#else

 for(int k = -2; k < = 2; k++){

  sumMeans + = data[smemPos + k];

  sumSquares + = data[smemPos + k] * data[smemPos + k];

  }

#endif

   // Writes the strip sums to preallocated global memory

  d_Means[rowStart + writePos] = sumMeans;

   d_Squares[rowStart + writePos] = sumSquares;

 }

}

////////////////////////////////////////////////////////////////////////////////////////

// Column convolution filter + Speckle Contrast + Speckle Flow Index (5×5)

////////////////////////////////////////////////////////////////////////////////////////

// dataW and dataH are the width and height of raw data, respectively

// d_MeansRow, d_SquaresRow are preallocated device memory for strip sums of

// mean and squares values obtained in the previous kernel, respectively

// d_MeansFull, d_SquaresFull are preallocated device memory for full sums of

// mean and squares values buffer, respectively

// d_SC, d_SFI are preallocated device memory for speckle contrast and speckle

// flow index values buffer, respectively

// expT is the exposure time of the raw images

// smemStride, gmemStride are the incremental position movement of the

// shared and globably memory processing positions, respectively

// WINDOW_RADIUS is NOT current implemented, 5×5 window filter is hardcoded

__global__ void convolutionColumnGPU5(

int *d_MeansRow,

int *d_SquaresRow,

int *d_MeansFull,

int *d_SquaresFull,

float *d_SC,

float *d_SFI,

int dataW,

int dataH,

int WINDOW_RADIUS,

float expT,

int smemStride,

int gmemStride

){

// Shared memory allocation

__shared__ int rowMeans[COLUMN_TILE_W * (2 + COLUMN_TILE_H + 2)];

__shared__ int rowSquares[COLUMN_TILE_W * (2 + COLUMN_TILE_H + 2)];

// Current tile and apron limits, in rows

const int tileStart = IMUL(blockIdx.y, COLUMN_TILE_H);

const int tileEnd = tileStart + COLUMN_TILE_H - 1;

const int apronStart = tileStart - 2;

const int apronEnd = tileEnd + 2;

// Clamp tile and apron limits by image borders

const int tileEndClamped = min(tileEnd, dataH - 1);

const int apronStartClamped = max(apronStart, 0);

const int apronEndClamped = min(apronEnd, dataH - 1);

// Current column index

const int columnStart = IMUL(blockIdx.x, COLUMN_TILE_W) + threadIdx.x;

const int num_WINDOW_ELEMENTS = IMUL(WINDOW_RADIUS+1+WINDOW_RADIUS,

WINDOW_RADIUS+1+WINDOW_RADIUS);

// Shared and global memory indices for current column

int smemPos = IMUL(threadIdx.y, COLUMN_TILE_W) + threadIdx.x;

int gmemPos = IMUL(apronStart + threadIdx.y, dataW) + columnStart;

// Cycle through the entire data cache

// Load global memory values, if indices are within the image borders,

// or initialize with zero otherwise

for(int y = apronStart + threadIdx.y; y < = apronEnd; y + = blockDim.y){

 rowMeans[smemPos] =

 ((y > = apronStartClamped) (y < =  apronEndClamped)) ?

 d_MeansRow[gmemPos] : 0;

 rowSquares[smemPos] =

 ((y > = apronStartClamped) (y < =  apronEndClamped)) ?

 d_SquaresRow[gmemPos] : 0;

 smemPos + = smemStride;

 gmemPos + = gmemStride;

}

// Ensure the completeness of the loading stage

// because results, emitted by each thread depend on the data,

// loaded by another threads

__syncthreads();

// Shared and global memory indices for current column

smemPos = IMUL(threadIdx.y + WINDOW_RADIUS, COLUMN_TILE_W) + threadIdx.x;

gmemPos = IMUL(tileStart + threadIdx.y, dataW) + columnStart;

// Cycle through the tile body, clamped by image borders

// Calculate the sum of the strips sums (and squares)

for(int y = tileStart + threadIdx.y; y < = tileEndClamped; y + = blockDim.y){

 int sumMeansFull = 0;

 int sumSquaresFull = 0;

// Loop unrolling section

// only uses unrolling templates (defined at beginning) if defined

// if not defined, uses for loop

#ifdef UNROLL_INNER

 sumMeansFull = convolutionColumnFull5<2 *  2>(rowMeans + smemPos);

 sumSquaresFull = convolutionColumnFull5<2 * 2>  (rowSquares + smemPos);

#else

for(int k = -2; k < = 2; k++){

 sumMeansFull + =

 rowMeans[smemPos + IMUL(k, COLUMN_TILE_W)];

 sumSquaresFull + =

 rowSquares[smemPos + IMUL(k, COLUMN_ TILE_W)];

}

#endif

 // Write sums to device memory

 d_MeansFull[gmemPos] = sumMeansFull;

 d_SquaresFull[gmemPos] = sumSquaresFull;

 // Calcuation of Speckle Contrast (K)

 // K = standard deviation / mean

 // Compute mean

 float mean = sumMeansFull/(float)num_WINDOW_ELEMENTS;

 // Compute standard deviation

 float coeff = __fdividef(1,num_WINDOW_ELEMENTS-1);

 float SD = sqrtf(coeff*(sumSquaresFull -

 (num_WINDOW_ELEMENTS*mean*mean)));

 float K = SD/mean;

 // Calculation of Speckle Flow Index (SFI)

 // SFI = 1 / (2 * tc * K * K)

 // tc = correlation time = 10ms

 float SFI = 1 / (2 * expT * K * K);

 // move calculated values into master Speckle Contrast matrix

 d_SC[gmemPos] = K;

 // move calculated values into master Speckle Flow matrix

 d_SFI[gmemPos] = SFI;

 smemPos + = smemStride;

 gmemPos + = gmemStride;

 }

}

Appendix B

////////////////////////////////////////////////////////////////////////////////////////

// ‘Roll’ Algorithm, Custom Integration based off Tom et al (Ref. 17)

////////////////////////////////////////////////////////////////////////////////////////

// width and height are the width and height of the raw image, respecively

// Raw is the raw speckle image

// SpeckleContrast and SpeckleFlowIndex are preallocated 1D matricies

// rollRow, rollColumn, rollRowSquares, rollColumnSquared are preallocated

// buffers for holding the strip sums in the row and column and strip sums

// of squared values in the row and column, respectively

// the size of these 4 buffers are width*1

// wR is the size of the sliding window radius

// t is the exposure time of the camera used to obtain the raw speckle images

void LSI_roll(int *Raw, float *SpeckleContrast, float *SpeckleFlowIndex,

 int *rollRow, int *rollColumn,

 int *rollRowSquared, int *rollColumnSquared,

 int width, int height,

 int wR, float t)

{

// Full window size

int w = wR*2 + 1;

// Number of elements within sliding window

int els = w*w;

// Inverse degrees of freedom

float iDoF = (float)1.0/(els*(els-1));

// Intialize counters

int i, j;

// Step 1) Calculate first accumulated sum to start

// Perform analysis on first accumulated row

for (i = 0; i<width; i++)

{

for (j = 0; j<w; j++)

{

rollRow[i] + = Raw[j*width+i];

rollRowSquared[i] + = (Raw[j*width+i]*Raw[j*width+i]);

}

if (i> = w)

{

rollColumn[i] = rollColumn[i-1] - rollRow[i-w] + rollRow[i];

rollColumnSquared[i] = rollColumnSquared[i-1] - rollRowSquared[i-w] +

rollRowSquared[i];

SpeckleContrast[wR*width+i-wR] = (els * sqrt((float)iDoF*(els*

rollColumnSquared[i]-(rollColumn[i]*rollColumn[i]))))/rollColumn[i];

SpeckleFlowIndex[wR*width+i-wR] = (float)1/(2*SpeckleContrast[wR*width+i-wR]*

SpeckleContrast[wR*width+i-wR]*t);

}

else

{

rollColumn[w-1] + = rollRow[i];

rollColumnSquared[w-1] + = rollRowSquared[i];

if(i = = (w-1))

{

SpeckleContrast[wR*width+i-wR] = (els * sqrt((float)iDoF* (els*

rollColumnSquared[w-1]-(rollColumn[w-1]*rollColumn[w-1]))))/rollColumn[w-1];

SpeckleFlowIndex[wR*width+i-wR] = (float)1/(2*SpeckleContrast[wR*width+iwR]*

SpeckleContrast[wR*width+i-wR]*t);

}

}

}

// Step 2) Perform optimized roll algorithm

// Cumulate rows and columns simultaneously

// Calculate SC and SFI immediately

for (j = w; j<height; j++)

{

rollColumn[w-1] = 0;

rollColumnSquared[w-1] = 0;

for (i = 0; i<width; i++)

{

rollRow[i] = rollRow[i] - Raw[(j-w)*width+i] + Raw[j*width+i];

rollRowSquared[i] = rollRowSquared[i] - (Raw[(j-w)*width+i]*Raw[(j-w)*width+

i]) + (Raw[j*width+i]*Raw[j*width+i]);

if (i> = w)

{

rollColumn[i] = rollColumn[i-1] - rollRow[i-w] + rollRow[i];

rollColumnSquared[i] = rollColumnSquared[i-1] –

rollRowSquared[i-w] + rollRowSquared[i];

SpeckleContrast[(j-wR)*width+i-wR] = (els * sqrt((float)iDoF*(els*

rollColumnSquared[i]-(rollColumn[i]* rollColumn[i]))))/rollColumn[i];

SpeckleFlowIndex[(j-wR)*width+i-wR] = (float)1/(2*SpeckleContrast[(j-wR)*

width+i-wR]*SpeckleContrast[(j-wR)*width+i-wR]*t);

}

else

{

rollColumn[w-1] + = rollRow[i];

rollColumnSquared[w-1] + = rollRowSquared[i];

if(i = = (w-1))

{

SpeckleContrast[(j-wR)*width+i-wR] = (els * sqrt((float)iDoF*(els*

rollColumnSquared[w-1]-(rollColumn[w-1]*rollColumn[w-1]))))/rollColumn[w-1];

SpeckleFlowIndex[(j-wR)*width+i-wR] = (float)1/(2*SpeckleContrast[(jwR)*

width+i-wR]*SpeckleContrast[(j-wR)*width+i-wR]*t);

}

}

}

}

}