diff --git a/Fraktalism/common.sh b/Fraktalism/common.sh index 224ccb3..5ce2996 100644 --- a/Fraktalism/common.sh +++ b/Fraktalism/common.sh @@ -1,4 +1,4 @@ -POVOPT=" -d +q9 -a +W1920 +H1080 -v +WT4" +POVOPT=" -d +q9 +a +W1920 +H1080 -v +WT4" diff --git a/Fraktalism/fraktals.f90 b/Fraktalism/fraktals.f90 index 76ff7c6..6707e7c 100644 --- a/Fraktalism/fraktals.f90 +++ b/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 diff --git a/Fraktalism/mkvoxvidz.sh b/Fraktalism/mkvoxvidz.sh index b263c0f..e4c2740 100755 --- a/Fraktalism/mkvoxvidz.sh +++ b/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' - numbers=$(tail -1 WS/camvox.log | \ - awk '{printf " %6.3f %6.3f %6.3f", $1, $2, $3}') - echo "$numbers" + titre='Voxelisation - tTh - Avril 2022' + numbers=$(tail -1 WS/camvox.log | \ + awk '{printf " K=%5d : %6.3f %6.3f %6.3f", \ + $1, $2, $3, $4}') + + echo "numbers " "$numbers" " txtidx " $txtidx convert ${TMPNG} \ - -fill Gray70 \ - -kerning 1 \ - -pointsize 24 \ + -fill Orange \ + -kerning 2 \ + -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 diff --git a/Fraktalism/showvoxels.pov b/Fraktalism/showvoxels.pov index d7e0f82..ee11d14 100644 --- a/Fraktalism/showvoxels.pov +++ b/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>, } } @@ -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 + /* + * un peu de calcul empirique ? + */ + #local TRK = DIMC/2.0000000; + translate <-TRK, -TRK, -TRK> + // rotate } //---------------------------------------------------------------- -// light_source { <-12, 45, -25> color Gray70 } -light_source { <-52, 5, -38> color Yellow*0.45 } -light_source { < 59, 5, 48> color Gray20 } -light_source { < 3, 59, 8> color Red*0.65 } -light_source { < 8, -61, 3> color Green*0.75 } +#declare TriAxe = object +{ +#local Sz = 300; +#local Ra = 0.20; +union { + cylinder { <-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); -#declare XCAM = ECAM * sin(radians(clock)); -#declare YCAM = -12; -#declare ZCAM = ECAM * cos(radians(clock)); -#declare ACAM = 65 + (48*NormClock); +plane { + <0, 1, 0>, -100 + texture { + pigment { color srgb <0.203, 0.185, 0.191> } + 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 { diff --git a/Fraktalism/voxelize.f90 b/Fraktalism/voxelize.f90 index 96e7625..ed345c0 100644 --- a/Fraktalism/voxelize.f90 +++ b/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 - double precision, dimension(4) :: coefs, KB - double precision :: dval + integer :: nbr_points, maxcube + double precision, dimension(4) :: KA, KB, KM + 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 - coefs(3) = -0.64 ; coefs(4) = -2.42 - - KB(1) = 1.51 ; KB(2) = -1.89 - KB(3) = 1.69 ; KB(4) = 0.79 - call compute_pickover(points, KB) - + KA(1) = -1.3402 ; KA(2) = 1.5245 + KA(3) = 1.0966 ; KA(4) = -2.3423 + KB(1) = -1.2100 ; KB(2) = 1.3685 + KB(3) = 1.1237 ; KB(4) = -2.1992 + call interp4dp(KA, KB, KM, delta) + 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)) - write(0, *) "--- maximum of the cube= ", dval + maxcube = MAXVAL(cube) + dmaxcube = DBLE(maxcube) + write(0, *) "--- maxval(cube) = ", maxcube - 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 + call spit_cube_as_union(filename, cube, & + maxcube/1000, dble(9000.00)) - 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) - character(*) :: fname - type(integer), dimension(:,:,:), intent(in) :: datas - integer :: notused +subroutine print_cube(cube, points, scaling) + type(integer), dimension(:,:,:), intent(in) :: cube + type(t_point3d), dimension(:), intent(in) :: points + double precision, intent(in) :: scaling - integer :: fd, errcode, foo + integer :: foo - open (newunit=fd, file='WS/k-pick.txt', & - status='unknown', position='append', & + 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)) / 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)