\ *****************************************************************
\                   Loading all necessary modules
\ *****************************************************************
\needs float import float also float
\ \needs locals| | include locals.fs
\ \needs $@ | include string.fs
[IFUNDEF] farray.fs include farray.fs [THEN]
[IFUNDEF] string.fs include string.fs [THEN]
\ only forth definitions also assembler
\ warning off
\ \needs locals| | include locals.fs
\ warning on
\ \needs $@ | include string.fs
\ previous also float
\ [IFUNDEF] farray.fs include farray.fs [THEN]
\ *****************************************************************
\                Words necessary to define variables
\ *****************************************************************
forth also
: file->dim ( addr u -- x y z l )
    r/o open-file throw 4 cells allocate throw
    { fid t_pad |
    t_pad 4 cells fid read-file throw drop
    t_pad dup @ swap cell+ dup @ swap cell+ dup @ swap cell+ @
    fid close-file throw
    t_pad free throw } ;
 previous
\ 0 value отладка

\ *****************************************************************
\                          Global variables
\ *****************************************************************
\ ---------------------------User-defined--------------------------
\ A file used by this program should be saven in the following
\ binary format:
\ First 32 bytes (4 32-bit words)
\ x y and z dimensions of the array
\ l - length of an element in bytes
\ The array that follows should be saved as
\ [a[:,1,1]..a[:,n,1]]..[a[:,1,m]..a[:,n,m]]
\ -----------------------------------------------------------------
variable 3dfile
s" brainsmall.dat"
3dfile $!

3dfile $@ file->dim
constant lformat
constant zdim
constant ydim
constant xdim

\ define fetching operation depending on the size of an element
lformat 1 = [IF]
    ' c@ alias fetch@
[THEN]
lformat 2 = [IF]
    ' w@ alias fetch@
[THEN]
lformat 4 = [IF]
    ' @ alias fetch@
[THEN]

\ 3d array that'll hold the final CT cube
xdim ydim zdim lformat * * * allocate throw constant cube{{{

\ find the closest power of 2 to the array size
: nearest2pow ( n -- nearest bigger 2^{c} )
    2 begin 2* 2dup <= until nip ;

xdim ydim max nearest2pow value texdimxy
zdim xdim max nearest2pow value texdimzx
zdim ydim max nearest2pow value texdimzy

\ reserve texture array
texdimxy dup 3 * * allocate throw constant imtext-xy 
texdimzy dup 3 * * allocate throw constant imtext-zy
texdimzx dup 3 * * allocate throw constant imtext-zx

\ calsulate texture offsets
texdimxy xdim - 2/ constant xy-offset-x
texdimxy ydim - 2/ constant xy-offset-y
texdimxy xy-offset-y * xy-offset-x + 3 * value xy-offset

texdimzx xdim - 2/ constant zx-offset-x
texdimzx zdim - 2/ constant zx-offset-z
texdimzx zx-offset-x * zx-offset-z + 3 * value zx-offset

texdimzy ydim - 2/ constant zy-offset-y
texdimzy zdim - 2/ constant zy-offset-z
texdimzy zy-offset-y * zy-offset-z + 3 * value zy-offset

forward gmax
\ *****************************************************************
\                              Words
\ *****************************************************************
\ maximal element of my 2D array
: max-arr ( array xdim ydim -- max )
    0 -rot * 0 do
	swap dup lformat i * + aligned fetch@ rot 2dup
	> if drop
	else nip then
    loop
    nip ;
\ minimal element of my 2D array
: min-arr ( array xdim ydim -- max )
    0 -rot * 0 do
	swap dup lformat i * + aligned fetch@ rot 2dup
	< if drop
	else nip then
    loop
    nip ;

: maxxz ( array xdim zdim -- max )
    \   0 -rot * 0 do
    drop 2drop gmax
    \    loop nip
;
: maxyz ( array ydim zdim -- max )
    \   0 -rot * 0 do
    drop 2drop gmax
    \    loop nip
;

: gulp->array ( array fid x y z l -- )
    { fd x y z l |
    4 cells extend fd reposition-file throw
    x y z l * * * fd read-file throw drop } ;

\ functions to work with palletes are written basen on
\ the source code of http://amide.sourceforge.net/
: temp1 ( datum max min -- f: tmp )
    rot >r dup r> swap - s>f - s>f f/ ;
: red-temp ( datum max min -- r g b )
    >r 2dup r@ temp1 r> 2dup - !255 s>f f/ 
    { datum max min f: temp f: scale |
    temp !1.0 f> if $FF $FF $FF exit then
    temp !0.0 f< if 0 0 0 exit then    
    temp !0.7 f>= if $FF else datum min - s>f scale f* !0.7 f/ f>s then
    temp !0.5 f>= if
	datum s>f min s>f !2 f/ f- max s>f !2 f/ f-
	scale f* !2 f* f>s else 0  then    
    temp !0.5 f>= if
    	datum s>f min s>f !2 f/ f- max s>f !2 f/ f-
	scale f* !2 f* f>s else 0 then } ;
: temp2 ( datum max min -- f: tmp )
    temp1 !500 f*
    fdup !255 f> if !511 fswap f- then
    fdup !0 f< if fdrop !0 then ;
: bwb-temp ( datum max min -- r g b )
    temp2 f>s dup dup ;
: blue-temp ( datum max min -- r g b )
    red-temp swap rot ;
: green-temp ( datum max min -- r g b )
    red-temp rot swap ;
: hot-metal-temp ( datum max min -- r g b )
    >r 2dup r@ temp1 r>
    { datum max min f: temp |
    temp !1.0 f> if $FF $FF $FF exit then
    temp !0.0 f< if 0 0 0 exit then
    temp [ !182 !255 f/ ] fliteral f>= if $FF else temp !255 f* [ !255 !182 f/ ] fliteral f* f>s then
    temp [ !128 !255 f/ ] fliteral f<
    if $00 else
	temp [ !219 !255 f/ ] fliteral f>=
	if
	    $FF
	else
	    temp [ !128 !255 f/ ] fliteral f- !255 f* [ !255 !91 f/ ] fliteral f* f>s
	then
    then
    temp [ !192 !255 f/ ] fliteral f>= if $FF else
	temp [ !192 !255 f/ ] fliteral f- !255 f* [ !255 !63 f/ ] fliteral f* f>s then } ;
: hot-blue-temp ( datum max min -- r g b )
    hot-metal-temp swap rot ;
: old-temp ( datum max min -- r g b )
    - swap s>f !255 f* s>f f/ fround f>s dup 0<> if 1- then dup dup ;    
: swap3 ( a b c -- c b a )
    -rot swap ;

' bwb-temp alias current-temp
\ array to image area in memory
: fillimagexy ( arr{{ xdim ydim graynorm itext -- )
    xy-offset +
    { arr{{ xd yd gn itext |
    yd 0 do
	xd 0 do
	    arr{{ xd yd * 1- i j xd * - - lformat * +
	    fetch@ gn 0 current-temp
	    itext j texdimxy 3 * * + i 3 * +
	    dup >r  c!
	    r@  1+ c!
	    r>  2+ c!
	loop
    loop
    } ;

: xzslice ( i j n -- addr-increment )
    xdim * + swap [ ydim xdim * ] literal * + lformat * ;
: fillimagezx ( n arr{{ xdim zdim graynorm itext -- )
    zx-offset +
    { n arr{{ xd zd gn itext |
    xd 0 do
	zd 0 do
	    arr{{ i j n xzslice +
	    fetch@ gn 0 current-temp
	    itext j texdimzx 3 * * + i 3 * +
	    dup >r  c!
	    r@  1+ c!
	    r>  2+ c!	    
	loop
    loop } ;

: yzslice ( i j n -- addr-increment )
    swap xdim * + swap [ ydim xdim * ] literal * + lformat * ;    
: fillimagezy ( n arr{{ ydim zdim graynorm itext -- )
    zy-offset +
    { n arr{{ yd zd gn itext |
    yd 0 do
	zd 0 do
	    arr{{ i j n yzslice +
	    fetch@ gn 0 current-temp
	    itext j texdimzy 3 * * + i 3 * +
	    dup >r  c!
	    r@  1+ c!
	    r>  2+ c!	    	    
	loop
    loop } ;

: xy->textr ( zn -- )
    dup zdim <= if
	cube{{{ xdim ydim lformat * * rot * +
	dup xdim ydim max-arr xdim ydim rot
	imtext-xy dup texdimxy dup 3 * * erase fillimagexy
    else drop ." index exceeds matrix dimension" cr then ;

: zy->textr ( xn -- )
    dup xdim <= if
	cube{{{ ydim zdim
	cube{{{ ydim zdim maxyz imtext-zy dup texdimzy dup 3 * * erase fillimagezy
    else drop ." index exceeds matrix dimension" cr then ;

: zx->textr ( yn -- )
    dup ydim <= if 
	cube{{{ xdim zdim
	cube{{{ xdim zdim maxxz imtext-zx dup texdimzx dup 3 * * erase fillimagezx
    else drop ." index exceeds matrix dimension" cr then ;

: textr->ppm ( itext xdim ydim -- )
    { im x y |
    ." P3" cr
    ." # CT" cr
    x . y . cr
    255 . cr
    y 0 do
	x 0 do
	    im j x 3 * * + i 3 * +
	    dup     c@ .
	    dup  1+ c@ .
	    2+ c@ .
	loop cr
    loop } ;


3dfile $@ r/o open-file throw value fid
!time
cube{{{ fid xdim ydim zdim lformat gulp->array
.time
fid close-file throw
cube{{{ xdim ydim * zdim max-arr constant gmax

[IFDEF] отладка
    gmax . cr
    cube{{{ xdim ydim * zdim min-arr . cr
    cr xdim . ydim . zdim . cr
    cr texdimzy . texdimxy . cr
    \            ydim 2/ 50  + zx->textr
    \            imtext-zx texdimxy texdimz textr->ppm
    \        xdim 2/  zy->textr
    \        imtext-zy texdimxy texdimz textr->ppm
    \    zdim 2/  xy->textr
    \    imtext-xy texdimxy dup  textr->ppm    
    bye    
[THEN]

\ *****************************************************************
\                              Action
\ *****************************************************************
[IFUNDEF] отладка
    zdim 2/ xy->textr    
    ydim 2/ zx->textr
    xdim 2/ zy->textr
    \ *****************************************************************
    \                             Cleanup
    \ *****************************************************************
[THEN]