essai du trolldi sur le voxel
This commit is contained in:
tth
@@ -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)
|
||||
|
||||
|
||||
Reference in New Issue
Block a user