Fractals/Computer graphic techniques/2D/algorithms


postprocessing = modification of the image = graphic algorithms = image processing

Algorithms by graphic type

  • raster algorithms
  • vector algorithms

Five classes of image processing algorithms:

  • image restoration
  • image analysis
  • image synthesis
  • image enhancement
  • image compression.


  • Morphological operations on binary images[6][7]
    • morphological closing = dilation followed by erosion
    • morphological opening = erosion followed by dilation


  • Two types of edge detection
  • Pseudo-3D projections
  • Star field generator
  • Random dot stereograms (aka Magic Eye)
  • Motion blur for animations
  • Interlacing
  • Embossing
  • Antialiasing
  • Palette emulation to allow color cycling on true-color displays
  • True color emulation that provides dithering on 256-color display




The image processing operators ( algorithms) can be classified into the 4 categories

  • pixel algorithms = pixel processing: the point operator acts on individual pixels : g = H(f)
    • add two images: g(x)=f1(x)+f2(x)
    • the histogram
    • look-up table (LUT)
    • windowing
  • local algorithms = kernaling = local processing: "A local operator calculates the pixel value g(x) based not only on the value in the same position in the input image, i.e. f(x) , but it used several values in nearby points f(x+y). Local operators are at the heart of almost all image processing tasks."[9]
    • spatial location filtering (convolution),
    • spatial frequency filtering using high- and low-pass digital filters,
    • the unsharp masking technique[10]
  • geometric processing = geometric transformations: "Whereas a point operator changes the value of all pixels a geometrical operator doesn’t change the value but instead it ‘moves’ a pixel to a new position."
  • Global Operators: "A global operator (potentially) needs al pixel values in the input image to calculate the value for just one pixel in the output image."


Digital topologyEdit


gamma correctionEdit

dense imageEdit

Dense image[11][12][13][14][15]

  • downscaling with gamma correction[16]
  • path finding[17]
  • aliasing [18]
  • changing algorithm ( representation function) from discrete to continous, like from level set method of escape time to continous ( DEM )
  "the denser the area, the more heavy the anti-aliasing have to be in order to make it look good."  knighty[19]
  "the details are smaller than pixel spacing, so all that remains is the bands of colour shift from period-doubling of features making it even denser"  claude[20]

path findingEdit

path finding





  • geometric operations for two-dimensional polygons[22][23]

How to tell whether a point is to the right or left side of a line ?Edit

  How to tell whether a point is to the right or left side of a line ?

  a, b = points
  line = ab
 pont to check = z

  position = sign((Bx - Ax) * (Y - Ay) - (By - Ay) * (X - Ax))
  It is 0 on the line, and +1 on one side, -1 on the other side.


double CheckSide(double Zx, double Zy, double Ax, double Ay, double Bx, double By)
  return ((Bx - Ax) * (Zy - Ay) - (By - Ay) * (Zx - Ax));


Testing if point is inside triangleEdit

c console program 
gcc t.c -Wall


# include <stdio.h>

// 3 points define triangle 
double Zax = -0.250000000000000;
double Zay = 0.433012701892219;
// left when y
double Zlx = -0.112538773749444;  
double Zly = 0.436719687479814 ;

double Zrx = -0.335875821657728;
double Zry = 0.316782798339332;

// points to test 
// = inside triangle 
double Zx = -0.209881783739630;
double Zy =   +0.4;

// outside triangle 
double Zxo = -0.193503885412548  ;
double Zyo = 0.521747636163664;

double Zxo2 = -0.338750000000000;
double Zyo2 = +0.440690927838329;

// ============
// In general, the simplest (and quite optimal) algorithm is checking on which side of the half-plane created by the edges the point is.
double sign (double  x1, double y1,  double x2, double y2, double x3, double y3)
    return (x1 - x3) * (y2 - y3) - (x2 - x3) * (y1 - y3);

int  PointInTriangle (double x, double y, double x1, double y1, double x2, double y2, double x3, double y3)
    double  b1, b2, b3;

    b1 = sign(x, y, x1, y1, x2, y2) < 0.0;
    b2 = sign(x, y, x2, y2, x3, y3) < 0.0;
    b3 = sign(x, y, x3, y3, x1, y1) < 0.0;

    return ((b1 == b2) && (b2 == b3));

int Describe_Position(double Zx, double Zy){
if (PointInTriangle( Zx, Zy, Zax, Zay, Zlx, Zly, Zrx, Zry))
  printf(" Z is inside \n");
  else printf(" Z is outside \n");

return 0;

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

int main(void){

Describe_Position(Zx, Zy);
Describe_Position(Zxo, Zyo);
Describe_Position(Zxo2, Zyo2);

return 0;

Orientation and area of the triangleEdit

Orientation and area of the triangle : how to do it ?

// gcc t.c -Wall
// ./a.out
# include <stdio.h>

double GiveTriangleArea(double xa, double ya, double xb, double yb, double xc, double yc)
return ((xb*ya-xa*yb)+(xc*yb-xb*yc)+(xa*yc-xc*ya))/2.0;


wiki Curve_orientation

The orientation of a triangle (clockwise/counterclockwise) is the sign of the determinant

matrix = { {1 , x1, y1}, {1 ,x2, y2} , {1,  x3, y3}}

(x_1,y_1), (x_2,y_2), (x_3,y_3)$ 
are the Cartesian coordinates of the three vertices of the triangle.

:<math>\mathbf{O} = \begin{bmatrix}

1 & x_{A} & y_{A} \\
1 & x_{B} & y_{B} \\
1 & x_{C} & y_{C}\end{bmatrix}.</math>

A formula for its determinant may be obtained, e.g., using the method of [[cofactor expansion]]:
\det(O) &= 1\begin{vmatrix}x_{B}&y_{B}\\x_{C}&y_{C}\end{vmatrix}
+y_{A}\begin{vmatrix}1&x_{B}\\1&x_{C}\end{vmatrix} \\
&= x_{B}y_{C}-y_{B}x_{C}-x_{A}y_{C}+x_{A}y_{B}+y_{A}x_{C}-y_{A}x_{B} \\
&= (x_{B}y_{C}+x_{A}y_{B}+y_{A}x_{C})-(y_{A}x_{B}+y_{B}x_{C}+x_{A}y_{C}).

If the determinant is negative, then the polygon is oriented clockwise.  If the determinant is positive, the polygon is oriented counterclockwise.  The determinant  is non-zero if points A, B, and C are non-[[collinear]].  In the above example, with points ordered A, B, C, etc., the determinant is negative, and therefore the polygon is clockwise.


double IsTriangleCounterclockwise(double xa, double ya, double xb, double yb, double xc, double yc)
{return  ((xb*yc + xa*yb +ya*xc) - (ya*xb +yb*xc + xa*yc)); }

int DescribeTriangle(double xa, double ya, double xb, double yb, double xc, double yc)
 double t = IsTriangleCounterclockwise( xa,  ya, xb,  yb,  xc,  yc);
 double a = GiveTriangleArea( xa,  ya, xb,  yb,  xc,  yc);
 if (t>0)  printf("this triangle is oriented counterclockwise,     determinent = %f ; area = %f\n", t,a);
 if (t<0)  printf("this triangle is oriented clockwise,            determinent = %f; area = %f\n", t,a);
 if (t==0) printf("this triangle is degenerate: colinear or identical points, determinent = %f; area = %f\n", t,a);

 return 0;

int main()
 // clockwise oriented triangles 
 DescribeTriangle(-94,   0,  92,  68, 400, 180); //
 DescribeTriangle(4.0, 1.0, 0.0, 9.0, 8.0, 3.0); // clockwise orientation
 //  counterclockwise oriented triangles
 DescribeTriangle(-50.00, 0.00, 50.00,  0.00, 0.00,  0.02); // a "cap" triangle. This example has an area of 1.
 DescribeTriangle(0.0,  0.0, 3.0,  0.0, 0.0,  4.0); // a right triangle. This example has an area of (?? 3 ??)  =  6
 DescribeTriangle(4.0, 1.0, 8.0, 3.0, 0.0, 9.0);  //
 DescribeTriangle(-0.5, 0.0,  0.5,  0.0, 0.0,  0.866025403784439); // an equilateral triangle. This triangle has an area of sqrt(3)/4.

 // degenerate triangles 
 DescribeTriangle(1.0, 0.0, 2.0, 2.0, 3.0, 4.0); // This triangle is degenerate: 3 colinear points.
 DescribeTriangle(4.0, 1.0, 0.0, 9.0, 4.0, 1.0); //2 identical points 
 DescribeTriangle(2.0, 3.0, 2.0, 3.0, 2.0, 3.0); // 3 identical points

 return 0; 

Testing if point is inside polygonEdit

Point_in_polygon_problem : image and source code

gcc p.c -Wall

----------- git --------------------
cd existing_folder
git init
git remote add origin
git add .
git commit
git push -u origin master


#include <stdio.h>

#define LENGTH 6


Argument	Meaning
nvert	Number of vertices in the polygon. Whether to repeat the first vertex at the end is discussed below.
vertx, verty	Arrays containing the x- and y-coordinates of the polygon's vertices.
testx, testy	X- and y-coordinate of the test point.
PNPOLY - Point Inclusion in Polygon Test
W. Randolph Franklin (WRF)


I run a semi-infinite ray horizontally (increasing x, fixed y) out from the test point, 
and count how many edges it crosses. 
At each crossing, the ray switches between inside and outside. 
This is called the Jordan curve theorem.
The case of the ray going thru a vertex is handled correctly via a careful selection of inequalities. 
Don't mess with this code unless you're familiar with the idea of Simulation of Simplicity. 
This pretends to shift the ray infinitesimally down so that it either clearly intersects, or clearly doesn't touch. 
Since this is merely a conceptual, infinitesimal, shift, it never creates an intersection that didn't exist before, 
and never destroys an intersection that clearly existed before.

The ray is tested against each edge thus:

Is the point in the half-plane to the left of the extended edge? and
Is the point's Y coordinate within the edge's Y-range?
Handling endpoints here is tricky.

I run a semi-infinite ray horizontally (increasing x, fixed y) out from the test point, 
and count how many edges it crosses. At each crossing, 
the ray switches between inside and outside. This is called the Jordan curve theorem.
The variable c is switching from 0 to 1 and 1 to 0 each time the horizontal ray crosses any edge. 
So basically it's keeping track of whether the number of edges crossed are even or odd. 
0 means even and 1 means odd.


int pnpoly(int nvert, double *vertx, double *verty, double testx, double testy)
  int i, j, c = 0;
  for (i = 0, j = nvert-1; i < nvert; j = i++) {
    if ( ((verty[i]>testy) != (verty[j]>testy)) &&
	 (testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
       c = !c;

  return c;

void CheckPoint(int nvert, double *vertx, double *verty, double testx, double testy){

int flag;

flag =  pnpoly(nvert, vertx, verty, testx, testy);

   case 0  : printf("outside\n"); break;
   case 1  : printf("inside\n"); break;
   default : printf(" ??? \n");

int main (){

// values from
// number from 0 to (LENGTH-1)
double zzx[LENGTH] = { 13.5,  6.0, 13.5, 42.5, 39.5, 42.5};
double zzy[LENGTH] = {100.0, 70.5, 41.5, 56.5, 69.5, 84.5};

CheckPoint(LENGTH, zzx, zzy, zzx[4]-0.001, zzy[4]);
CheckPoint(LENGTH, zzx, zzy, zzx[4]+0.001, zzy[4]);

return 0;



  • closed/open
  • with/without multiple points


  • tracing /drawing : Generating Discrete Curves
  • sketching
  • sampling
  • Pathfinding or pathing is the plotting, by a computer application, of the shortest route between two points
  • clipping
  • Approximation of digitized curves (with cubic Bézier splines )
  • curve fitting[29]
  • Mending broken lines [30]
  • edge detection


curve simplifyingEdit

reduce the number of points, but still keep the overall shape of the curve = PolyLine Reduction

thick line drawingEdit

thick line

curve drawingEdit

neighborhood variants of 2D algorithms:

  • 8-way stepping (8WS) for 8-direction neighbors of a pixel p(x, y)
  • 4-way stepping (4WS) 4-direction neighbors of a pixel p(x, y)

curve samplingEdit

  • uniform = gives equidistant points
  • adaptive. " an adaptive method for sampling the domain with respect to local curvature. Samples concentration is in proportion to this curvature, resulting in a more efficient approximation—in the limit, a flat curve is approximated by merely two endpoints." [38]

Field linesEdit

Field line[39]


Tracing a curve on a triangular grid

Tracing curve[40]


  • general, (analytic or systematic) = curve sketching[41]
  • local method

trace an external ray for angle t in turns, which means ( description by Claude Heiland-Allen)[42]

  • starting at a large circle (e.g. r = 65536, x = r cos(2 pi t), y = r sin(2 pi t))
  • following the escape lines (like the edges of binary decomposition colouring with large escape radius) in towards the Mandelbrot set.

this algorithm is O(period^2), which means doubling the period takes 4x as long, so it's only really feasible for periods up to a few 1000.

Three Curve Tracing Models[43]

  • Pixel-by-Pixel tracing
  • The bipartite receptive field operator
  • The zoom lens operator




Curve rasterisationEdit


Ray can be parametrised with radius (r)

Closed curveEdit

Simple closed curve ("a connected curve that does not cross itself and ends at the same point where it begins"[45] = having no endpoints) can be parametrized with angle (t).

Border, boundary, contour, edgeEdit

contour modelsEdit

  • snakes = active contour models[46]

border tracingEdit

Filling contourEdit

Edge detectionEdit

Sobel filterEdit

Short introductionEdit

Sobel filter G consist of 2 kernels (masks):

  • Gh for horizontal changes.
  • Gv for vertical changes.
Sobel kernelsEdit
8-point neighborhood on a 2D grid
2D Convolution Animation

The Sobel kernel contains weights for each pixel from the 8-point neighbourhood of a tested pixel. These are 3x3 kernels.

There are 2 Sobel kernels, one for computing horizontal changes and other for computing vertical changes. Notice that a large horizontal change may indicate a vertical border, and a large vertical change may indicate a horizontal border. The x-coordinate is here defined as increasing in the "right"-direction, and the y-coordinate is defined as increasing in the "down"-direction.

The Sobel kernel for computing horizontal changes is:


The Sobel kernel for computing vertical changes is:


Note that:

  • sum of weights of kernels are zero



  • One kernel is simply the other rotated by 90 degrees.[52]
  • 3 weights in each kernal are zero.
Pixel kernelEdit

Pixel kernel A containing central pixel   with its 3x3 neighbourhood:


Other notations for pixel kernel:



unsigned char ul, // upper left
unsigned char um, // upper middle
unsigned char ur, // upper right
unsigned char ml, // middle left
unsigned char mm, // middle = central pixel
unsigned char mr, // middle right
unsigned char ll, // lower left
unsigned char lm, // lower middle
unsigned char lr, // lower right
Pixel 3x3 neighbourhood (with Y axis directed down)

In array notation it is:[54]


In geographic notation usede in cellular aotomats[check spelling] it is central pixel of Moore neighbourhood.

So central (tested) pixel is:


Sobel filtersEdit

Compute Sobel filters (where   here denotes the 2-dimensional convolution operation not matrix multiplication). It is a sum of products of pixel and its weights:


Because 3 weights in each kernal are zero so there are only 6 products.[55]

short Gh = ur + 2*mr + lr - ul - 2*ml - ll;
short Gv = ul + 2*um + ur - ll - 2*lm - lr;

Result is computed (magnitude of gradient):


It is a color of tested pixel.

One can also approximate result by sum of 2 magnitudes:


which is much faster to compute.[56]

  • choose pixel and its 3x3 neighberhood A
  • compute Sobel filter for horizontal Gh and vertical lines Gv
  • compute Sobel filter G
  • compute color of pixel
Sobel filters (2 filters 3x3): image and full c code
Skipped pixel - some points from its neighbourhood are out of the image

Lets take array of 8-bit colors (image) called data. To find borders in this image simply do:

     Gv= - data[iY-1][iX-1] - 2*data[iY-1][iX] - data[iY-1][iX+1] + data[iY+1][iX-1] + 2*data[iY+1][iX] + data[iY+1][iX+1];
     Gh= - data[iY+1][iX-1] + data[iY-1][iX+1] - 2*data[iY][iX-1] + 2*data[iY][iX+1] - data[iY-1][iX-1] + data[iY+1][iX+1];
     G = sqrt(Gh*Gh + Gv*Gv);
     if (G==0) {edge[iY][iX]=255;} /* background */
         else {edge[iY][iX]=0;}  /* boundary */

Note that here points on borders of array (iY= 0, iY = iYmax, iX=0, iX=iXmax) are skipped

Result is saved to another array called edge (with the same size).

One can save edge array to file showing only borders, or merge 2 arrays:

    for(iX=1;iX<iXmax-1;++iX){ if (edge[iY][iX]==0) data[iY][iX]=0;}}

to have new image with marked borders.

Above example is for 8-bit or indexed color. For higher bit colors "the formula is applied to all three color channels separately" (from RoboRealm doc).

Other implementations:

Bad edge position seen in the middle of image. Lines are not meeting in good points, like z = 0

Edge position:

In ImageMagick as "you can see, the edge is added only to areas with a color gradient that is more than 50% white! I don't know if this is a bug or intentional, but it means that the edge in the above is located almost completely in the white parts of the original mask image. This fact can be extremely important when making use of the results of the "-edge" operator."[57]

The result is:

  • doubling edges; "if you are edge detecting an image containing an black outline, the "-edge" operator will 'twin' the black lines, producing a weird result."[58]
  • lines are not meeting in good points.

See also new operators from 6 version of ImageMagick: EdgeIn and EdgeOut from Morphology[59]

Edge thickeningEdit


convert $tmp0 -convolve "1,1,1,1,1,1,1,1,1" -threshold 0 $outfile

SDF Signed Distance FunctionEdit

test external tangency of 2 circlesEdit

 distance between 2 points 
  z1 = x1 + y1*I
  z2 = x2 + y2*I
double GiveDistance(int x1, int y1, int x2, int y2){
  return sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));

mutually and externally tangent circles
Two circles are mutually and externally tangent if distance between their centers is equal to the sum of their radii


double TestTangency(int x1, int y1, int r1, int x2, int y2, int r2){
  double distance;

  distance = GiveDistance(x1, y1, x2, y2);
  return ( distance - (r1+r2));
// return should be zero


  1. Michael Abrash's Graphics Programming Black Book Special Edition
  2. geometrictools Documentation
  3. list-of-algorithms by Ali Abbas
  4. Mitch Richling: 3D Mandelbrot Sets
  5. fractalforums : antialiasing-fractals-how-best-to-do-it/
  6. Matlab : images/morphological-filtering
  7. Tim Warburton : morphology in matlab
  8. A Cheritat wiki : see image showing gamma-correct downscale of dense part of Mandelbropt set
  9. Image processing by Rein van den Boomgaard.
  10. Seeram E, Seeram D. Image Postprocessing in Digital Radiology-A Primer for Technologists. J Med Imaging Radiat Sci. 2008 Mar;39(1):23-41. doi: 10.1016/j.jmir.2008.01.004. Epub 2008 Mar 22. PMID: 31051771.
  11. wikipedi : dense_set
  12. mathoverflow question : is-there-an-almost-dense-set-of-quadratic-polynomials-which-is-not-in-the-inte/254533#254533
  13. fractalforums : dense-image
  14. : m andelbrot-set-various-structures
  15. : technical-challenge-discussion-the-lichtenberg-figure
  16. A Cheritat wiki : see image showing gamma-correct downscale of dense part of Mandelbropt set
  17. fractal forums : pathfinding-in-the-mandelbrot-set/
  18. serious_statistics_aliasing by Guest_Jim
  19. : newton-raphson-zooming
  20. fractalforums : gerrit-images
  21. 5-ways-to-find-the-shortest-path-in-a-graph by Johannes Baum
  22. d3-polygon - Geometric operations for two-dimensional polygons from
  23. github repo for d3-polygon
  24. accurate-point-in-triangle-test by Cedric Jules
  25. stackoverflow question  : how-to-determine-a-point-in-a-2d-triangle
  26. js code
  27. stackoverflow questions : How can I determine whether a 2D Point is within a Polygon?
  28. Punkt wewnątrz wielokąta - W Muła
  29. Matlab examples : curvefitting
  30. Mending broken lines by Alan Gibson.
  31. The Beauty of Bresenham's Algorithm by Alois Zingl
  32. bresenhams-drawing-algorithms
  33. Bresenham’s Line Drawing Algorithm by Peter Occil
  34. Peter Occil
  35. Algorytm Bresenhama by Wojciech Muła
  36. wolfram : NumberTheoreticConstructionOfDigitalCircles
  37. Drawing Thick Lines
  38. IV.4 - Adaptive Sampling of Parametric Curves by Luiz Henrique deFigueiredo
  39. wikipedia: Field line
  40. Curve sketching in wikipedia
  42. help-with-ideas-for-fractal-art
  43. Predicting the shape of distance functions in curve tracing: Evidence for a zoom lens operator by PETER A. McCORMICK and PIERRE JOLICOEUR
  44. stackoverflow question: line-tracking-with-matlab
  45. mathwords: simple_closed_curve
  46. INSA : active-contour-models
  47. Tracing Boundaries in 2D Images by V. Kovalevsky
  48. Calculating contour curves for 2D scalar fields in Julia
  49. Smooth Contours from
  50. d3-contour from
  51. matrixlab - line-detection
  52. Sobel Edge Detector by R. Fisher, S. Perkins, A. Walker and E. Wolfart.
  53. NVIDIA Forums, CUDA GPU Computing discussion by kr1_karin
  54. Sobel Edge by RoboRealm
  55. nvidia forum: Sobel Filter Don't understand a little thing in the SDK example
  56. Sobel Edge Detector by R. Fisher, S. Perkins, A. Walker and E. Wolfart.
  57. ImageMagick doc
  58. Edge operator from ImageMagick docs
  59. ImageMagick doc: morphology / EdgeIn
  60. dilation at HIPR2 by Robert Fisher, Simon Perkins, Ashley Walker, Erik Wolfart
  61. ImageMagick doc: morphology, dilate
  62. Fred's ImageMagick Scripts