# How to handle floating point operations in HLSL?

I'm trying to write a perona malik anisotropic diffusion filter for the GPU.

I'm basing my shader off a matlab implementation of the filter. I'm running into trouble because of what I suspect is floating point arithmetic issues.

I have a line in my shader

float cN = 1 / pow( 1 + (nablaN / k), 2);

cN always evaluates to zero, whereas the matlab version does not. I think I am mishandling floating point numbers.

How should I handle this properly?

EDIT

Nabla is a float computed by convolution of the texture with a kernel hN.

float4 nablaN = Convolution(psInput, hN);

The kernel is

static const float3x3 hN =
{
0,  1, 0,
0, -1, 0,
0,  0, 0
};

and convolution is

{
float4 pixel = float4(0.0f, 0.0f, 0.0f, 0.0f);

for (int i = -1; i <= 1; ++i)
{
for (int j = -1; j <= 1; ++j)
{
pixel += kernel[i+1][j+1] * tex2D(Input0Sampler, input.TextureCoordinate + float2(i,j));
};
};

return pixel;
}

EDIT

If I return nablaN from the pixel shader I get an image like this

http://imageshack.us/photo/my-images/819/dw69.png/ (nablaN/2 is similar)

If i return ( 1 + (nablaN / k)) * ( 1 + (nablaN / k)) from the pixel shader I get

something looking like this http://imageshack.us/photo/my-images/822/s49z.png/

This looks similar to an edge detection filter. Most of the pixel values are zero. So you can see why 1/(( 1 + (nablaN / k)) * ( 1 + (nablaN / k)) ) just gives me a black screen.

EDIT

OK I have three versions of the program now in three languages

FREEMAT

function diff_im = anisodiff2D(im, num_iter, delta_t, kappa, option)
%ANISODIFF2D Conventional anisotropic diffusion
%   DIFF_IM = ANISODIFF2D(IM, NUM_ITER, DELTA_T, KAPPA, OPTION) perfoms
%   conventional anisotropic diffusion (Perona & Malik) upon a gray scale
%   image. A 2D network structure of 8 neighboring nodes is considered for
%   diffusion conduction.
%
%       ARGUMENT DESCRIPTION:
%               IM       - gray scale image (MxN).
%               NUM_ITER - number of iterations.
%               DELTA_T  - integration constant (0 <= delta_t <= 1/7).
%                          Usually, due to numerical stability this
%                          parameter is set to its maximum value.
%               KAPPA    - gradient modulus threshold that controls the conduction.
%               OPTION   - conduction coefficient functions proposed by Perona & Malik:
%                          1 - c(x,y,t) = exp(-(nablaI/kappa).^2),
%                              privileges high-contrast edges over low-contrast ones.
%                          2 - c(x,y,t) = 1./(1 + (nablaI/kappa).^2),
%                              privileges wide regions over smaller ones.
%
%       OUTPUT DESCRIPTION:
%                DIFF_IM - (diffused) image with the largest scale-space parameter.
%
%   Example
%   -------------
%   s = phantom(512) + randn(512);
%   num_iter = 15;
%   delta_t = 1/7;
%   kappa = 30;
%   option = 2;
%   figure, subplot 121, imshow(s,[]), subplot 122, imshow(ad,[])
%

% References:
%   P. Perona and J. Malik.
%   Scale-Space and Edge Detection Using Anisotropic Diffusion.
%   IEEE Transactions on Pattern Analysis and Machine Intelligence,
%   12(7):629-639, July 1990.
%
%   G. Grieg, O. Kubler, R. Kikinis, and F. A. Jolesz.
%   Nonlinear Anisotropic Filtering of MRI Data.
%   IEEE Transactions on Medical Imaging,
%   11(2):221-232, June 1992.
%
%   MATLAB implementation based on Peter Kovesi's anisodiff(.):
%   P. D. Kovesi. MATLAB and Octave Functions for Computer Vision and Image Processing.
%   School of Computer Science & Software Engineering,
%   The University of Western Australia. Available from:
%   <http://www.csse.uwa.edu.au/~pk/research/matlabfns/>.
%
% Credits:
% Daniel Simoes Lopes
% ICIST
% Instituto Superior Tecnico - Universidade Tecnica de Lisboa
% danlopes (at) civil ist utl pt
% http://www.civil.ist.utl.pt/~danlopes
%
% May 2007 original version.

% Convert input image to double.
im = double(im);

% PDE (partial differential equation) initial condition.
diff_im = im;

% Center pixel distances.
dx = 1;
dy = 1;
dd = sqrt(2);

% 2D convolution masks - finite differences.
hN = [0 1 0; 0 -1 0; 0 0 0];
hS = [0 0 0; 0 -1 0; 0 1 0];
hE = [0 0 0; 0 -1 1; 0 0 0];
hW = [0 0 0; 1 -1 0; 0 0 0];
hNE = [0 0 1; 0 -1 0; 0 0 0];
hSE = [0 0 0; 0 -1 0; 0 0 1];
hSW = [0 0 0; 0 -1 0; 1 0 0];
hNW = [1 0 0; 0 -1 0; 0 0 0];

nablaN = zeros(1072,1912);
nablaS = zeros(1072,1912);
nablaW = zeros(1072,1912);
nablaE = zeros(1072,1912);
nablaNE = zeros(1072,1912);
nablaSE = zeros(1072,1912);
nablaSW = zeros(1072,1912);
nablaNW = zeros(1072,1912);

% Anisotropic diffusion.
for t = 1:num_iter

% Finite differences. [imfilter(.,.,'conv') can be replaced by conv2(.,.,'same')]
% nablaN = conv2(diff_im,hN)(2:1081,2:1921);
nablaN = conv2(diff_im,hN)(2:1073,2:1913);
nablaS = conv2(diff_im,hS)(2:1073,2:1913);
nablaW = conv2(diff_im,hW)(2:1073,2:1913);
nablaE = conv2(diff_im,hE)(2:1073,2:1913);
nablaNE = conv2(diff_im,hNE)(2:1073,2:1913);
nablaSE = conv2(diff_im,hSE)(2:1073,2:1913);
nablaSW = conv2(diff_im,hSW)(2:1073,2:1913);
nablaNW = conv2(diff_im,hNW)(2:1073,2:1913);

% Diffusion function.
if option == 1
cN = exp(-(nablaN/kappa).^2);
cS = exp(-(nablaS/kappa).^2);
cW = exp(-(nablaW/kappa).^2);
cE = exp(-(nablaE/kappa).^2);
cNE = exp(-(nablaNE/kappa).^2);
cSE = exp(-(nablaSE/kappa).^2);
cSW = exp(-(nablaSW/kappa).^2);
cNW = exp(-(nablaNW/kappa).^2);
elseif option == 2
cN = 1./(1 + (nablaN/kappa).^2);
cS = 1./(1 + (nablaS/kappa).^2);
cW = 1./(1 + (nablaW/kappa).^2);
cE = 1./(1 + (nablaE/kappa).^2);
cNE = 1./(1 + (nablaNE/kappa).^2);
cSE = 1./(1 + (nablaSE/kappa).^2);
cSW = 1./(1 + (nablaSW/kappa).^2);
cNW = 1./(1 + (nablaNW/kappa).^2);
end

% Discrete PDE solution.
diff_im = diff_im + ...
delta_t*(...
(1/(dy^2))*cN.*nablaN + (1/(dy^2))*cS.*nablaS + ...
(1/(dx^2))*cW.*nablaW + (1/(dx^2))*cE.*nablaE + ...
(1/(dd^2))*cNE.*nablaNE + (1/(dd^2))*cSE.*nablaSE + ...
(1/(dd^2))*cSW.*nablaSW + (1/(dd^2))*cNW.*nablaNW );

% Iteration warning.
fprintf('\rIteration %d\n',t);
end

HLSL

texture2D Input0;
sampler2D Input0Sampler = sampler_state
{
Texture = <Input0>;
MinFilter = Point;
MagFilter = Point;
MipFilter = Point;
};

{
float4 Position : POSITION0;
float2 TextureCoordinate : TEXCOORD0;
};

{
float4 Position : POSITION0;
float2 TextureCoordinate : TEXCOORD0;
};

{
// TODO: Optionally add/remove output indices to match GPUProcessor.numOutputs
float4 Index0 : COLOR0;
};

// input texture dimensions
static float w = 1920 - 8;
static float h = 1080 - 8;

static const float2 pixel = float2(1.0 / w, 1.0 / h);
static const float2 halfPixel = float2(pixel.x / 2, pixel.y / 2);

static const float3x3 hN =
{
0,  1, 0,
0, -1, 0,
0,  0, 0
};
static const float3x3 hS =
{
0,  0, 0,
0, -1, 0,
0,  1, 0
};
static const float3x3 hE =
{
0,  0, 0,
0, -1, 1,
0,  0, 0
};
static const float3x3 hW =
{
0,  0, 0,
1, -1, 0,
0,  0, 0
};
static const float3x3 hNE =
{
0,  0, 1,
0, -1, 0,
0,  0, 0
};
static const float3x3 hSE =
{
0,  0, 0,
0, -1, 0,
0,  0, 1
};
static const float3x3 hSW =
{
0,  0, 0,
0, -1, 0,
1,  0, 0
};
static const float3x3 hNW =
{
1,  0, 0,
0, -1, 0,
0,  0, 0
};

{

//output.Position = vsInput.Position;
//output.TextureCoordinate = vsInput.TextureCoordinate;

vsInput.Position.x =  vsInput.Position.x - 2*halfPixel.x;
vsInput.Position.y =  vsInput.Position.y + 2*halfPixel.y;
output.Position = vsInput.Position;
output.TextureCoordinate = vsInput.TextureCoordinate;
return output;

//return output;
}

{
//float4 pixel = float4(0.0f, 0.0f, 0.0f, 0.0f);
float pixel = 0.0L;

for (int i = -1; i <= 1; ++i)
{
for (int j = -1; j <= 1; ++j)
{
pixel += (kernel[i+1][j+1] * tex2D(Input0Sampler, input.TextureCoordinate + float2(i,j))).x;
};
};

return pixel;
}

{
//output.Index0 = float4(0,0,0,0);//tex2D(Input0Sampler, psInput.TextureCoordinate);
//output.Index0 = tex2D(Input0Sampler, psInput.TextureCoordinate);

float big = 100000.0f;
float total = tex2D(Input0Sampler, psInput.TextureCoordinate);
float dx, dy, dd;
dx = 1; dy = 1; dd = pow(2.0L, 0.5L);
float delta_t = 1/7;//0.25f;
float k = 30.0;//0.015f;

float nablaN = Convolution(psInput, hN);
float nablaS = Convolution(psInput, hS);
float nablaW = Convolution(psInput, hW);
float nablaE = Convolution(psInput, hE);
float nablaNE = Convolution(psInput, hNE);
float nablaSE = Convolution(psInput, hSE);
float nablaSW = Convolution(psInput, hSW);
float nablaNW = Convolution(psInput, hNW);

float cN  = ( 1 + (nablaN / k)) * ( 1 + (nablaN / k));
float cS  = ( 1 + (nablaS / k)) * ( 1 + (nablaS / k));
float cW  = ( 1 + (nablaW / k)) * ( 1 + (nablaW / k));
float cE  = ( 1 + (nablaE / k)) * ( 1 + (nablaE / k));
float cNE = ( 1 + (nablaNE / k)) * ( 1 + (nablaNE / k));
float cSE = ( 1 + (nablaSE / k)) * ( 1 + (nablaSE / k));
float cSW = ( 1 + (nablaSW / k)) * ( 1 + (nablaSW / k));
float cNW = ( 1 + (nablaNW / k)) * ( 1 + (nablaNW / k));

/*
float4 cN  = exp(-(nablaN/k)*(nablaN/k));
float4 cS  = exp(-(nablaS/k)*(nablaS/k));
float4 cW  = exp(-(nablaW/k)*(nablaW/k));
float4 cE  = exp(-(nablaE/k)*(nablaE/k));
float4 cNE = exp(-(nablaNE/k)*(nablaNE/k));
float4 cSE = exp(-(nablaSE/k)*(nablaSE/k));
float4 cSW = exp(-(nablaSW/k)*(nablaSW/k));
float4 cNW = exp(-(nablaNW/k)*(nablaNW/k));
*/

total += delta_t *
(
//(mul(cN, nablaN)/(dy*dy)) + (mul(cS, nablaS)/(dy*dy)) + (mul(cW, nablaW)/(dx*dx)) + (mul(cE, nablaE)/(dx*dx)) + (dd*dd)*(mul(cNE, nablaNE) + mul(cSE, nablaSE) + mul(cSW, nablaSW) + mul(cNW, nablaNW))
(nablaN/cN) + (nablaS/cS) + (nablaW/cW) + (nablaE/cE) + 2*((nablaNE/cNE) + (nablaSE/cSE) + (nablaSW/cSW) + (nablaNW/cNW))
);

output.Index0 = total;

return output;
}

technique PeronaMalik
{
pass pass1
{
}
}

OpenCL

float4 Convolution(__read_only image2d_t srcImg, int2 point, float * kern)
{
const sampler_t smp = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_LINEAR;
float4 sum = (float4)(0.0f,0.0f,0.0f,0.0f);
{
{
}
}
return sum;
}

__kernel void imagingTest(__read_only  image2d_t srcImg, __write_only image2d_t dstImg)
{
float k = 30;
float delta_t = 1/7;

float hN[9];
hN[0] = 0; hN[1] = 1; hN[2] = 0;
hN[3] = 0; hN[4] =-1; hN[5] = 0;
hN[6] = 0; hN[7] = 0; hN[8] = 0;

float hS[9];
hS[0] = 0; hS[1] = 0; hS[2] = 0;
hS[3] = 0; hS[4] =-1; hS[5] = 0;
hS[6] = 0; hS[7] = 1; hS[8] = 0;

float hE[9];
hE[0] = 0; hE[1] = 0; hE[2] = 0;
hE[3] = 0; hE[4] =-1; hE[5] = 1;
hE[6] = 0; hE[7] = 0; hE[8] = 0;

float hW[9];
hW[0] = 0; hW[1] = 0; hW[2] = 0;
hW[3] = 1; hW[4] =-1; hW[5] = 0;
hW[6] = 0; hW[7] = 0; hW[8] = 0;

float hNE[9];
hNE[0] = 0; hNE[1] = 0; hNE[2] = 1;
hNE[3] = 0; hNE[4] =-1; hNE[5] = 0;
hNE[6] = 0; hNE[7] = 0; hNE[8] = 0;

float hSE[9];
hSE[0] = 0; hSE[1] = 0; hSE[2] = 0;
hSE[3] = 0; hSE[4] =-1; hSE[5] = 0;
hSE[6] = 0; hSE[7] = 0; hSE[8] = 1;

float hSW[9];
hSW[0] = 0; hSW[1] = 0; hSW[2] = 0;
hSW[3] = 0; hSW[4] =-1; hSW[5] = 0;
hSW[6] = 1; hSW[7] = 0; hSW[8] = 0;

float hNW[9];
hNW[0] = 1; hNW[1] = 0; hNW[2] = 0;
hNW[3] = 0; hNW[4] =-1; hNW[5] = 0;
hNW[6] = 0; hNW[7] = 0; hNW[8] = 0;

const sampler_t smp = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_LINEAR;

int2 coord = (int2)(get_global_id(0), get_global_id(1));

uint4 bgra = read_imageui(srcImg, smp, coord);

float4 nablaN = Convolution(srcImg, coord, hN);
float4 nablaS = Convolution(srcImg, coord, hS);
float4 nablaE = Convolution(srcImg, coord, hE);
float4 nablaW = Convolution(srcImg, coord, hW);

float4 nablaNE = Convolution(srcImg, coord, hNE);
float4 nablaNW = Convolution(srcImg, coord, hNW);
float4 nablaSE = Convolution(srcImg, coord, hSE);
float4 nablaSW = Convolution(srcImg, coord, hSW);

float4 cN  = exp(-(nablaN/k)*(nablaN/k));
float4 cS  = exp(-(nablaS/k)*(nablaS/k));
float4 cW  = exp(-(nablaW/k)*(nablaW/k));
float4 cE  = exp(-(nablaE/k)*(nablaE/k));
float4 cNE = exp(-(nablaNE/k)*(nablaNE/k));
float4 cSE = exp(-(nablaSE/k)*(nablaSE/k));
float4 cSW = exp(-(nablaSW/k)*(nablaSW/k));
float4 cNW = exp(-(nablaNW/k)*(nablaNW/k));

float4 sum = delta_t * ((nablaN/cN) + (nablaS/cS) + (nablaW/cW) + (nablaE/cE) + 2*((nablaNE/cNE) + (nablaSE/cSE) + (nablaSW/cSW) + (nablaNW/cNW)));

bgra.x = bgra.y = bgra.z = convert_int(sum.x);

bgra.w = 255;
write_imageui(dstImg, coord, bgra);
}

It works fine in freemat, its just really slow. The OpenCL and HLSL versions are really fast but they don't give the right results.

• The matlab code is here mathworks.com.au/matlabcentral/fileexchange/… in case anyone is interested
– sav
Commented Jul 30, 2013 at 6:33
• The original image looks like this imageshack.us/photo/my-images/197/ema9.png
– sav
Commented Aug 2, 2013 at 1:42
• I can replace float with double in my OpenCL code and still get a black image as a result.
– sav
Commented Aug 12, 2013 at 5:19
• It turns out that there was trouble with how I was encoding pixel values in the image to begin with. Take a look at emgu.com/wiki/files/1.3.0.0/html/…
– sav
Commented Jun 11, 2014 at 2:09

Is nablaN an integer? Because if it is, then the expression 1 / pow( 1 + (nablaN / k), 2); is considered the division of an integer by an integer, which is going to be zero if nablaN is greater than 2. These are just the rules of the C programming language.

• nablaN is a float
– sav
Commented Aug 1, 2013 at 4:29
• Then print out the individual terms: what is the value of nablaN/2? What is the value of pow( 1 + (nablaN / k), 2)? Commented Aug 1, 2013 at 10:49

Looks like I have figured it out, at least for OpenCL ...

#ifdef cl_khr_fp64
#pragma OPENCL EXTENSION cl_khr_fp64 : enable
#elif defined(cl_amd_fp64)
#pragma OPENCL EXTENSION cl_amd_fp64 : enable
#else
#error "Double precision doubleing point not supported by OpenCL implementation."
#endif

double4 Convolution(__read_only image2d_t srcImg, int2 point, double * kern)
{
const sampler_t smp = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_LINEAR;
double4 sum = (double4)(0.0f,0.0f,0.0f,0.0f);
{
{
int2 pos = (int2)(i,j);
sum += kern[(delta.y*3) + delta.x] * convert_double4(read_imageui(srcImg, smp, point + pos));
}
}
return sum;
}

__kernel void imagingTest(__read_only  image2d_t srcImg, __write_only image2d_t dstImg)
{

double k = 30.0L;
double delta_t = 0.14285714285714285714285714285714L; // 1/7

double hN[9];
hN[0] = 0; hN[1] = 1; hN[2] = 0;
hN[3] = 0; hN[4] =-1; hN[5] = 0;
hN[6] = 0; hN[7] = 0; hN[8] = 0;

double hS[9];
hS[0] = 0; hS[1] = 0; hS[2] = 0;
hS[3] = 0; hS[4] =-1; hS[5] = 0;
hS[6] = 0; hS[7] = 1; hS[8] = 0;

double hE[9];
hE[0] = 0; hE[1] = 0; hE[2] = 0;
hE[3] = 0; hE[4] =-1; hE[5] = 1;
hE[6] = 0; hE[7] = 0; hE[8] = 0;

double hW[9];
hW[0] = 0; hW[1] = 0; hW[2] = 0;
hW[3] = 1; hW[4] =-1; hW[5] = 0;
hW[6] = 0; hW[7] = 0; hW[8] = 0;

double hNE[9];
hNE[0] = 0; hNE[1] = 0; hNE[2] = 1;
hNE[3] = 0; hNE[4] =-1; hNE[5] = 0;
hNE[6] = 0; hNE[7] = 0; hNE[8] = 0;

double hSE[9];
hSE[0] = 0; hSE[1] = 0; hSE[2] = 0;
hSE[3] = 0; hSE[4] =-1; hSE[5] = 0;
hSE[6] = 0; hSE[7] = 0; hSE[8] = 1;

double hSW[9];
hSW[0] = 0; hSW[1] = 0; hSW[2] = 0;
hSW[3] = 0; hSW[4] =-1; hSW[5] = 0;
hSW[6] = 1; hSW[7] = 0; hSW[8] = 0;

double hNW[9];
hNW[0] = 1; hNW[1] = 0; hNW[2] = 0;
hNW[3] = 0; hNW[4] =-1; hNW[5] = 0;
hNW[6] = 0; hNW[7] = 0; hNW[8] = 0;

const sampler_t smp = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_LINEAR;

int2 coord = (int2)(get_global_id(0), get_global_id(1));

uint4 bgra = read_imageui(srcImg, smp, coord);

double4 nablaN = Convolution(srcImg, coord, hN);
double4 nablaS = Convolution(srcImg, coord, hS);
double4 nablaE = Convolution(srcImg, coord, hE);
double4 nablaW = Convolution(srcImg, coord, hW);

double4 nablaNE = Convolution(srcImg, coord, hNE);
double4 nablaNW = Convolution(srcImg, coord, hNW);
double4 nablaSE = Convolution(srcImg, coord, hSE);
double4 nablaSW = Convolution(srcImg, coord, hSW);

double4 cN  = exp(-(nablaN /k) * (nablaN /k));
double4 cS  = exp(-(nablaS /k) * (nablaS /k));
double4 cW  = exp(-(nablaW /k) * (nablaW /k));
double4 cE  = exp(-(nablaE /k) * (nablaE /k));
double4 cNE = exp(-(nablaNE/k) * (nablaNE/k));
double4 cSE = exp(-(nablaSE/k) * (nablaSE/k));
double4 cSW = exp(-(nablaSW/k) * (nablaSW/k));
double4 cNW = exp(-(nablaNW/k) * (nablaNW/k));

double4 sum = 0.5 * (nablaNE * cNE) + (nablaSE * cSE) + (nablaSW * cSW) + (nablaNW * cNW);
sum += (nablaN * cN) + (nablaS * cS) + (nablaW * cW) + (nablaE * cE);
sum *= delta_t;

bgra.x = bgra.y = bgra.z = convert_int(sum.x);

bgra.w = 255;
write_imageui(dstImg, coord, bgra);
}

This seems to give good results. I'm not sure whether this is a limitation of HLSL.