# problems generating signed distance field

## Recommended Posts

Hi! I'm trying to implemente Valve's paper about improved alpha testing (http://www.valvesoftware.com/publications/2007/SIGGRAPH2007_AlphaTestedMagnification.pdf) but unfortunately I'm not getting the results I was expecting. Have a look at the following image:

[attachment=35993:Screenshot_3.png]

I borrowed the image from here [http://www.codersnotes.com/notes/signed-distance-fields/]. At the right we have the source image. At the left we have the expected result. And in the middle my results. At first look it may seem as if the middle image is just overall darker than the correct one, but if you look closely you will notice some harsh transitions where the edge (the 'middle value', 0.5) should be. I can't really tell if it's a problem of how the value mapping is being done, or how the distance are being calculated. Below I leave a few of the key functions that I'm using:

Here is the main loop, InputData is an array of floats representing the grey values of an image (monochrome single channel).

void GenerateSDF()
{
for(int i=0; i< InputData.Length; i++)
{

int y = i/image_width;
int x = i % image_width;
float dist = findClosest(x,y);

//this writes the computed distance to an array of floats
set_output(x,y, dist);
}

}


Here is the function that calculates the minimum distance to a texel of the oposite color (as the paper indicates)

float findClosest(int current_x, int current_y)
{

float min_distance = 9999f;

//sample_raw_image just samples the InputData array, converting <x,y> coordinates into array index.
bool texel_outside = sample_raw_image(current_x, current_y) < 0.5f;
// texel_outside = true     = black
// texel_outside = false    = white

for(int i=0; i< image_width * image_height; i++)
{
int y = i /image_width;
int x = i % image_width;

// This function determines if the current texel being processed is of the oposite color
// than the reference pixel from the main loop. If it is then we calculate the nearest distance.
if(texel_is_oposite(texel_outside, x,y))
{
min_distance = get_min_distance(current_x, current_y, x, y, min_distance);
}

}

// if the texel is outside (is black) of the shape  then we use positive distance.
return texel_outside? min_distance : -min_distance;
}


For clarity I also leave here the texel_is_oposite and the get_min_distance functions:

    bool texel_is_oposite(bool texel_is_outside, int x, int y)
{
if(texel_is_outside && sample_raw_image(x,y)> 0.5f) //outside = black
return true;
if(!texel_is_outside && sample_raw_image(x,y)< 0.5f)//inside = white
return true;

return false;
}
float get_min_distance(float x1, float y1,  float x2, float y2, float current_min)
{
float x = x2 - x1;
float y = y2 - y1;
float length = (float)Math.Sqrt(x * x + y * y);
return Math.Min(length, current_min);

}


Then after the output array is calculated, I proceed to calculte the minimum and maximum values of the entire output array, and map the values from the range minimum..maximum to -1..1 and then 0..1

float[] remap(float[] values)
{
float minimum = 9999f;
float maximum = -9999f;
float[] new_values = new float[values.Length];

for(int i=0; i< values.Length; i++)
{
minimum = Mathf.Min(values[i], minimum);
maximum = Mathf.Max(values[i], maximum);
}

for(int i=0; i< values.Length; i++)
{
// map converts a value from an initial range to a target range,
// in this case from minimum..maximum to -1..1
new_values[i] = map(minimum, maximum, -1f, 1f, values[i]) * 0.5f + 0.5f;
}
return new_values;
}


All that remains is to create an image with that array and display it on screen. Is there anything wrong with the code? Any help is greatly apreciated.

Also note that I'm using a brute-force approach, this is mainly because I want to understand with my own hands how the algorithm works, that's why I haven't used the algorithm described in the link I posted which is more efficient.

Cheers

Edited by ramiro_fages

##### Share on other sites
to visualize it for your eyes, you might need to map your image to sRGB.

##### Share on other sites
You need to develop debugging skills.

Pick a particular pixel for which the output doesn't match the reference output. Then compute by hand what you think the output should be. If your manual computation matches the reference code, the problem is with your code; otherwise, the problem is with your understanding of the paper (or with the paper itself). Go from there. Edited by Ã�lvaro

##### Share on other sites

Your calculations seem correct. Definitely sub-optimal, but ok for a "reference" implementation as you mentioned.

I think the problem is with the final step where it converts the distance data to the bitmap for visualisation.

Your technique scales the entire distance field based on the maximum and minimum distances.  I've found that clamping the values to a certain range produces better results.

For SDF's the important information is generally around 0, so I kept my distance values to the range (-64 to 64 pixels), anything higher or lower is clamped. Then you map that range of distances to your grayscale.

Another thing I just thought of, is that your algorithm may not permit a minimum range of 0, since its always looking for texels that are opposite to the current one. Therefore the minimum distance may always be 1, which would explain the discontinuity as the gradient would go directly from -1 to 1.

Edited by Postie

##### Share on other sites

Try this code that implements the very fast 4SED algorithm. This code computes the SDF for a large (e.g. 1024x1024) image in tens of milliseconds. It takes a single-channel grayscale image, applies a threshold operation to compute a binary image, then applies the 4SED algorithm twice to compute the distance for the foreground/background. The final image is the subtraction of the two fields, that is then normalized and clamped to a maximum distance range. The only part missing would be conversion to fixed point [0,255], if desired.

typedef math::Vector2f Offset;

Bool DistanceFilter:: processFrame( const ImageFrame& inputFrame, ImageFrame& outputFrame )
{
if ( inputFrame.getImageCount() == 0 || outputFrame.getImageCount() == 0 )
return false;

const ImageBuffer* inputImage = inputFrame.getImage(0);
ImageBuffer* outputImage = outputFrame.getImage(0);

if ( inputImage == NULL || outputImage == NULL )
return false;

const PixelFormat& inputPixelFormat = inputImage->getPixelFormat();

// Make sure the input image is 2D and has one channel.
if ( inputImage->getDimensionCount() != 2 || inputPixelFormat.getChannelCount() != 1 )
return false;

// Determine the size of the output image.
const Size2D inputSize = inputImage->getSize2D();
const Size2D outputSize = inputSize + padding*2;
PixelFormat outputPixelFormat( ColorSpace::LINEAR_GRAY, ScalarType(ScalarType::FLOAT32) );
outputImage->setFormat( outputPixelFormat, outputSize.x, outputSize.y );
Float32* const outputPixels = outputImage->getPixels();

// Allocate temporary storage.
const Size outputPixelCount = outputSize.x*outputSize.y;
offsetMap.allocate( outputPixelCount );
Offset* const offsetPixels = offsetMap.getPointer();

//********************************************************************************
// Compute the distance map for the background area.

computeDistanceMap( inputImage->getPixels(), inputSize, padding, threshold, offsetPixels, outputSize, invert );

// Distance evaluation.
{
Offset* offset = offsetPixels;
const Offset* const offsetEnd = offsetPixels + outputPixelCount;
Float32* output = outputPixels;

for ( ; offset != offsetEnd; offset++, output++ )
*output = (*offset).getMagnitude();
}

//********************************************************************************
// Compute the distance map for the foreground area.

if ( signedDistance )
{
computeDistanceMap( inputImage->getPixels(), inputSize, padding, threshold, offsetPixels, outputSize, !invert );

// Distance evaluation.
{
Offset* offset = offsetPixels;
const Offset* const offsetEnd = offsetPixels + outputPixelCount;
Float32* output = outputPixels;

for ( ; offset != offsetEnd; offset++, output++ )
*output -= (*offset).getMagnitude();
}
}

//********************************************************************************
// Normalization.

if ( normalize )
{
Float32 maxDistance = range;

if ( maxDistance == 0.0f )
{
maxDistance = math::max( math::abs(math::min( outputPixels, outputPixelCount )),
math::abs(math::max( outputPixels, outputPixelCount )) );
}

if ( maxDistance != 0.0f )
{
const Float32 outputCenter = signedDistance ? outputThreshold : 0.0f;
const Float32 invRange = (Float32(1.0) - outputCenter)/maxDistance;

Float32* output = outputPixels;
const Float32* const outputEnd = outputPixels + outputPixelCount;

for ( ; output != outputEnd; output++ )
*output = ((*output))*invRange + outputCenter;
}
}

return true;
}
/**
* This implementation uses the 4SED distance mapping algorithm proposed in:
* Danielsson, P. "Euclidean Distance Mapping" (1980)
*
* The algorithm is slightly modified to handle special cases at the image boundaries,
* where the original algorithm would produce incorrect results if the foreground touched
* the edge of the image.
*
* TODO: Implement the corrections proposed in:
* Cuisenaire, O. and Macq, B.
* "Fast and Exact Signed Euclidean Distance Transformation with Linear Complexity" (1999)
*/
void DistanceFilter:: computeDistanceMap( const Float32* inputPixels, Size2D inputSize, Size2D inputPaddding,
Float32 inputThreshold, Offset* offsetPixels, Size2D outputSize, Bool invert )
{
const Size outputPixelCount = outputSize.x*outputSize.y;

//********************************************************************************
// Initialize the offset image.

const Float32 maxPossibleDistance = math::Vector2f(outputSize).getMagnitude();
const Float32 backgroundInitial = invert ? 0.0f : maxPossibleDistance;
const Float32 foregroundInitial = invert ? maxPossibleDistance : 0.0f;

{
Offset* offset = offsetPixels;
const Offset* const offsetInputEnd = offsetRowPaddingEnd + outputSize.x*inputSize.y;
const Offset* const offsetEnd = offsetPixels + outputPixelCount;
const Float32* input = inputPixels;

// The padding area above the image.
for ( ; offset != offsetRowPaddingEnd; offset++ )
*offset = Offset(backgroundInitial);

// The rows containing the input image.
while ( offset != offsetInputEnd )
{
const Offset* const inputRowEnd = paddingEnd + inputSize.x;
const Offset* const offsetRowEnd = offset + outputSize.x;

for ( ; offset != paddingEnd; offset++ )
*offset = Offset(backgroundInitial);

for ( ; offset != inputRowEnd; offset++, input++ )
*offset = Offset( (*input < inputThreshold) ? backgroundInitial : foregroundInitial );

for ( ; offset != offsetRowEnd; offset++ )
*offset = Offset(backgroundInitial);
}

// The padding area below the image.
for ( ; offset != offsetEnd; offset++ )
*offset = Offset(backgroundInitial);
}

//********************************************************************************
// First scan.

for ( Index y = 0; y < outputSize.y; y++ )
{
if ( y == 0 )
{
// Handle first row outside bounds of image.
for ( Index x = 0; x < outputSize.x; x++ )
{
Offset& L_xy = offsetPixels[y*outputSize.x + x];
Offset Lshift = Offset(backgroundInitial) + Offset(0,1);
if ( Lshift.getMagnitudeSquared() < L_xy.getMagnitudeSquared() )
L_xy = Lshift;
}
}
else
{
for ( Index x = 0; x < outputSize.x; x++ )
{
Offset& L_xy = offsetPixels[y*outputSize.x + x];
Offset Lshift = offsetPixels[(y-1)*outputSize.x + x] + Offset(0,1);
if ( Lshift.getMagnitudeSquared() < L_xy.getMagnitudeSquared() )
L_xy = Lshift;
}
}

// first column
{
Offset& L_xy = offsetPixels[y*outputSize.x + 0];
Offset Lshift = Offset(backgroundInitial) + Offset(1,0);
if ( Lshift.getMagnitudeSquared() < L_xy.getMagnitudeSquared() )
L_xy = Lshift;
}

for ( Index x = 1; x < outputSize.x; x++ )
{
Offset& L_xy = offsetPixels[y*outputSize.x + x];
Offset Lshift = offsetPixels[y*outputSize.x + (x-1)] + Offset(1,0);
if ( Lshift.getMagnitudeSquared() < L_xy.getMagnitudeSquared() )
L_xy = Lshift;
}

// last column
{
Offset& L_xy = offsetPixels[y*outputSize.x + (outputSize.x-1)];
Offset Lshift = Offset(backgroundInitial) + Offset(1,0);
if ( Lshift.getMagnitudeSquared() < L_xy.getMagnitudeSquared() )
L_xy = Lshift;
}

for ( Index x = outputSize.x - 2; x < outputSize.x; x-- )
{
Offset& L_xy = offsetPixels[y*outputSize.x + x];
Offset Lshift = offsetPixels[y*outputSize.x + (x+1)] + Offset(1,0);
if ( Lshift.getMagnitudeSquared() < L_xy.getMagnitudeSquared() )
L_xy = Lshift;
}
}

//********************************************************************************
// Second scan.

const Index lastRow = outputSize.y - 1;

for ( Index y = lastRow; y < outputSize.y; y-- )
{
if ( y == lastRow )
{
// Handle last row outside bounds of image.
for ( Index x = 0; x < outputSize.x; x++ )
{
Offset& L_xy = offsetPixels[y*outputSize.x + x];
Offset Lshift = Offset(backgroundInitial) + Offset(0,1);
if ( Lshift.getMagnitudeSquared() < L_xy.getMagnitudeSquared() )
L_xy = Lshift;
}
}
else
{
for ( Index x = 0; x < outputSize.x; x++ )
{
Offset& L_xy = offsetPixels[y*outputSize.x + x];
Offset Lshift = offsetPixels[(y+1)*outputSize.x + x] + Offset(0,1);
if ( Lshift.getMagnitudeSquared() < L_xy.getMagnitudeSquared() )
L_xy = Lshift;
}
}

// first column
{
Offset& L_xy = offsetPixels[y*outputSize.x + 0];
Offset Lshift = Offset(backgroundInitial) + Offset(1,0);
if ( Lshift.getMagnitudeSquared() < L_xy.getMagnitudeSquared() )
L_xy = Lshift;
}

for ( Index x = 1; x < outputSize.x; x++ )
{
Offset& L_xy = offsetPixels[y*outputSize.x + x];
Offset Lshift = offsetPixels[y*outputSize.x + (x-1)] + Offset(1,0);
if ( Lshift.getMagnitudeSquared() < L_xy.getMagnitudeSquared() )
L_xy = Lshift;
}

// last column
{
Offset& L_xy = offsetPixels[y*outputSize.x + (outputSize.x-1)];
Offset Lshift = Offset(backgroundInitial) + Offset(1,0);
if ( Lshift.getMagnitudeSquared() < L_xy.getMagnitudeSquared() )
L_xy = Lshift;
}

for ( Index x = outputSize.x - 2; x < outputSize.x; x-- )
{
Offset& L_xy = offsetPixels[y*outputSize.x + x];
Offset Lshift = offsetPixels[y*outputSize.x + (x+1)] + Offset(1,0);
if ( Lshift.getMagnitudeSquared() < L_xy.getMagnitudeSquared() )
L_xy = Lshift;
}
}
}

Edited by Aressera

##### Share on other sites

Looks like a gamma problem, e.g., your algorithm expects linear pixel data and you are using non-linear sRGB pixel data.

##### Share on other sites

Thanks to everyone who commented on this thread. Here are my findings:

1) Applying gamma correction (as Krypt0n and BFG suggested) fixed the brightness issue, but as you can se below, the "harsh edge" issue still persisted.
[attachment=36011:Screenshot_1.png]

2) Thanks to Alvaro's suggestion I decided to make the example case easier, and tested with an image of 20x1 pixels, with the first 10 pixels white, and printed the integer distances before doing the range map. I noticed that (as Postie mentioned) I was having trouble with the distances near 0, because the distances were jumping from -1 to 1 without going through 0. I fixed it this way ( I don't know if it's the proper way to do it...but it works  :P).

    float findClosest(int current_x, int current_y)
{

int min_distance = 9999;
bool texel_outside = sample_raw_image(current_x, current_y) < 0.5f;

for(int i=0; i< image_width * image_height; i++)
{
int y = i /image_width;
int x = i % image_width;

if(texel_is_oposite(texel_outside, x,y))
{
min_distance = get_min_distance(current_x, current_y, x, y, min_distance);
}

}

// FIX :)
return texel_outside? min_distance -1 : -min_distance;
}


This yielded the following results:

[attachment=36012:face 256 SDF.png]

Unfortunately as you can see, it's not the same as the reference image, and when rendering I have to use a value like 0.6 as the cutoff instead of 0.5 as the paper mentions. Having said that, it kinda works. I will give it a few more tries and then will move to a better implementation, as the algorithm that Aressera mentioned (thanks for it! ).

Cheers!

Edited by ramiro_fages

##### Share on other sites

Your remapping from minimum/maximum is going to shift the middle point, because it's not symmetrical.
e.g. if your edge distance is 0, your maximum distance is 10 and your minimum distance is -1, then after remapping:
min = 0.0 = 0/255
edge = ~0.09 = 23/255 -- should be 0.5 or 127.5/255!!!
max = 1.0 = 255/255

You could make it symmetrical with something like:

        for(int i=0; i< values.Length; i++)
{
minimum = min(values[i], minimum);
maximum = max(values[i], maximum);
}
float m = max(abs(minimum), abs(maximum));
maximum = m;
minimum = -m;

##### Share on other sites

Oh you're right, I can't believe I didn't notice it earlier. Unfortunately that seems to 'bias' things to one side or the other depending on which side has higher range than the other. Here's the result with your code:

[attachment=36028:carita 256 SDF_5.png]

But I can see that there's definitely a problem with how I'm mapping the values. Unfortunately the paper makes no mention about this at all, so I'll keep trying. Thanks for your time!

Edited by ramiro_fages

## Create an account

Register a new account