# 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 Convolution(VertexShaderOutput input, float3x3 kernel)
{
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.

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
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)? 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.