more work on voxelize

Fraktalism/Makefile

@@ -1,4 +1,6 @@
 
+all:	voxelize evolvopick pickover julia lorentz
+
 GFOPT  = -Wall -Wextra -time -O -Imods/
 
 # ---------------------------------------------
@@ -37,6 +39,6 @@ lorentz.pgm:	lorentz Makefile
 	./lorentz $@ > /dev/null
 
 pickover.pgm:	pickover Makefile
-	time ./pickover $@ > /dev/null
+	./pickover $@ > /dev/null
 
 # ---------------------------------------------

Fraktalism/WS/.gitignore

@@ -1 +1,2 @@
 *.inc
+*.err

Fraktalism/common.sh

@@ -0,0 +1,4 @@
+
+
+POVOPT=" -d +q9 +a +W1024 +H768 -v +WT4"
+

Fraktalism/fraktals.f90

@@ -152,7 +152,6 @@ end subroutine lorentz_0
 !           -- some support functions --
 !-----------------------------------------------------------
 
-
 !-----------------------------------------------------------
 
 function dist0 (x, y)

Fraktalism/mkvoxvidz.sh

@@ -0,0 +1,33 @@
+#!/bin/bash
+
+echo
+source "./common.sh"
+
+make voxelize
+if [ $? -ne 0 ] ; then
+	echo
+	echo "Make error " $?
+	exit 1
+fi
+
+./voxelize > WS/voxels.dat
+./vox2inc.awk < WS/voxels.dat > WS/voxels.inc
+
+for idx in $(seq 0 359)
+do
+
+	dst=$(printf "frames/voxel/%05d.png" $idx)
+	echo "renderbox work on " $dst
+	povray		-ishowvoxels.pov -K$idx ${POVOPT}	\
+			-O${dst} 2> WS/toto.err
+	if [ $? -ne 0 ] ; then
+		echo "ERROR ERROR ERROR ERROR"
+		tail -20 WS/toto.err
+		exit 1
+	fi
+	grep 'Parse Time' WS/toto.err
+	grep 'Trace Time' WS/toto.err
+
+done
+
+./encode.sh frames/voxel/ foo.mp4

Fraktalism/showvoxels.pov

@@ -0,0 +1,49 @@
+/*
+ *		SHOW VOXELS
+ *
+ *  see also : vox2inc.awk voxelize.f90
+ */
+
+#version 3.7;
+
+global_settings {
+                ambient_light rgb <0.12, 0.04, 0.04>
+                assumed_gamma 1.0
+                }
+//----------------------------------------------------------------
+
+#include  "colors.inc"
+
+#declare VOXEL = object
+{
+// sphere { 0, 1.18 }
+#local D = 0.880;
+box { <-D, -D, -D>, <D, D, D> }
+}
+
+#include "WS/voxels.inc"
+object {
+	Voxels  pigment { color White }
+	translate <-Bary_X, -Bary_Y, -Bary_Z>
+	rotate <40, 0, 29>
+	}
+
+//----------------------------------------------------------------
+light_source { <-12, 45, -25> color Gray90 }
+light_source { < 52, -5,  25> color Gray90 }
+
+#declare ECAM = 92;
+#declare XCAM =  ECAM * sin(radians(clock));
+#declare YCAM =  -20;
+#declare ZCAM =  ECAM * cos(radians(clock));
+
+camera  {
+	location   	<XCAM, YCAM, ZCAM>
+	// look_at    	<Bary_X, Bary_Y, Bary_Z>
+	look_at		<0, 0, 0>
+        right           x*image_width/image_height
+        angle           64
+        }
+
+//----------------------------------------------------------------
+//----------------------------------------------------------------

Fraktalism/vox2inc.awk

@@ -0,0 +1,32 @@
+#!/usr/bin/awk -f
+
+BEGIN	{
+	maxcount = 0
+	nbrvox = 0
+	bx = by = bz = 0.0
+	print "// generated file, don't touch it bastard !"
+	print "#declare Voxels = object {"
+	print "union {"
+	}
+
+$4 > 120 {
+	count = $4
+	value = $5
+	if (count > maxcount) 
+		{ maxcount = count }
+	nbrvox++;
+	bx += $1
+	by += $2
+	bz += $3
+	printf "object { VOXEL scale %f translate <%f, %f, %f> }    // %d \n",	\
+					value, $1, $2, $3, count
+	}
+
+END	{
+	print "} } // done, ",     NR, " records read"
+	print "#declare VoxMaxcount = ", maxcount, ";"
+	print "#declare NbrVox   = ",   nbrvox, ";"
+	print "#declare Bary_X   = ", bx/nbrvox, ";"
+	print "#declare Bary_Y   = ", by/nbrvox, ";"
+	print "#declare Bary_Z   = ", bz/nbrvox, ";"
+	}

Fraktalism/voxelize.f90

@@ -1,27 +1,29 @@
 !-----------------------------------------------------
 !           VOXELIZE
 !           ========
-!    this is the main program
+!    this is the main program, see also mkvoxvidz.sh
+!    showvoxels.pov and vox2inc.awk
 !-----------------------------------------------------
 program voxelize
   use fraktals
 
-  integer, parameter                          :: DIM = 500
+  integer, parameter                          :: DIM = 100
   integer, dimension(:,:,:), allocatable      :: cube
   type(t_point3d), dimension(:), allocatable  :: points
-  integer                                     :: errcode, foo
+  integer                                     :: errcode, foo, maxi
   integer                                     :: ix, iy, iz
   double precision, dimension(4)              :: coefs
+  double precision                            :: dval
 
-  foo = (DIM*DIM*DIM) / (1024)
+  ! XXX foo = (DIM*DIM*DIM) / (1024)
-  PRINT *, "memory request for cube (in Kwords) ", foo
+  ! XXX PRINT *, "memory request for cube (in Kwords) ", foo
 
   allocate (cube(DIM,DIM,DIM), stat=errcode)
   if (0 .NE. errcode) then
     STOP " : NO ENOUGH MEMORY FOR CUBE"
   endif
 
-  nbr_points = 99999
+  nbr_points = 3500000
   allocate(points(nbr_points), stat=errcode)
   if (0 .NE. errcode) then
     STOP " : NO ENOUGH MEMORY FOR POINTS"
@@ -33,14 +35,64 @@ program voxelize
 
   call clear_cube(cube)
 
+  ! and now, we loop over all the pre-computed
+  ! points of the attractor
+  ! 
   do foo=1, nbr_points
-    ix = nint(points(foo)%x * dble(DIM))
+    call fcoor2icoor(points(foo)%x, ix)
-    iy = nint(points(foo)%y * dble(DIM))
+    call fcoor2icoor(points(foo)%y, iy)
-    iz = nint(points(foo)%z * dble(DIM))
+    call fcoor2icoor(points(foo)%z, iz)
+    cube(ix,iy,iz) = cube(ix,iy,iz) + 1
   enddo
+
+  dval = DBLE(MAXVAL(cube))
+
+  write(0, *) "--- cube maximum = ", dval
+
+  do foo=1, nbr_points
+    call fcoor2icoor(points(foo)%x, ix)
+    call fcoor2icoor(points(foo)%y, iy)
+    call fcoor2icoor(points(foo)%z, iz)
+    print *, ix, iy, iz,                   &
+             cube(ix,iy,iz),               &
+             DBLE(cube(ix,iy,iz)) / dval
+  enddo
+
 !-----------------------------------------------------
 contains
 !-----------------------------------------------------
+  ! or maybe, we can write a function ?
+subroutine fcoor2icoor(in, out)
+  double precision, intent(in)       :: in
+  integer, intent(out)               :: out
+  double precision                   :: invalue
+  integer                            :: outvalue
+
+  invalue = (in + 2.0) / 2.0
+  outvalue = int(invalue * real(DIM/2))
+
+  ! add molly-guard here
+  out = outvalue
+  if (outvalue .LT. 1)      out = 1
+  if (outvalue .GE. DIM)    out = DIM-1
+
+end subroutine
+!-----------------------------------------------------
+
+subroutine clear_cube(cube)
+  type(integer), dimension(:,:,:), intent(out)  :: cube
+
+  integer        :: i, j, k
+
+  do i=lbound(cube, 1), ubound(cube, 1)
+    do j=lbound(cube, 2), ubound(cube, 2)
+      do k=lbound(cube, 3), ubound(cube, 3)
+        cube(i, j, k) = 0
+      enddo
+    enddo
+  enddo
+
+end subroutine
 !-----------------------------------------------------
 end program voxelize
 !-----------------------------------------------------