essai du trolldi sur le voxel

Fraktalism/.gitignore

@@ -5,9 +5,11 @@ lorentz
 voxelize
 evolvopick
 henon
+essai
 
 WS/*.dat
 WS/*.txt
+WS/*.inc
 
 *.pgm
 *.gif

Fraktalism/Makefile

@@ -1,7 +1,7 @@
 
-all:	voxelize evolvopick pickover julia lorentz
+all:	voxelize evolvopick pickover julia lorentz essai
 
-GFOPT  = -Wall -Wextra -time -O -Imods/
+GFOPT  = -Wall -Wextra -time -g -O -Imods/
 
 # ---------------------------------------------
 
@@ -18,6 +18,9 @@ OBJS  = mods/spitpgm.o mods/points3d.o fraktals.o
 
 # ---------------------------------------------
 
+essai:		essai.f90 Makefile $(OBJS)
+	gfortran $(GFOPT) $<  $(OBJS)    -o $@
+
 henon:		henon.f90 Makefile $(OBJS)
 	gfortran $(GFOPT) $<  $(OBJS)    -o $@
 

Fraktalism/common.sh

@@ -1,4 +1,4 @@
 
 
-POVOPT=" -d +q9 -a +W1920 +H1080 -v +WT4"
+POVOPT=" -d +q9 +a +W1920 +H1080 -v +WT4"
 

Fraktalism/essai.f90

@@ -0,0 +1,9 @@
+!-----------------------------------------------------
+program essai
+  use spitpgm
+
+  implicit none
+
+
+!-----------------------------------------------------
+end program

Fraktalism/fraktals.f90

@@ -151,7 +151,7 @@ end subroutine lorentz_0
 !-----------------------------------------------------------
 !           -- some support functions --
 !-----------------------------------------------------------
-!       usage : evolvopick
+!       usage : evolvopick & voxelize
   subroutine interp4dp (ina, inb, out, dpk)
     double precision, dimension(4), intent(in)   :: ina, inb
     double precision, dimension(4), intent(out)  :: out

Fraktalism/mkvoxvidz.sh

@@ -10,20 +10,22 @@ if [ $? -ne 0 ] ; then
 	exit 1
 fi
 
-./voxelize > WS/voxels.dat
-./vox2inc.awk < WS/voxels.dat > WS/voxels.inc
-grep NbrVox WS/voxels.inc
 
 TMPNG="/dev/shm/voxvidz.png"
-NBIMG=360
+NBIMG=1500
+
+printf "#declare NbImg = %d;\n" $NBIMG | tee WS/nbimg.inc
 
 for idx in $( seq   0 $(( NBIMG - 1)) )
 do
 
 	dst=$(printf "frames/voxel/%05d.png" $idx)
-	delta=$( echo "scale=6 ; ${idx}/360" | bc -l)
+	delta=$( echo "scale=6 ; ${idx}/${NBIMG}" | bc -l)
 	echo "Renderbox work on "$dst"     delta = "  $delta
 
+	./voxelize "WS/voxels.inc" $delta
+	grep 'NbrVox' "WS/voxels.inc"
+
 	povray		-ishowvoxels.pov -K$idx ${POVOPT}	\
 			-O${TMPNG} 2> WS/toto.err
 	if [ $? -ne 0 ] ; then
@@ -33,24 +35,30 @@ do
 		#exit 1
 	fi
 
-	titre='Voxelisation - tTh - avril 2022'
+	titre='Voxelisation - tTh - Avril 2022'
-	numbers=$(tail -1 WS/camvox.log |		\
+	numbers=$(tail -1 WS/camvox.log |			\
-		  awk '{printf "  %6.3f   %6.3f   %6.3f", $1, $2, $3}')
+		  awk '{printf "  K=%5d : %6.3f   %6.3f   %6.3f",	\
-	echo "$numbers"
+					$1, $2, $3, $4}')
+
+	echo "numbers " "$numbers" "     txtidx " $txtidx
 
 	convert		${TMPNG}			\
-			-fill  Gray70			\
+			-fill  Orange			\
-			-kerning 1			\
+			-kerning 2			\
-			-pointsize 24			\
+			-pointsize 32			\
 			-font AvantGarde-Book		\
 			-gravity South-West		\
-			-annotate +20+32 "$titre"	\
+			-annotate +20+45 "$titre"	\
+			-pointsize 24			\
 			-annotate +20+10 "$numbers"	\
+							\
 			$dst
 
 	grep 'Parse Time' WS/toto.err
 	grep 'Trace Time' WS/toto.err
 
+	echo ; sleep 15
+
 done
 
-./encode.sh frames/voxel/ voxel-2.mp4
+./encode.sh frames/voxel/ voxel-3.mp4

Fraktalism/showvoxels.pov

@@ -1,7 +1,7 @@
 /*
  *		SHOW VOXELS
  *
- *  see also : vox2inc.awk voxelize.f90
+ *  see also : vox2inc.awk and voxelize.f90
  */
 
 #version 3.7;
@@ -17,7 +17,7 @@ global_settings {
 #declare VOXEL = object
 {
 // sphere { 0, 1.18 }
-#local D = 1.0111;
+#local D = 1.92;
 box { <-D, -D, -D>, <D, D, D> }
 }
 
@@ -26,33 +26,66 @@ object {
 	Voxels
 	texture {
 		pigment { color White }
-		finish { phong 0.8 specular 0.8}
+		finish { phong 0.6 specular 0.8 }
 		}
-	translate <-Bary_X, -Bary_Y, -Bary_Z>
+	/*
-	rotate <clock*2, 0, clock*1.414>
+	 * un peu de calcul empirique ?
+	 */
+	#local TRK = DIMC/2.0000000;
+	translate <-TRK, -TRK, -TRK>
+	// rotate <clock*0.22, 0, clock*0.17>
 	}
 
 //----------------------------------------------------------------
-// light_source { <-12,  45, -25> color Gray70 }
+#declare TriAxe = object
-light_source { <-52,  5, -38>   color Yellow*0.45 	}
+{
-light_source { < 59,  5,  48>   color Gray20		}
+#local Sz = 300;
-light_source { <  3,  59,  8>   color Red*0.65 		}
+#local Ra = 0.20;
-light_source { <  8,  -61,  3>   color Green*0.75 	}
+union	{
+	cylinder { <-Sz, 0, 0>, <Sz, 0, 0>, Ra pigment { color Red }    }
+	cylinder { <0, -Sz, 0>, <0, Sz, 0>, Ra pigment { color Green }  }
+	cylinder { <0, 0, -Sz>, <0, 0, Sz>, Ra pigment { color Blue }   }
+	}
+finish { phong 0.6 specular 0.8 }
+}
 
-#declare NormClock = (clock/360);
+object { TriAxe }
 
-#declare ECAM  =  88 - (73*NormClock);
+plane   {
-#declare XCAM  =  ECAM * sin(radians(clock));
+        <0, 1, 0>, -100
-#declare YCAM  =  -12;
+        texture {
-#declare ZCAM  =  ECAM * cos(radians(clock));
+                pigment { color srgb <0.203, 0.185, 0.191> } 
-#declare ACAM  =  65 + (48*NormClock);
+                finish { phong 0.18 metallic 0.55 }
+                }
+        }
+
+light_source { <-29,  45, -27> color Gray70 }
+light_source { <-52,  5, -48>   color Yellow*0.45 	}
+light_source { < 59,  45,  48>   color Gray20		}
+light_source { < 59, -45,  48>   color Gray20		}
+light_source { <  9,  59,  18>   color Red*0.65 		}
+light_source { <  8,  -48,  3>   color Green*0.75 	}
+
+#include "WS/nbimg.inc"
+#declare NormClock = (clock/NbImg);
+
+#declare ECAM    =  190 - (100*NormClock);
+#declare CKsmall =  NormClock * 77.20;
+#declare Offset  =  0.10;
+#declare XCAM  =  ECAM * (sin(radians(CKsmall)) + Offset);
+#declare YCAM  =  16;
+#declare ZCAM  =  ECAM * (cos(radians(CKsmall)) + Offset);
+#declare ACAM  =  65 + (53*NormClock);
+
+// #declare XCAM = ECAM * 0.8;
+// #declare ZCAM = ECAM * 0.35;
 
 #if (0 = clock)
   #fopen CL "WS/camvox.log" write
 #else 
   #fopen CL "WS/camvox.log" append
 #end
-#write (CL, NormClock, " ", ECAM, " ", ACAM, "\n")
+#write (CL, clock, " ", NormClock, " ", ECAM, " ", ACAM, "\n")
 #fclose CL
 
 camera  {

Fraktalism/voxelize.f90

@@ -9,37 +9,47 @@ program voxelize
 
   implicit none
 
-  integer, parameter                          :: DIM = 100
+  integer, parameter                          :: DIMC = 200
   integer, dimension(:,:,:), allocatable      :: cube
   type(t_point3d), dimension(:), allocatable  :: points
-  integer                                     :: errcode, foo
+  integer                                     :: errcode, foo, argc
   integer                                     :: ix, iy, iz
-  integer                                     :: nbr_points
+  integer                                     :: nbr_points, maxcube
-  double precision, dimension(4)              :: coefs, KB
+  double precision, dimension(4)              :: KA, KB, KM
-  double precision                            :: dval
+  double precision                            :: dmaxcube, delta
+  character(200)                              :: filename, string
 
   write(0, *) "--- start of voxelize"
 
-  allocate (cube(DIM,DIM,DIM), stat=errcode)
+  argc = IARGC()
+  if (2 .NE. argc) then
+    STOP ": VOXELIZE NEED PARAMETERS !"
+  endif
+
+  call getarg(1, filename)
+  call getarg(2, string)       ;   read (string, *) delta
+  write(0, "( ' --- delta = ', F11.6)") delta
+
+  allocate (cube(DIMC,DIMC,DIMC), stat=errcode)
   if (0 .NE. errcode) then
     STOP " : NO ENOUGH MEMORY FOR CUBE"
   endif
 
-  nbr_points = 2800000
+  nbr_points = 7000000
   allocate(points(nbr_points), stat=errcode)
   if (0 .NE. errcode) then
     STOP " : NO ENOUGH MEMORY FOR POINTS"
   endif
 
-  coefs(1) =  2.23     ;     coefs(2) =  0.42
+  KA(1) =  -1.3402     ;     KA(2) =   1.5245
-  coefs(3) = -0.64     ;     coefs(4) = -2.42
+  KA(3) =   1.0966     ;     KA(4) =  -2.3423
-
+  KB(1) =  -1.2100     ;     KB(2) =   1.3685
-  KB(1) =  1.51        ;     KB(2) = -1.89
+  KB(3) =   1.1237     ;     KB(4) =  -2.1992
-  KB(3) =  1.69        ;     KB(4) =  0.79
+  call interp4dp(KA, KB, KM, delta)
-  call compute_pickover(points, KB)
+  write(0, "(' --- coefs = ', 4F11.6)") KM
-
+  call compute_pickover(points, KM)
   call clear_cube(cube)
-
+  ! 
   ! and now, we loop over all the pre-computed
   ! points of the attractor
   ! 
@@ -50,19 +60,13 @@ program voxelize
     cube(ix,iy,iz) = cube(ix,iy,iz) + 1
   enddo
 
-  dval = DBLE(MAXVAL(cube))
+  maxcube = MAXVAL(cube)
-  write(0, *) "--- maximum of the cube= ", dval
+  dmaxcube = DBLE(maxcube)
+  write(0, *) "--- maxval(cube) = ", maxcube
 
-  do foo=1, nbr_points
+  call spit_cube_as_union(filename,   cube,       &
-    call fcoor2icoor(points(foo)%x, ix)
+                          maxcube/1000, dble(9000.00))
-    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
 
-  dval = DBLE(MAXVAL(cube))
   write(0, *) "--- end of voxelize"
 
 !-----------------------------------------------------
@@ -76,15 +80,17 @@ subroutine fcoor2icoor(in, out)
   integer                            :: outvalue
 
   invalue = (in + 2.0) / 2.0
-  outvalue = int(invalue * real(DIM/2))
+  outvalue = int(invalue * real(DIMC/2))
 
   ! add molly-guard here
   out = outvalue
   if (outvalue .LT. 1)      out = 1
-  if (outvalue .GE. DIM)    out = DIM-1
+  if (outvalue .GE. DIMC)    out = DIMC-1
 
 end subroutine
 !------------------------------------------------------------
+!           USELESS USE OF LOOPS !
+
 subroutine clear_cube(cube)
   type(integer), dimension(:,:,:), intent(out)  :: cube
 
@@ -100,20 +106,93 @@ subroutine clear_cube(cube)
 
 end subroutine
 !------------------------------------------------------------
-subroutine spit_cube_as_union(fname, datas, notused)
+subroutine print_cube(cube, points, scaling)
-  character(*)                                 :: fname
+  type(integer), dimension(:,:,:), intent(in)  :: cube
-  type(integer), dimension(:,:,:), intent(in)  :: datas
+  type(t_point3d), dimension(:), intent(in)    :: points
-  integer                                      :: notused
+  double precision, intent(in)                 :: scaling
 
-  integer             :: fd, errcode, foo
+  integer                     :: foo
 
-  open (newunit=fd, file='WS/k-pick.txt',            &
+  do foo=1, nbr_points
-        status='unknown', position='append',         &
+    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)) / scaling
+  enddo
+
+end subroutine
+
+!------------------------------------------------------------
+!           generation Povray STL source file !
+
+subroutine spit_cube_as_union(fname, voxels, limit, scaling)
+  character(*), intent(in)                     :: fname
+  type(integer), dimension(:,:,:), intent(in)  :: voxels
+  integer, intent(in)                          :: limit
+  double precision, intent(in)                 :: scaling
+
+  integer             :: fd, errcode
+  integer             :: ix, iy, iz, maxv
+  integer             :: nbrvox = 0
+  double precision    :: bx, by, bz, valeur
+  character(200)      :: chaine
+
+  !   molly-guard
+  maxv = limit
+  if (maxv .LT. 2)   maxv = 2
+
+  open (newunit=fd, file=trim(fname),           &
+        status='replace',                       &
         action='write', iostat=errcode)
   if (0 .NE. errcode) then
-    STOP ' : SPITUNION : FAIL OPEN OUTPUT FILE'
+    STOP ' : SPIT UNION, FAIL OPEN OUTPUT FILE'
   endif
 
+  write(fd, *) "// generated file, don't touch it bastard !"
+  write(fd, *) "// version 2.09"
+  write(fd, *) "#declare DIMC = ", DIMC, ";"
+  write(fd, *) "#declare Voxels = object {"
+  write(fd, *) "union {"
+
+  bx = 0.0   ;   by = 0.0   ;   bz = 0.0
+
+  do ix=lbound(voxels,1), ubound(voxels,1)
+    do iy=lbound(voxels,2), ubound(voxels,2)
+      do iz=lbound(voxels,3), ubound(voxels,3)
+
+        if (cube(ix,iy,iz) .LT. maxv) then
+          ! print *, "foo = ", foo, cube(ix,iy,iz)
+          cycle                                   ! REDO FROM START
+        endif
+
+        nbrvox  = nbrvox + 1
+        bx = bx + dble(ix)
+        by = by + dble(iy)
+        bz = bz + dble(iz)
+        valeur = DBLE(cube(ix,iy,iz)) / scaling
+        ! XXX
+        if (valeur .GT. 1.0) then
+          valeur = 1.0
+        endif
+
+        write(chaine, "( 'translate <', I4, ',', I4, ',', I4, '> ' )")   &
+                            ix, iy, iz
+        write(unit=fd,                                                   &
+              fmt="( 'object { VOXEL scale ', F11.6, 1X, A, ' }'  )",    &
+              iostat=errcode)                                            &
+              valeur, trim(chaine)
+      enddo
+    enddo
+  enddo
+
+  write(fd, *) "} }"
+  write(fd, *) "// limit = ", limit
+  write(fd, "( '#declare NbrVox = ', I9,    ';' )")   nbrvox
+  write(fd, "( '#declare Bary_X = ', F11.6, ';' )")   bx / dble(nbrvox)
+  write(fd, "( '#declare Bary_Y = ', F11.6, ';' )")   by / dble(nbrvox)
+  write(fd, "( '#declare Bary_Z = ', F11.6, ';' )")   bz / dble(nbrvox)
 
   close(fd)