In a further exploration of using PovRay to do fast calculation of Voronoi polygons, let’s look to a real stream system as an example.  Here’s where the magic comes out, and where medial axes are found.

Here’s the povray code:

#include "transforms.inc"
#version 3.6;

#include "edge_coords1.inc"

#declare Camera_Location = <2258076, 10, 659923>;
#declare Camera_Lookat   = <2258076, 0, 659923>;
 #declare Camera_Size = 1800;

camera {
orthographic
location Camera_Location
look_at Camera_Lookat
right Camera_Size*x
up Camera_Size*y
}

background {color <0.5, 0.5, 0.5>}

//union {

#declare Rnd_1 = seed (1153);

#declare LastIndex = dimension_size(edge_coords, 1)-2;
#declare Index = 0;
#while(Index <= LastIndex)

cone {
<0.0, -5, 0.0>, 30,  <0.0, 5, 0.0>, 0.0
     hollow
material {
texture {
pigment {
rgb <rand(Rnd_1)*2 + 1, rand(Rnd_1)*5 + 2, rand(Rnd_1)*5 +5>
}
finish {
diffuse 0.6
                 brilliance 1.0
}
}
interior {
ior 1.3
}
}
translate edge_coords[Index]
}


#declare Index = Index + 1;

#end

#declare edge_coords = array[16842]
{<2256808.71000000000, 0.00000000000, 660094.46988100000> ,
<2256810.06170000000, 0.00000000000, 660092.99580200000> ,
<2256811.41308000000, 0.00000000000, 660091.52172400000> ,
  <2256811.68965000000, 0.00000000000, 660091.21988700000> ,
<2256813.46393000000, 0.00000000000, 660090.29698900000> ,
<2256815.23820000000, 0.00000000000, 660089.37376200000> ,
<2256817.01248000000, 0.00000000000, 660088.45086400000> ,
  <2256818.78675000000, 0.00000000000, 660087.52763700000> ,
<2256820.56103000000, 0.00000000000, 660086.60473900000> ,
<2256822.33530000000, 0.00000000000, 660085.68151200000> ,
<2256824.10958000000, 0.00000000000, 660084.75861400000> ,
  <2256825.88352000000, 0.00000000000, 660083.83538800000> ,
<2256827.65780000000, 0.00000000000, 660082.91248900000> ,
<2256829.43207000000, 0.00000000000, 660081.98926300000> ,
<2256831.20635000000, 0.00000000000, 660081.06636400000> ,