3D Geometry , basics

Post your FreeBASIC source, examples, tips and tricks here. Please don’t post code without including an explanation.
Luxan
Posts: 222
Joined: Feb 18, 2009 12:47
Location: New Zealand

3D Geometry , basics

Post by Luxan »

These are the basic elements used when constructing simple 3D shapes.
Mostly intended for use with wire frame structures.

This code isn't heavily dependent upon external libraries.
Compatible with Windows , Linux and DOS.

Code: Select all

'
' -----------------------------------------------------------------------------
'
'   My graf 3d 
'
'    (c) copyright 2015 , sciwise@ihug.co.nz ,
'
'             Edward.Q.Montigue.  [ alias]
'
'
'
'
'
' -----------------------------------------------------------------------------
'
type point
         x as single
         y as single
         z as single
         u as single '  possible extension for special coord system
end type
'
const Pi=4*atn(1)
'
dim as single x1,y1,z1,x2,y2,z2
dim as integer i,j,k
'
'
'
dim as point p1(1 to 8)
dim as integer edge(1 to 12,0 to 1)
'
'                  Looking at a cube .
'
'             -1,1 _______<_______  1,1    start         z = -1
'                   |                              |             back face.
'                   |                              |
'                  v                             ^
'                   |                              |
'                   |_______________|
'                -1,-1        >        1,-1
'                
'
' -----------------------------------------------------------------------------
'
declare function rotx(q as point,angx as single) as point
declare function roty(q as point,angy as single) as point
declare function rotz(q as point,angz as single) as point
declare function tranx(q as point,movx as single) as point
declare function trany(q as point,movy as single) as point
declare function tranz(q as point,movz as single) as point
declare function persp(q as point,d as single) as point
'
declare function Trall( p1() as point,n as integer,edge() as integer, div as integer ) as integer
'
' ==============================
'
restore store1
for i=1 to 8
   read p1(i).x
   read p1(i).y
   read p1(i).z
next i
'
restore store2
for i=1 to 12
   read edge(i,0)
   read edge(i,1)
next i
'
' -----------------------------------------------------------------------------
'
screen 12
window (-1.5,-1.5)-(1.5,1.5)
line (-1.4,-1.4)-(1.4,1.4),11,b
'
'cls
Trall( p1() ,8,edge(), 32 )
sleep
end
'
' ===================================
'
'     vertex data , easier to keep track of
'  data when we use multiple data statements.
'
store1: 
data  1,1,1
data -1,1,1
data-1,-1,1
data 1,-1,1
data 1,1,-1
data -1,1,-1
data -1,-1,-1
data 1,-1,-1
'
'  edge data 
'
store2:
data 1,2
data 1,4
data 1,5
data 2,3
data 2,6
data 3,4
data 3,7
data 4,8
data 5,6
data 5,8
data 6,7
data 7,8
'
' -------------------------------------------------------------------------------
'
function rotx(q as point,angx as single) as point
'
'
'
static as point p
'
             p.x = q.x
             p.y= q.y*cos(angx)-sin(angx)*p.z
             p.z= q.z*cos(angx)+sin(angx)*p.y
'
              return p
'
end function 
'
' -----------------------------------------------------------------------------
'
function roty(q as point,angy as single) as point
'
'
'
static as point p
'
            p.x = sin(angy)*q.z + cos(angy)*q.x
            p.y = q.y
            p.z = cos(angy)*q.z -sin(angy)*q.x
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function rotz(q as point,angz as single) as point
'
'                         Rotate around z axis .
'
static as point p
'
            p.x = sin(angz)*q.y + cos(angz)*q.x
            p.y = cos(angz)*q.y-sin(angz)*q.x
            p.z = q.z
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function tranx(q as point,movx as single) as point
'
'              Translate point along x axis
'
static as point p
'
              p.x=q.x + movx
              p.y=q.y 
              p.z=q.z 
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function trany(q as point,movy as single) as point
'
'
'
static as point p
'
              p.x=q.x
              p.y=q.y + movy
              p.z=q.z 
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function tranz(q as point,movz as single) as point
'
'
'
static as point p
'
              p.x=q.x
              p.y=q.y
              p.z=q.z + movz
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function persp(q as point,d as single) as point
'
'     3d  perspective .  
'
'    Add 2 to the numerator when using any negative z value.
'
static as point p
'
     p.x = d*q.x/(q.z+2)
     p.y = d*q.y/(q.z+2)
     p.z = d
'
     return p
'
end function
'
' -----------------------------------------------------------------------------
'
function Trall( p1() as point,n as integer,edge() as integer, div as integer ) as integer
'
'  Translate and rotate all vertices .
'     as an animation ,  for  n  cycles .
'
'   With  div number of angle divisions .
'
static as point p2(1 to 8)
static as single theta,thi,x1,y1,z1,x2,y2,z2
static as integer i,j,k
static as integer i1,j1,k1
'
theta = Pi/div
'
for i=1 to n
  for j = 0 to div
  cls
       thi = j*theta
   for k = 1 to 8
     p2(k) = roty(p1(k),thi)
     p2(k)=persp(p2(k),0.8)
   next k     
'
for i1 = 1 to 12
      j1 = edge(i1,0)
     k1 = edge(i1,1)
   x1 = p2(j1).x
   y1 = p2(j1).y
'   z1 = p2(j1).z    
   x2 = p2(k1).x
   y2 = p2(k1).y
'   z2 = p2(k1).z    
line(x1,y1)-(x2,y2),14 
next i1     
'
sleep 100
  next j
next i
'
     return 0
'
end function
'
' --------------------------------------------------------------------------------
'






BasicCoder2
Posts: 3906
Joined: Jan 01, 2009 7:03
Location: Australia

Re: 3D Geometry , basics

Post by BasicCoder2 »

Have played around a bit myself although Dodicat has it down to a tee.
http://www.freebasic.net/forum/viewtopi ... it=3d+cube
http://www.freebasic.net/forum/viewtopi ... it=3d+cube
.
Luxan
Posts: 222
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: 3D Geometry , basics

Post by Luxan »

Thanks BasicCoder2 , that's an interesting piece of code.

At the moment my programming momentum is in a slightly different direction ;
however I shall return to this example at some stage.

I've come across other Free BASIC 3D code references , some however don't
lead to a valid link. In particular I'm seeking an efficient code for removing back faces;
other than the painters algorithm.

You will note my use of data statements for vertexes and edges , this is a primary
part of my code and differs from the example you provided.
Luxan
Posts: 222
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: 3D Geometry , basics

Post by Luxan »

This is a minor update , with the number of points, np , being
defined by a data value and the number of edges , ne , being
defined by a data value.

Code: Select all

'
' -----------------------------------------------------------------------------
'
'   My graf 3d 
'
'    (c) copyright 2015 , sciwise@ihug.co.nz ,
'
'             Edward.Q.Montague.  [ alias]
'
'
'
'
'
' -----------------------------------------------------------------------------
'
type point
         x as single
         y as single
         z as single
         u as single '  possible extension for special coord system
end type
'
const Pi=4*atn(1)
'
dim as single x1,y1,z1,x2,y2,z2
dim as integer i,j,k
'
' ----------------------------------------------------------------------------
'
'
'                  Looking at a cube .
'
'             -1,1 _______<_______  1,1    start         z = -1
'                   |                              |             back face.
'                   |                              |
'                  v                             ^
'                   |                              |
'                   |_______________|
'                -1,-1        >        1,-1
'                
'
' -----------------------------------------------------------------------------
'
declare function rotx(q as point,angx as single) as point
declare function roty(q as point,angy as single) as point
declare function rotz(q as point,angz as single) as point
declare function tranx(q as point,movx as single) as point
declare function trany(q as point,movy as single) as point
declare function tranz(q as point,movz as single) as point
declare function persp(q as point,d as single) as point
'
declare function Trall( p1() as point,n as integer,edge() as integer, div as integer ,np as integer,ne as integer) as integer
'
' =============================================
'
dim as integer np , ne
'
restore  storeA
read np
'
dim as point p1(1 to np)
'
restore store1
for i =1 to np
   read p1(i).x
   read p1(i).y
   read p1(i).z
next i
'
'
restore  storeB
read ne
'
dim as integer edge(1 to ne,0 to 1)
'
restore store2
for i = 1 to ne
   read edge(i,0)
   read edge(i,1)
next i
'
' -----------------------------------------------------------------------------
'
screen 12
window (-1.5,-1.5)-(1.5,1.5)
line (-1.4,-1.4)-(1.4,1.4),11,b
'
'cls
k=Trall( p1() ,3,edge() , 32 ,np ,ne ) 
sleep
end
'
' ===================================
'
'    number of vertices
'
storeA:
data 8  
'
'    vertex data , easier to keep track of
'  data when we use multiple data statements.
'
store1: 
data  1,1,1
data  -1,1,1
data  -1,-1,1
data 1,-1,1
data  1,1,-1
data  -1,1,-1
data  -1,-1,-1
data  1,-1,-1
'
'   Number of edges.
'
storeB:
data 12 
'
'  edge data 
'
store2:
data 1,2
data 1,4
data 1,5
data 2,3
data 2,6
data 3,4
data 3,7
data 4,8
data 5,6
data 5,8
data 6,7
data 7,8
'
' -------------------------------------------------------------------------------
'
function rotx(q as point,angx as single) as point
'
'              Rotate around x axis .
'
static as point p
'
             p.x = q.x
             p.y= q.y*cos(angx)-sin(angx)*q.z
             p.z= q.z*cos(angx)+sin(angx)*q.y
'
              return p
'
end function 
'
' -----------------------------------------------------------------------------
'
function roty(q as point,angy as single) as point
'
'              Rotate around y axis .
'
static as point p
'
            p.x = sin(angy)*q.z + cos(angy)*q.x
            p.y = q.y
            p.z = cos(angy)*q.z -sin(angy)*q.x
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function rotz(q as point,angz as single) as point
'
'              Rotate around z axis .
'
static as point p
'
            p.x = sin(angz)*q.y + cos(angz)*q.x
            p.y = cos(angz)*q.y-sin(angz)*q.x
            p.z = q.z
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function tranx(q as point,movx as single) as point
'
'             Translate point along x axis
'
static as point p
'
              p.x=q.x + movx
              p.y=q.y 
              p.z=q.z 
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function trany(q as point,movy as single) as point
'
'              Translate point along y axis .
'
static as point p
'
              p.x=q.x
              p.y=q.y + movy
              p.z=q.z 
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function tranz(q as point,movz as single) as point
'
'              Translate point along z axis .
'
static as point p
'
              p.x=q.x
              p.y=q.y
              p.z=q.z + movz
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function persp(q as point,d as single) as point
'
'              3d  perspective .  
'
'    Add 2 to the numerator when using any negative z value.
'
'     In this instance    -1 <= z  <= 1  ,  unit cube .
'
'    Therefore 2 is appropriate .
'
static as point p
'
     p.x = d*q.x/(q.z+2)
     p.y = d*q.y/(q.z+2)
     p.z = d
'
     return p
'
end function
'
' -----------------------------------------------------------------------------
'
function Trall( p1() as point,n as integer,edge() as integer, div as integer ,np as integer,ne as integer) as integer
'
'              Translate and rotate all vertices .
'     as an animation ,  for  n  cycles .
'
'    np number of points .
'    ne number of integers .
'
'
'   With  div number of angle divisions .
'
static as point p2(1 to np)
static as single theta,thi,x1,y1,z1,x2,y2,z2
static as integer i,j,k
static as integer i1,j1,k1
'
theta = Pi/div
'
for i=1 to n
  for j = 0 to div
  cls
       thi = j*theta
   for k = 1 to np
     p2(k) = roty(p1(k),thi)
     p2(k)=persp(p2(k),0.8)
   next k     
'
for i1 = 1 to ne
      j1 = edge(i1,0)
     k1 = edge(i1,1)
   x1 = p2(j1).x
   y1 = p2(j1).y
'   z1 = p2(j1).z    
   x2 = p2(k1).x
   y2 = p2(k1).y
'   z2 = p2(k1).z    
line(x1,y1)-(x2,y2),14 
next i1     
'
sleep 100
  next j
next i
'
     return 0
'
end function
'
' --------------------------------------------------------------------------------
'
Last edited by Luxan on Jan 26, 2016 6:37, edited 1 time in total.
BasicCoder2
Posts: 3906
Joined: Jan 01, 2009 7:03
Location: Australia

Re: 3D Geometry , basics

Post by BasicCoder2 »

Luxan wrote:You will note my use of data statements for vertexes and edges , this is a primary part of my code and differs from the example you provided.
So it is called the "painter's algorithm" I didn't know. Dodicat has some very fast 3D pixel setting of points in 3D space and some neat examples but much of his stuff is beyond me.

I assume this is the kind of thing you are looking at using?
https://en.wikipedia.org/wiki/Polygon_mesh

Not sure how far you can go with FreeBasic without using something like openGL to do the heavy lifting?

.
Luxan
Posts: 222
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: 3D Geometry , basics

Post by Luxan »

I went to that link and downloaded the page for future reference.

There are a number of directions that I want to go using the basic geometry
functions.

I need to change to another operating system to upload an example of one
such direction.
Luxan
Posts: 222
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: 3D Geometry , basics

Post by Luxan »

Okay , so this is a variation upon the theme of 3d geometry.
In this instance I'm graphing , with perspective , a large number
of points ; the maximum size being 128x128x128 .
A conditional statement is used to choose a range of magnitudes ;
the complete set of points is then rotated.

So what you say , however consider he number of points involved ;
few graphing packages can display this many and also animate
the data.
Adapting this to use 3d array values is also possible , in that instance we
need only store the magnitudes ; most likely in integer format.

A few aspects of this code require modification , the cube around the
data for instance.

Code: Select all

'
' --------------------------------------------------------------------
'
'               fn3d2.bas
'
'    (c) Copyright 2016  sciwise@ihug.co.nz
'
'        Luxan
'
'        Edward.Q.Montague  [ alias ]
'
'
'                 3d function 
'
'
'         Screen pages are used  .
'
'
' -------------------------------------------------------------------
'
type point
         x as single
         y as single
         z as single
         u as single '  possible extension for special coord system
end type
'
' -------------------------------------------------------------------
'
declare function roty(q as point,angy as single) as point
declare function persp(q as point,d as single) as point
'
declare function f3(x as single,y as single,z as single) as single
'
' =====================================
'
const Pi = 4*atn(1)

dim as point p1,p4

dim as integer n,m,p
dim as single x,y,z,u,maxa,thi,theta,div
'
 n = 128
m = 128
 p = 128
'
dim as single a1(0 to n-1,0 to m-1 , 0 to p-1)
dim as integer i,j,k,jt
'
' -------------------------------------------------------------
'
dim as long pal(15)
'
'     re defined 15 colour mode .
'
     pal(0)=5
     pal(1)=1
     pal(2)=9
     pal(3)=11
     pal(4)=2
     pal(5)=10
     pal(6)=14
     pal(7)=6
     pal(8)=4
     pal(9)=12
     pal(10)=13
     pal(11)=15
     pal(12)=0
     pal(13)=0
     pal(14)=0
     pal(15)=0
'
' ----------------------------------------------------------
'
dim as integer np , ne ,i1,j1,k1
dim as single x1,y1,x2,y2
'
restore  storeA
read np
'
dim as point p1a(1 to np),p2(1 to np)
'
restore store1
for i =1 to np
   read p1a(i).x
   read p1a(i).y
   read p1a(i).z
next i
'
'
restore  storeB
read ne
'
dim as integer edge(1 to ne,0 to 1)
'
restore store2
for i = 1 to ne
   read edge(i,0)
   read edge(i,1)
next i
'
' -----------------------------------------------------------
'
screen 12,8,2
'' image for working page #1 (visible page #0)
ScreenSet 1, 0 
Cls
window (-1.5,-1.5)-(1.5,1.5)
line (-1.4,-1.4)-(1.4,1.4),11,b
ScreenCopy 1, 0
'
' -----------------------------------------------------------
'
'   Draw cube frame and 
'
'   Generate from 3d function for each rotation .
'
' -----------------------------------------------------------
'
dim as single c1,d1,c2,d2
'
c1=1.2
d1=0.4
c2=1.35
d2=1.3
'
   div = 32
theta = Pi/div
'
for jt = 0 to div
  cls
'
'
for i = 0 to 11
     y1=(d2-d1)*i/15+d1
     y2=(d2-d1)*(i+1)/15+d1
      line(c1,y1)-(c2,y2),pal(i),bf    
     line(c1,y1)-(c2,y2),11,b
next i
'  
       thi = jt*theta
'       
   for k = 1 to np
     p2(k) = roty(p1a(k),thi)
     p2(k)=persp(p2(k),0.9)
   next k     
'
for i1 = 1 to ne
      j1 = edge(i1,0)
     k1 = edge(i1,1)
   x1 = p2(j1).x
   y1 = p2(j1).y
'   z1 = p2(j1).z    
   x2 = p2(k1).x
   y2 = p2(k1).y
'   z2 = p2(k1).z    
line(x1,y1)-(x2,y2),14 
next i1           
'
for k = 0 to p-1
      z = -1+2*k/(p-1)
 for j = 0 to m-1
      y = -1 + 2*j/(m-1)
  for i = 0 to n-1
       x = -1 + 2*i/(n-1)
       u =  f3(x,y,z)
       u = abs(u)
'
'   Set level where point becomes visible .
'
'   rotate and apply perspective .
'
if (u >= 0.1 and u<=0.15) then 
      p1.x = x
      p1.y = y
      p1.z = z
     p4 = roty(p1,thi)
     p4=persp(p4,0.9)      
     pset(p4.x,p4.y),pal(u*11)
end if
'
  next i
 next j
next  k
'
ScreenCopy 1, 0
'
'sleep 300
 next jt
'
end
'
' ====================================
'
'
'    number of vertices
'
storeA:
data 8  
'
'    vertex data , easier to keep track of
'  data when we use multiple data statements.
'
store1: 
data  1,1,1
data  -1,1,1
data -1,-1,1
data 1,-1,1
data 1,1,-1
data  -1,1,-1
data  -1,-1,-1
data  1,-1,-1
'
'   Number of edges.
'
storeB:
data 12 
'
'  edge data 
'
store2:
data 1,2
data 1,4
data 1,5
data 2,3
data 2,6
data 3,4
data 3,7
data 4,8
data 5,6
data 5,8
data 6,7
data 7,8
'
' -----------------------------------------------------------------------------
'
function roty(q as point,angy as single) as point
'
'              Rotate , a single point , around y axis .
'
static as point p
'
            p.x = sin(angy)*q.z + cos(angy)*q.x
            p.y = q.y
            p.z = cos(angy)*q.z -sin(angy)*q.x
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function persp(q as point,d as single) as point
'
'              3d  perspective .  
'
'     Applied to a single point .
'
'    Add 2 to the numerator when using any negative z value.
'
'     In this instance    -1 <= z  <= 1  ,  unit cube .
'
'    Therefore 2 is appropriate .
'
static as point p
'
     p.x = d*q.x/(q.z+2.5)
     p.y = d*q.y/(q.z+2.5)
     p.z = d
'
     return p
'
end function
'
' ---------------------------------------------------------------
'
function f3(x as single,y as single,z as single) as single
'
'        3d  sinc  function  
'
'
static as single u,v,w
'
    u=1
    if (abs(x)>0) then u= sin(4*x*Pi)/(x*4*Pi)
    v=1
    if (abs(y)>0) then v= sin(4*y*Pi)/(4*y*Pi)
    w=1
    if (abs(z)>0) then w= sin(4*z*Pi)/(4*z*Pi)
'
    return u*v*w
'
end function
'
' -----------------------------------------------------------
'
Luxan
Posts: 222
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: 3D Geometry , basics

Post by Luxan »

There were a few coding errors with the last post I sent.

Also I did a tidy up of the structure of the code , with the
use of more functions , my favorite way of getting some
understanding of what I'm attempting to do.

The function rotx was mis coded the rhs contained a
reference to the p variable , where q was expected .

The function pesp now has a larger numerical constant
in the divisor , this results in less extreme perspective.
I also increased the value of d , when calling the function
persp , to bring the object closer to the viewer.

I changed the edges of the cube that are displayed ,
the back faces are what are mostly dominant ;
however more coding is required with this.

The function that I'm displaying is now slightly different
with a different number of periods in the x,y,z directions.

I commented out the array a1() as this isn't used presently
and is something to implement in the future ; along
with a few other possible improvements.

Code: Select all

'
' --------------------------------------------------------------------
'
'               fn3d3.bas
'
'    (c) Copyright 2016  sciwise@ihug.co.nz
'
'        Luxan
'
'        Edward.Q.Montague  [ alias ]
'
'
'                 3d function 
'
'
'         Screen pages are used  .
'
'
' -------------------------------------------------------------------
'
type point
         x as single
         y as single
         z as single
         u as single '  possible extension for special coord system
end type
'
' -------------------------------------------------------------------
'
declare function roty(q as point,angy as single) as point
declare function rotx(q as point,angx as single) as point
declare function rotz(q as point,angz as single) as point
declare function persp(q as point,d as single) as point
declare function cproj(n as integer,m as integer,p as integer,p1a() as point ,np as integer,edge() as integer,ne as integer ,pal() as long , ng as integer , div as integer) as integer
declare  function pally(pal() as long) as integer
declare function cube( p1a() as point, ByRef np as integer, edge() as integer,ByRef ne as integer) as integer
'
declare function f3(x as single,y as single,z as single) as single
'
' =====================================
'
const Pi = 4*atn(1)
dim as integer n,m,p,fg
'
 n = 128
m = 128
 p = 128
'
'dim as single a1(0 to n-1,0 to m-1 , 0 to p-1)
' ----------------------------------------------------------
'
dim pal(15) as long
fg = pally(pal() ) 
'
' -----------------------------------------------------------
'
dim as integer np,ne
dim as point p1a()
dim as integer edge()
fg = cube(p1a() ,np ,edge() ,ne ) 
'
' -----------------------------------------------------------
'
screen 12,8,2
' image for working page #1 (visible page #0)
ScreenSet 1, 0 
Cls
window (-1.5,-1.5)-(1.5,1.5)
line (-1.4,-1.4)-(1.4,1.4),11,b
ScreenCopy 1, 0
'
fg = cproj(n,m,p,p1a(),np,edge() ,ne  ,pal()  , 4  , 32 ) 
'
end
'
' ====================================
'
'
'    number of vertices
'
storeA:
data 8  
'
'    vertex data , easier to keep track of
'  data when we use multiple data statements.
'
store1: 
data  1,1,1
data  -1,1,1
data -1,-1,1
data 1,-1,1
data 1,1,-1
data  -1,1,-1
data  -1,-1,-1
data  1,-1,-1
'
'   Number of edges.
'
storeB:
data 12 
'
'  edge data 
'
store2:
data 1,2
data 1,4
data 1,5
data 2,3
data 2,6
data 3,4
data 3,7
data 4,8
data 5,6
data 5,8
data 6,7
data 7,8
'
'   Number of edges.
'
storeC:
data 9 
'
'  edge data  , back faces .
'
store3:
data   1 , 2
 data   1 , 4
 data   2 , 3
 data   3 , 4
 data   2 , 6
 data   3  , 7
data    6 , 7
 data   7 , 8
 data   4 , 8
'
' -------------------------------------------------------------------------------
'
function rotx(q as point,angx as single) as point
'
'              Rotate around x axis .
'
static as point p
'
             p.x = q.x
             p.y= q.y*cos(angx)-sin(angx)*q.z
             p.z= q.z*cos(angx)+sin(angx)*q.y
'
              return p
'
end function 
'
' -----------------------------------------------------------------------------
'
function roty(q as point,angy as single) as point
'
'              Rotate , a single point , around y axis .
'
static as point p
'
            p.x = sin(angy)*q.z + cos(angy)*q.x
            p.y = q.y
            p.z = cos(angy)*q.z -sin(angy)*q.x
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function rotz(q as point,angz as single) as point
'
'              Rotate around z axis .
'
static as point p
'
            p.x = sin(angz)*q.y + cos(angz)*q.x
            p.y = cos(angz)*q.y-sin(angz)*q.x
            p.z = q.z
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function persp(q as point,d as single) as point
'
'              3d  perspective .  
'
'     Applied to a single point .
'
'    Add 2 to the numerator when using any negative z value.
'
'     In this instance    -1 <= z  <= 1  ,  unit cube .
'
'    Therefore 2 is appropriate .
'
static as point p
'
     p.x = d*q.x/(q.z+3.5)
     p.y = d*q.y/(q.z+3.5)
     p.z = d
'
     return p
'
end function
'
' ---------------------------------------------------------------
'
function pally(pal() as long) as integer
'
'          Assign  palette information .
'

'
'     re defined 15 colour mode .
'
     pal(0)=5
     pal(1)=1
     pal(2)=9
     pal(3)=11
     pal(4)=2
     pal(5)=10
     pal(6)=14
     pal(7)=6
     pal(8)=4
     pal(9)=12
     pal(10)=13
     pal(11)=15
     pal(12)=0
     pal(13)=0
     pal(14)=0
     pal(15)=0
'
'
        return 0
'
end function
'
' --------------------------------------------------------------
'
function cube( p1a() as point, ByRef np as integer, edge() as integer,ByRef ne as integer) as integer
'
'         Read Cube coordinate points and edge sequence data . 
'
static as integer i
'
restore  storeA
read np
'
redim  p1a(1 to np) 
'
restore store1
for i =1 to np
   read p1a(i).x
   read p1a(i).y
   read p1a(i).z
next i
'
'
restore  storeC
read ne
'
redim  edge(1 to ne,0 to 1)
'
restore store3
for i = 1 to ne
   read edge(i,0)
   read edge(i,1)
next i
'
     return  0
'
end function
'
' --------------------------------------------------------------
'
function cproj(n as integer,m as integer,p as integer,p1a() as point ,np as integer,edge() as integer,ne as integer, pal() as long , ng as integer , div as integer) as integer
'
' -----------------------------------------------------------
'
'   Draw cube frame and 
'
'   n  ,  number of points x direction
'   m ,  number of points y direction
'   p  ,  number of points z direction
'
'   p1a()    point array for cube 
'   np number of points for cube array , this is usually a constant value .
'   edge()   edge array for cube
'   ne number of edge pairs for edge() array , this is usually a constant value.
' 
'  pal() , palette array ,  this is a constant length .
'  
'  ng , the number of cycles for the complete 360 deg rotation. 
'
' div , sets the  size of each  discrete rotation  ; 2*Pi/div .    
'
'
'   Generate values from 3d function for each rotation .
'
' -----------------------------------------------------------
'
static as single c1,d1,c2,d2,theta
dim as single x,y,z,u,maxa,thi
dim as single x1,y1,z1,u1,x2,y2,z2,u2
static as integer i,j,k,ni,jt,i1,j1,k1
static as point p1,p4,p2(1 to np)
'
c1=1.2
d1=0.4
c2=1.35
d2=1.3
'
theta = Pi/div
'
for ni=0 to ng
'
for jt = 0 to div
  cls
'
'  ....   Draw color palette ,,,,,,,
'
for i = 0 to 11
     y1=(d2-d1)*i/15+d1
     y2=(d2-d1)*(i+1)/15+d1
      line(c1,y1)-(c2,y2),pal(i),bf    
     line(c1,y1)-(c2,y2),11,b
next i
'  
'  .......... Draw outer cube .............
'
       thi = jt*theta
'       
   for k = 1 to np
     p2(k) = roty(p1a(k),thi)
     p2(k) = rotx(p2(k),-0.75) 
     p2(k)=persp(p2(k),2)
   next k     
'
for i1 = 1 to ne
      j1 = edge(i1,0)
     k1 = edge(i1,1)
   x1 = p2(j1).x
   y1 = p2(j1).y
'   z1 = p2(j1).z    
   x2 = p2(k1).x
   y2 = p2(k1).y
'   z2 = p2(k1).z    
line(x1,y1)-(x2,y2),14 
next i1 
'
'  .................... Generate 3d point values .................         
'
for k = 0 to p-1
      z = -1+2*k/(p-1)
 for j = 0 to m-1
      y = -1 + 2*j/(m-1)
  for i = 0 to n-1
       x = -1 + 2*i/(n-1)
       u =  f3(x,y,z)
       u = abs(u)
'
'   Set level where point becomes visible .
'
'   rotate and apply perspective .
'
if (u >= 0.1 and u<=0.15) then 
      p1.x = x
      p1.y = y
      p1.z = z
     p4 = roty(p1,thi)
     p4 = rotx(p4,-0.75) 
     p4 = persp(p4,2)      
     pset(p4.x,p4.y),pal(u*11)
end if
'
  next i
 next j
next  k
'
ScreenCopy 1, 0
'
'sleep 300
 next jt
 '
 '    repeat cycle (ng - ni)  imes
 '
 next ni

'
      return 0
'
end function
'
' ----------------------------------------------------------------
'
function f3(x as single,y as single,z as single) as single
'
'        3d  sinc  function  
'
'
static as single u,v,w
'
    u=1
    if (abs(x)>0) then u= sin(3*x*Pi)/(x*3*Pi)
    v=1
    if (abs(y)>0) then v= sin(4*y*Pi)/(4*y*Pi)
    w=1
    if (abs(z)>0) then w= sin(5*z*Pi)/(5*z*Pi)
'
    return u*v*w
'
end function
'
' -----------------------------------------------------------
'

Luxan
Posts: 222
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: 3D Geometry , basics

Post by Luxan »

Another update , this occasion a selective data filter results in a significant
speed up when displaying a predefined data level .

Code: Select all


'
' --------------------------------------------------------------------
'
'               fn3d4.bas
'
'    (c) Copyright 2016  sciwise@ihug.co.nz
'
'        Luxan
'
'        Edward.Q.Montague  [ alias ]
'
'
'                 3d function 
'
'
'         Screen pages are used  .
'
'
' -------------------------------------------------------------------
'
type point
         x as single
         y as single
         z as single
         u as single '  possible extension for special coord system
end type
'
' -------------------------------------------------------------------
'
declare function roty(q as point,angy as single) as point
declare function rotx(q as point,angx as single) as point
declare function rotz(q as point,angz as single) as point
declare function persp(q as point,d as single) as point
declare function cproj(n as integer,m as integer,p as integer,p1a() as point ,np as integer,edge() as integer,ne as integer ,pal() as long , ng as integer , div as integer) as integer
declare  function pally(pal() as long) as integer
declare function cube( p1a() as point, ByRef np as integer, edge() as integer,ByRef ne as integer) as integer
'
declare function f3(x as single,y as single,z as single) as single
'
'
declare function dfilter(a2() as point , n as integer , m as integer , p as integer) as integer
declare function dproj(a2() as point , count as integer , p1a() as point ,np as integer,edge() as integer,ne as integer, pal() as long , ng as integer , div as integer) as integer
'
' =====================================
'
const Pi = 4*atn(1)
dim as integer n,m,p,fg
'
 n = 128
m = 128
 p = 128
'
'dim as single a1(0 to n-1,0 to m-1 , 0 to p-1)
' ----------------------------------------------------------
'
dim pal(15) as long
fg = pally(pal() ) 
'
' -----------------------------------------------------------
'
dim as integer np,ne
dim as point p1a()
dim as integer edge()
fg = cube(p1a() ,np ,edge() ,ne ) 
'
' -----------------------------------------------------------
'
screen 12,8,2
' image for working page #1 (visible page #0)
ScreenSet 1, 0 
Cls
window (-1.5,-1.5)-(1.5,1.5)
line (-1.4,-1.4)-(1.4,1.4),11,b
ScreenCopy 1, 0
'
redim a2(0 to 4) as point 

fg = dfilter(a2()  , n  , m  , p ) 

fg= dproj(a2()  , fg , p1a()  ,np ,edge() ,ne , pal()  , 4  , 32 ) 



'fg = cproj(n,m,p,p1a(),np,edge() ,ne  ,pal()  , 4  , 32 ) 
'


'redim a2(4) as point
end
'
' ====================================
'
'
'    number of vertices
'
storeA:
data 8  
'
'    vertex data , easier to keep track of
'  data when we use multiple data statements.
'
store1: 
data  1,1,1
data  -1,1,1
data -1,-1,1
data 1,-1,1
data 1,1,-1
data  -1,1,-1
data  -1,-1,-1
data  1,-1,-1
'
'   Number of edges.
'
storeB:
data 12 
'
'  edge data 
'
store2:
data 1,2
data 1,4
data 1,5
data 2,3
data 2,6
data 3,4
data 3,7
data 4,8
data 5,6
data 5,8
data 6,7
data 7,8
'
'   Number of edges.
'
storeC:
data 9 
'
'  edge data  , back faces .
'
store3:
data   1 , 2
 data   1 , 4
 data   2 , 3
 data   3 , 4
 data   2 , 6
 data   3  , 7
data    6 , 7
 data   7 , 8
 data   4 , 8
'
' -------------------------------------------------------------------------------
'
function rotx(q as point,angx as single) as point
'
'              Rotate around x axis .
'
static as point p
'
             p.x = q.x
             p.y= q.y*cos(angx)-sin(angx)*q.z
             p.z= q.z*cos(angx)+sin(angx)*q.y
'
              return p
'
end function 
'
' -----------------------------------------------------------------------------
'
function roty(q as point,angy as single) as point
'
'              Rotate , a single point , around y axis .
'
static as point p
'
            p.x = sin(angy)*q.z + cos(angy)*q.x
            p.y = q.y
            p.z = cos(angy)*q.z -sin(angy)*q.x
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function rotz(q as point,angz as single) as point
'
'              Rotate around z axis .
'
static as point p
'
            p.x = sin(angz)*q.y + cos(angz)*q.x
            p.y = cos(angz)*q.y-sin(angz)*q.x
            p.z = q.z
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function persp(q as point,d as single) as point
'
'              3d  perspective .  
'
'     Applied to a single point .
'
'    Add 2 to the numerator when using any negative z value.
'
'     In this instance    -1 <= z  <= 1  ,  unit cube .
'
'    Therefore 2 is appropriate .
'
static as point p
'
     p.x = d*q.x/(q.z+3.5)
     p.y = d*q.y/(q.z+3.5)
     p.z = d
'
     return p
'
end function
'
' ---------------------------------------------------------------
'
function pally(pal() as long) as integer
'
'          Assign  palette information .
'

'
'     re defined 15 colour mode .
'
     pal(0)=5
     pal(1)=1
     pal(2)=9
     pal(3)=11
     pal(4)=2
     pal(5)=10
     pal(6)=14
     pal(7)=6
     pal(8)=4
     pal(9)=12
     pal(10)=13
     pal(11)=15
     pal(12)=0
     pal(13)=0
     pal(14)=0
     pal(15)=0
'
'
        return 0
'
end function
'
' --------------------------------------------------------------
'
function cube( p1a() as point, ByRef np as integer, edge() as integer,ByRef ne as integer) as integer
'
'         Read Cube coordinate points and edge sequence data . 
'
static as integer i
'
restore  storeA
read np
'
redim  p1a(1 to np) 
'
restore store1
for i =1 to np
   read p1a(i).x
   read p1a(i).y
   read p1a(i).z
next i
'
'
restore  storeC
read ne
'
redim  edge(1 to ne,0 to 1)
'
restore store3
for i = 1 to ne
   read edge(i,0)
   read edge(i,1)
next i
'
     return  0
'
end function
'
' --------------------------------------------------------------
'
function cproj(n as integer,m as integer,p as integer,p1a() as point ,np as integer,edge() as integer,ne as integer, pal() as long , ng as integer , div as integer) as integer
'
' -----------------------------------------------------------
'
'   Draw cube frame and 
'
'   n  ,  number of points x direction
'   m ,  number of points y direction
'   p  ,  number of points z direction
'
'   p1a()    point array for cube 
'   np number of points for cube array , this is usually a constant value .
'   edge()   edge array for cube
'   ne number of edge pairs for edge() array , this is usually a constant value.
' 
'  pal() , palette array ,  this is a constant length .
'  
'  ng , the number of cycles for the complete 360 deg rotation. 
'
' div , sets the  size of each  discrete rotation  ; 2*Pi/div .    
'
'
'   Generate values from 3d function for each rotation .
'
' -----------------------------------------------------------
'
static as single c1,d1,c2,d2,theta
dim as single x,y,z,u,maxa,thi
dim as single x1,y1,z1,u1,x2,y2,z2,u2
static as integer i,j,k,ni,jt,i1,j1,k1
static as point p1,p4,p2(1 to np)
'
c1=1.2
d1=0.4
c2=1.35
d2=1.3
'
theta = Pi/div
'
for ni=0 to ng
'
for jt = 0 to div
  cls
'
'  ....   Draw color palette ,,,,,,,
'
for i = 0 to 11
     y1=(d2-d1)*i/15+d1
     y2=(d2-d1)*(i+1)/15+d1
      line(c1,y1)-(c2,y2),pal(i),bf    
     line(c1,y1)-(c2,y2),11,b
next i
'  
'  .......... Draw outer cube .............
'
       thi = jt*theta
'       
   for k = 1 to np
     p2(k) = roty(p1a(k),thi)
     p2(k) = rotx(p2(k),-0.75) 
     p2(k)=persp(p2(k),2)
   next k     
'
for i1 = 1 to ne
      j1 = edge(i1,0)
     k1 = edge(i1,1)
   x1 = p2(j1).x
   y1 = p2(j1).y
'   z1 = p2(j1).z    
   x2 = p2(k1).x
   y2 = p2(k1).y
'   z2 = p2(k1).z    
line(x1,y1)-(x2,y2),14 
next i1 
'
'  .................... Generate 3d point values .................         
'
for k = 0 to p-1
      z = -1+2*k/(p-1)
 for j = 0 to m-1
      y = -1 + 2*j/(m-1)
  for i = 0 to n-1
       x = -1 + 2*i/(n-1)
       u =  f3(x,y,z)
       u = abs(u)
'
'   Set level where point becomes visible .
'
'   rotate and apply perspective .
'
if (u >= 0.1 and u<=0.15) then 
      p1.x = x
      p1.y = y
      p1.z = z
     p4 = roty(p1,thi)
     p4 = rotx(p4,-0.75) 
     p4 = persp(p4,2)      
     pset(p4.x,p4.y),pal(u*11)
end if
'
  next i
 next j
next  k
'
ScreenCopy 1, 0
'
'sleep 300
 next jt
 '
 '    repeat cycle (ng - ni)  imes
 '
 next ni

'
      return 0
'
end function
'
' ----------------------------------------------------------------
'
function dfilter(a2() as point , n as integer , m as integer , p as integer) as integer
'
'    Filter data  from function and selectively store in array a2()
'
'
static as integer i,j,k,count,cn
static as single x,y,z,u
static as point p1
'
'
'  .................... Generate 3d point values .................         
'
count = 0
for k = 0 to p-1
      z = -1+2*k/(p-1)
 for j = 0 to m-1
      y = -1 + 2*j/(m-1)
  for i = 0 to n-1
       x = -1 + 2*i/(n-1)
       u =  f3(x,y,z)
       u = abs(u)
'
'   Set level where point becomes visible .
'
'   rotate and apply perspective .
'
if (u >= 0.1 and u<=0.15) then 
     count = count + 1
'      p1.x = x
'      p1.y = y
'      p1.z = z

end if
'
  next i
 next j
next  k
'
if (count > 0 ) then
'
     redim as point  a2(count)
'
    cn =  0
for k = 0 to p-1
      z = -1+2*k/(p-1)
 for j = 0 to m-1
      y = -1 + 2*j/(m-1)
  for i = 0 to n-1
       x = -1 + 2*i/(n-1)
       u =  f3(x,y,z)
       u = abs(u)
'
'   Set level where point becomes visible .
'
'   rotate and apply perspective .
'
if (u >= 0.1 and u<=0.15) then   
      p1.x = x
      p1.y = y
      p1.z = z
  a2(cn) = p1
        cn = cn + 1
end if
'
  next i
 next j
next  k
'
end if
'
      return count
'
end function
'
' ----------------------------------------------------------------
'
function dproj(a2() as point , count as integer , p1a() as point ,np as integer,edge() as integer,ne as integer, pal() as long , ng as integer , div as integer) as integer
'
' -----------------------------------------------------------
'
'   Draw cube frame and 
'
'   count  ,  number of points  in array a2()
'  
'
'   p1a()    point array for cube 
'   np number of points for cube array , this is usually a constant value .
'   edge()   edge array for cube
'   ne number of edge pairs for edge() array , this is usually a constant value.
' 
'  pal() , palette array ,  this is a constant length .
'  
'  ng , the number of cycles for the complete 360 deg rotation. 
'
' div , sets the  size of each  discrete rotation  ; 2*Pi/div .    
'
'
'   Generate values from 3d function for each rotation .
'
' -----------------------------------------------------------
'
static as single c1,d1,c2,d2,theta
dim as single x,y,z,u,maxa,thi
dim as single x1,y1,z1,u1,x2,y2,z2,u2
static as integer i,j,k,ni,jt,i1,j1,k1
static as point p1,p4,p2(1 to np)
'
c1=1.2
d1=0.4
c2=1.35
d2=1.3
'
theta = Pi/div
'
for ni=0 to ng
'
for jt = 0 to div
  cls
'
'  ....   Draw color palette ,,,,,,,
'
for i = 0 to 11
     y1=(d2-d1)*i/15+d1
     y2=(d2-d1)*(i+1)/15+d1
      line(c1,y1)-(c2,y2),pal(i),bf    
     line(c1,y1)-(c2,y2),11,b
next i
'  
'  .......... Draw outer cube .............
'
       thi = jt*theta
'       
   for k = 1 to np
     p2(k) = roty(p1a(k),thi)
     p2(k) = rotx(p2(k),-0.75) 
     p2(k)=persp(p2(k),2)
   next k     
'
for i1 = 1 to ne
      j1 = edge(i1,0)
     k1 = edge(i1,1)
   x1 = p2(j1).x
   y1 = p2(j1).y
'   z1 = p2(j1).z    
   x2 = p2(k1).x
   y2 = p2(k1).y
'   z2 = p2(k1).z    
line(x1,y1)-(x2,y2),14 
next i1 
'
'  ....................  3d point values from array a2()  .................         
'
for k=0 to count
      p1 = a2(k)
     p4 = roty(p1,thi)
     p4 = rotx(p4,-0.75) 
     p4 = persp(p4,2)      
     pset(p4.x,p4.y),pal(u*11)
next k
'
ScreenCopy 1, 0
'
'sleep 300
 next jt
 '
 '    repeat cycle (ng - ni)  imes
 '
 next ni
'
      return 0
'
end function
'
' ----------------------------------------------------------------
'
function f3(x as single,y as single,z as single) as single
'
'        3d  sinc  function  
'
'
static as single u,v,w
'
    u=1
    if (abs(x)>0) then u= sin(3*x*Pi)/(x*3*Pi)
    v=1
    if (abs(y)>0) then v= sin(4*y*Pi)/(4*y*Pi)
    w=1
    if (abs(z)>0) then w= sin(5*z*Pi)/(5*z*Pi)
'
    return u*v*w
'
end function
'
' -----------------------------------------------------------
'




Luxan
Posts: 222
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: 3D Geometry , basics

Post by Luxan »

This is a minor update , as such I'm calling the code fn3d4a.bas.

In this instance I'm using the .u component of the point variable
array a2() to store the magnitude value of the point at (x,y,z).

The levels , which select a range of values for the magnitude , are now defined by
two variables l1 and l2 ; these are passed to the function dfilter ,

The function cproj has been removed in this version as dproj appears to
be more relevant.

Code: Select all

'
' --------------------------------------------------------------------
'
'               fn3d4a.bas
'
'    (c) Copyright 2016  sciwise@ihug.co.nz
'
'        Luxan
'
'        Edward.Q.Montague  [ alias ]
'
'
'                 3d function 
'
'
'         Screen pages are used  .
'
'
' -------------------------------------------------------------------
'
type point
         x as single
         y as single
         z as single
         u as single '  possible extension for special coord system
end type
'
' -------------------------------------------------------------------
'
declare function roty(q as point,angy as single) as point
declare function rotx(q as point,angx as single) as point
declare function rotz(q as point,angz as single) as point
declare function persp(q as point,d as single) as point
declare  function pally(pal() as long) as integer
declare function cube( p1a() as point, ByRef np as integer, edge() as integer,ByRef ne as integer) as integer
'
declare function f3(x as single,y as single,z as single) as single
'
'
declare function dfilter(a2() as point , n as integer , m as integer , p as integer,l1 as single,l2 as single) as integer
declare function dproj(a2() as point , count as integer , p1a() as point ,np as integer,edge() as integer,ne as integer, pal() as long , ng as integer , div as integer) as integer
'
' =====================================
'
const Pi = 4*atn(1)
dim as integer n,m,p,fg
'
 n = 128
m = 128
 p = 128
'
'dim as single a1(0 to n-1,0 to m-1 , 0 to p-1)
' ----------------------------------------------------------
'
dim pal(15) as long
fg = pally(pal() ) 
'
' -----------------------------------------------------------
'
dim as integer np,ne
dim as point p1a()
dim as integer edge()
dim as single l1,l2   '   levels
fg = cube(p1a() ,np ,edge() ,ne ) 
'
' -----------------------------------------------------------
'
screen 12,8,2
' image for working page #1 (visible page #0)
ScreenSet 1, 0 
Cls
window (-1.5,-1.5)-(1.5,1.5)
line (-1.4,-1.4)-(1.4,1.4),11,b
ScreenCopy 1, 0
'
redim a2(0 to 4) as point 

l1=0.1
l2=0.15
fg = dfilter(a2()  , n  , m  , p , l1 , l2 ) 

fg= dproj(a2()  , fg , p1a()  ,np ,edge() ,ne , pal()  , 4  , 32 ) 

redim a2(4) as point
end
'
' ====================================
'
'
'    number of vertices
'
storeA:
data 8  
'
'    vertex data , easier to keep track of
'  data when we use multiple data statements.
'
store1: 
data  1,1,1
data  -1,1,1
data -1,-1,1
data 1,-1,1
data 1,1,-1
data  -1,1,-1
data  -1,-1,-1
data  1,-1,-1
'
'   Number of edges.
'
storeB:
data 12 
'
'  edge data 
'
store2:
data 1,2
data 1,4
data 1,5
data 2,3
data 2,6
data 3,4
data 3,7
data 4,8
data 5,6
data 5,8
data 6,7
data 7,8
'
'   Number of edges.
'
storeC:
data 9 
'
'  edge data  , back faces .
'
store3:
data   1 , 2
 data   1 , 4
 data   2 , 3
 data   3 , 4
 data   2 , 6
 data   3  , 7
data    6 , 7
 data   7 , 8
 data   4 , 8
'
' -------------------------------------------------------------------------------
'
function rotx(q as point,angx as single) as point
'
'              Rotate around x axis .
'
static as point p
'
             p.x = q.x
             p.y= q.y*cos(angx)-sin(angx)*q.z
             p.z= q.z*cos(angx)+sin(angx)*q.y
'
              return p
'
end function 
'
' -----------------------------------------------------------------------------
'
function roty(q as point,angy as single) as point
'
'              Rotate , a single point , around y axis .
'
static as point p
'
            p.x = sin(angy)*q.z + cos(angy)*q.x
            p.y = q.y
            p.z = cos(angy)*q.z -sin(angy)*q.x
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function rotz(q as point,angz as single) as point
'
'              Rotate around z axis .
'
static as point p
'
            p.x = sin(angz)*q.y + cos(angz)*q.x
            p.y = cos(angz)*q.y-sin(angz)*q.x
            p.z = q.z
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function persp(q as point,d as single) as point
'
'              3d  perspective .  
'
'     Applied to a single point .
'
'    Add 2 to the numerator when using any negative z value.
'
'     In this instance    -1 <= z  <= 1  ,  unit cube .
'
'    Therefore 2 is appropriate .
'
static as point p
'
     p.x = d*q.x/(q.z+3.5)
     p.y = d*q.y/(q.z+3.5)
     p.z = d
'
     return p
'
end function
'
' ---------------------------------------------------------------
'
function pally(pal() as long) as integer
'
'          Assign  palette information .
'

'
'     re defined 15 colour mode .
'
     pal(0)=5
     pal(1)=1
     pal(2)=9
     pal(3)=11
     pal(4)=2
     pal(5)=10
     pal(6)=14
     pal(7)=6
     pal(8)=4
     pal(9)=12
     pal(10)=13
     pal(11)=15
     pal(12)=0
     pal(13)=0
     pal(14)=0
     pal(15)=0
'
'
        return 0
'
end function
'
' --------------------------------------------------------------
'
function cube( p1a() as point, ByRef np as integer, edge() as integer,ByRef ne as integer) as integer
'
'         Read Cube coordinate points and edge sequence data . 
'
static as integer i
'
restore  storeA
read np
'
redim  p1a(1 to np) 
'
restore store1
for i =1 to np
   read p1a(i).x
   read p1a(i).y
   read p1a(i).z
next i
'
'
restore  storeC
read ne
'
redim  edge(1 to ne,0 to 1)
'
restore store3
for i = 1 to ne
   read edge(i,0)
   read edge(i,1)
next i
'
     return  0
'
end function
'
' --------------------------------------------------------------
'

'
' ----------------------------------------------------------------
'
function dfilter(a2() as point , n as integer , m as integer , p as integer,l1 as single,l2 as single) as integer
'
'    Filter data  from function and selectively store in array a2()
'
'
static as integer i,j,k,count,cn
static as single x,y,z,u
static as point p1
'
'
u=l1
if l1 > l2 then
   l1=l2
   l2=u
end if
'
'
'  .................... Generate 3d point values .................         
'
count = 0
for k = 0 to p-1
      z = -1+2*k/(p-1)
 for j = 0 to m-1
      y = -1 + 2*j/(m-1)
  for i = 0 to n-1
       x = -1 + 2*i/(n-1)
       u =  f3(x,y,z)
       u = abs(u)
'
'   Set level where point becomes visible .
'
'   rotate and apply perspective .
'
'if (u >= 0.1 and u<=0.15) then 
if (u >= l1 and u<=l2) then 
     count = count + 1
'      p1.x = x
'      p1.y = y
'      p1.z = z

end if
'
  next i
 next j
next  k
'
if (count > 0 ) then
'
     redim as point  a2(count)
'
    cn =  0
for k = 0 to p-1
      z = -1+2*k/(p-1)
 for j = 0 to m-1
      y = -1 + 2*j/(m-1)
  for i = 0 to n-1
       x = -1 + 2*i/(n-1)
       u =  f3(x,y,z)
       u = abs(u)
'
'   Set level where point becomes visible .
'
'   rotate and apply perspective .
'
if (u >= 0.1 and u<=0.15) then   
      p1.x = x
      p1.y = y
      p1.z = z
      p1.u = u
  a2(cn) = p1
        cn = cn + 1
end if
'
  next i
 next j
next  k
'
end if
'
      return count
'
end function
'
' ----------------------------------------------------------------
'
function dproj(a2() as point , count as integer , p1a() as point ,np as integer,edge() as integer,ne as integer, pal() as long , ng as integer , div as integer) as integer
'
' -----------------------------------------------------------
'
'   Draw cube frame and 
'
'   count  ,  number of points  in array a2()
'  
'
'   p1a()    point array for cube 
'   np number of points for cube array , this is usually a constant value .
'   edge()   edge array for cube
'   ne number of edge pairs for edge() array , this is usually a constant value.
' 
'  pal() , palette array ,  this is a constant length .
'  
'  ng , the number of cycles for the complete 360 deg rotation. 
'
' div , sets the  size of each  discrete rotation  ; 2*Pi/div .    
'
'
'   Generate values from 3d function for each rotation .
'
' -----------------------------------------------------------
'
static as single c1,d1,c2,d2,theta
dim as single x,y,z,u,maxa,thi
dim as single x1,y1,z1,u1,x2,y2,z2,u2
static as integer i,j,k,ni,jt,i1,j1,k1
static as point p1,p4,p2(1 to np)
'
c1=1.2
d1=0.4
c2=1.35
d2=1.3
'
theta = Pi/div
'
for ni=0 to ng
'
for jt = 0 to div
  cls
'
'  ....   Draw color palette ,,,,,,,
'
for i = 0 to 11
     y1=(d2-d1)*i/15+d1
     y2=(d2-d1)*(i+1)/15+d1
      line(c1,y1)-(c2,y2),pal(i),bf    
     line(c1,y1)-(c2,y2),11,b
next i
'  
'  .......... Draw outer cube .............
'
       thi = jt*theta
'       
   for k = 1 to np
     p2(k) = roty(p1a(k),thi)
     p2(k) = rotx(p2(k),-0.75) 
     p2(k)=persp(p2(k),2)
   next k     
'
for i1 = 1 to ne
      j1 = edge(i1,0)
     k1 = edge(i1,1)
   x1 = p2(j1).x
   y1 = p2(j1).y
'   z1 = p2(j1).z    
   x2 = p2(k1).x
   y2 = p2(k1).y
'   z2 = p2(k1).z    
line(x1,y1)-(x2,y2),14 
next i1 
'
'  ....................  3d point values from array a2()  .................         
'
for k=0 to count
      p1 = a2(k)
      u = p1.u
     p4 = roty(p1,thi)
     p4 = rotx(p4,-0.75) 
     p4 = persp(p4,2)      
     pset(p4.x,p4.y),pal(u*11)
next k
'
ScreenCopy 1, 0
'
sleep 100
 next jt
 '
 '    repeat cycle (ng - ni)  imes
 '
 next ni
'
      return 0
'
end function
'
' ----------------------------------------------------------------
'
function f3(x as single,y as single,z as single) as single
'
'        3d  sinc  function  
'
'
static as single u,v,w
'
    u=1
    if (abs(x)>0) then u= sin(3*x*Pi)/(x*3*Pi)
    v=1
    if (abs(y)>0) then v= sin(4*y*Pi)/(4*y*Pi)
    w=1
    if (abs(z)>0) then w= sin(5*z*Pi)/(5*z*Pi)
'
    return u*v*w
'
end function
'
' -----------------------------------------------------------
'

Luxan
Posts: 222
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: 3D Geometry , basics

Post by Luxan »

Now we use the routines thus far developed to select and display data
from a 3d array .

I consider this as being a more than minor update , therefore I'm calling
the code fn3d5.bas .

I removed the function dfilter , as this isn't compatible with the data
generated for array a1() .

I also made certain that the level selection is made within both iteration loops
of efilter .

Code: Select all


'
' --------------------------------------------------------------------
'
'               fn3d5.bas
'
'    (c) Copyright 2016  sciwise@ihug.co.nz
'
'        Luxan
'
'        Edward.Q.Montague  [ alias ]
'
'
'                 3d function 
'
'
'         Screen pages are used  .
'
'
' -------------------------------------------------------------------
'
type point
         x as single
         y as single
         z as single
         u as single '  possible extension for special coord system
end type
'
' -------------------------------------------------------------------
'
declare function roty(q as point,angy as single) as point
declare function rotx(q as point,angx as single) as point
declare function rotz(q as point,angz as single) as point
declare function persp(q as point,d as single) as point
declare  function pally(pal() as long) as integer
declare function cube( p1a() as point, ByRef np as integer, edge() as integer,ByRef ne as integer) as integer
'
declare function f3(x as single,y as single,z as single) as single
declare function  Generate(a1() as single,n as integer,m as integer,p as integer) as integer
declare function efilter(a1() as single , a2() as point , n as integer , m as integer , p as integer,l1 as single,l2 as single) as integer
'
declare function dproj(a2() as point , count as integer , p1a() as point ,np as integer,edge() as integer,ne as integer, pal() as long , ng as integer , div as integer) as integer
'
' =====================================
'
const Pi = 4*atn(1)
dim as integer n,m,p,fg
'
 n = 128
m = 128
 p = 128
'
redim as single a1(0 to 4,0 to 4 , 0 to 4)
' ----------------------------------------------------------
'
dim pal(15) as long
fg = pally(pal() ) 
'
' -----------------------------------------------------------
'
dim as integer np,ne
dim as point p1a()
dim as integer edge()
dim as single l1,l2   '   levels
fg = cube(p1a() ,np ,edge() ,ne ) 
'
' -----------------------------------------------------------
'
screen 12,8,2
' image for working page #1 (visible page #0)
ScreenSet 1, 0 
Cls
window (-1.5,-1.5)-(1.5,1.5)
line (-1.4,-1.4)-(1.4,1.4),11,b
ScreenCopy 1, 0
'
redim a2(0 to 4) as point 
'
' --------------------------------------------------------------
'
fg = Generate(a1(),n ,m ,p ) 
'
l1=0.1
l2=0.15
fg= efilter(a1()  , a2()  , n  , m  , p ,l1 ,l2 ) 
'
fg= dproj(a2()  , fg , p1a()  ,np ,edge() ,ne , pal()  , 4  , 32 ) 
'
redim a2(4) as point
redim as single a1(0 to 4,0 to 4 , 0 to 4)
'
sleep
'
end
'
' ====================================
'
'
'    number of vertices
'
storeA:
data 8  
'
'    vertex data , easier to keep track of
'  data when we use multiple data statements.
'
store1: 
data  1,1,1
data  -1,1,1
data -1,-1,1
data 1,-1,1
data 1,1,-1
data  -1,1,-1
data  -1,-1,-1
data  1,-1,-1
'
'   Number of edges.
'
storeB:
data 12 
'
'  edge data 
'
store2:
data 1,2
data 1,4
data 1,5
data 2,3
data 2,6
data 3,4
data 3,7
data 4,8
data 5,6
data 5,8
data 6,7
data 7,8
'
'   Number of edges.
'
storeC:
data 9 
'
'  edge data  , back faces .
'
store3:
data   1 , 2
 data   1 , 4
 data   2 , 3
 data   3 , 4
 data   2 , 6
 data   3  , 7
data    6 , 7
 data   7 , 8
 data   4 , 8
'
' -------------------------------------------------------------------------------
'
function rotx(q as point,angx as single) as point
'
'              Rotate around x axis .
'
static as point p
'
             p.x = q.x
             p.y= q.y*cos(angx)-sin(angx)*q.z
             p.z= q.z*cos(angx)+sin(angx)*q.y
'
              return p
'
end function 
'
' -----------------------------------------------------------------------------
'
function roty(q as point,angy as single) as point
'
'              Rotate , a single point , around y axis .
'
static as point p
'
            p.x = sin(angy)*q.z + cos(angy)*q.x
            p.y = q.y
            p.z = cos(angy)*q.z -sin(angy)*q.x
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function rotz(q as point,angz as single) as point
'
'              Rotate around z axis .
'
static as point p
'
            p.x = sin(angz)*q.y + cos(angz)*q.x
            p.y = cos(angz)*q.y-sin(angz)*q.x
            p.z = q.z
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function persp(q as point,d as single) as point
'
'              3d  perspective .  
'
'     Applied to a single point .
'
'    Add 2 to the numerator when using any negative z value.
'
'     In this instance    -1 <= z  <= 1  ,  unit cube .
'
'    Therefore 2 is appropriate .
'
static as point p
'
     p.x = d*q.x/(q.z+3.5)
     p.y = d*q.y/(q.z+3.5)
     p.z = d
'
     return p
'
end function
'
' ---------------------------------------------------------------
'
function pally(pal() as long) as integer
'
'          Assign  palette information .
'

'
'     re defined 15 colour mode .
'
     pal(0)=5
     pal(1)=1
     pal(2)=9
     pal(3)=11
     pal(4)=2
     pal(5)=10
     pal(6)=14
     pal(7)=6
     pal(8)=4
     pal(9)=12
     pal(10)=13
     pal(11)=15
     pal(12)=0
     pal(13)=0
     pal(14)=0
     pal(15)=0
'
'
        return 0
'
end function
'
' --------------------------------------------------------------
'
function cube( p1a() as point, ByRef np as integer, edge() as integer,ByRef ne as integer) as integer
'
'         Read Cube coordinate points and edge sequence data . 
'
static as integer i
'
restore  storeA
read np
'
redim  p1a(1 to np) 
'
restore store1
for i =1 to np
   read p1a(i).x
   read p1a(i).y
   read p1a(i).z
next i
'
'
restore  storeC
read ne
'
redim  edge(1 to ne,0 to 1)
'
restore store3
for i = 1 to ne
   read edge(i,0)
   read edge(i,1)
next i
'
     return  0
'
end function
'
' ----------------------------------------------------------------
'
function dproj(a2() as point , count as integer , p1a() as point ,np as integer,edge() as integer,ne as integer, pal() as long , ng as integer , div as integer) as integer
'
' -----------------------------------------------------------
'
'   Draw cube frame and 
'
'   count  ,  number of points  in array a2()
'  
'
'   p1a()    point array for cube 
'   np number of points for cube array , this is usually a constant value .
'   edge()   edge array for cube
'   ne number of edge pairs for edge() array , this is usually a constant value.
' 
'  pal() , palette array ,  this is a constant length .
'  
'  ng , the number of cycles for the complete 360 deg rotation. 
'
' div , sets the  size of each  discrete rotation  ; 2*Pi/div .    
'
'
'   Generate values from 3d function for each rotation .
'
' -----------------------------------------------------------
'
static as single c1,d1,c2,d2,theta
dim as single x,y,z,u,maxa,thi
dim as single x1,y1,z1,u1,x2,y2,z2,u2
static as integer i,j,k,ni,jt,i1,j1,k1
static as point p1,p4,p2(1 to np)
'
c1=1.2
d1=0.4
c2=1.35
d2=1.3
'
theta = Pi/div
'
for ni=0 to ng
'
for jt = 0 to div
  cls
'
'  ....   Draw color palette ,,,,,,,
'
for i = 0 to 11
     y1=(d2-d1)*i/15+d1
     y2=(d2-d1)*(i+1)/15+d1
      line(c1,y1)-(c2,y2),pal(i),bf    
     line(c1,y1)-(c2,y2),11,b
next i
'  
'  .......... Draw outer cube .............
'
       thi = jt*theta
'       
   for k = 1 to np
     p2(k) = roty(p1a(k),thi)
     p2(k) = rotx(p2(k),-0.75) 
     p2(k)=persp(p2(k),2)
   next k     
'
for i1 = 1 to ne
      j1 = edge(i1,0)
     k1 = edge(i1,1)
   x1 = p2(j1).x
   y1 = p2(j1).y
'   z1 = p2(j1).z    
   x2 = p2(k1).x
   y2 = p2(k1).y
'   z2 = p2(k1).z    
line(x1,y1)-(x2,y2),14 
next i1 
'
'  ....................  3d point values from array a2()  .................         
'
for k=0 to count
      p1 = a2(k)
      u = p1.u
     p4 = roty(p1,thi)
     p4 = rotx(p4,-0.75) 
     p4 = persp(p4,2)      
     pset(p4.x,p4.y),pal(u*11)
next k
'
ScreenCopy 1, 0
'
sleep 100
 next jt
 '
 '    repeat cycle (ng - ni)  imes
 '
 next ni
'
      return 0
'
end function
'
' ----------------------------------------------------------------
'
function  Generate(a1() as single,n as integer,m as integer,p as integer) as integer
'
'   Generate data for array a1() from function f3  .
'
static as integer i,j,k
static as single x,y,z,u
'
redim as single a1(0 to n-1,0 to m-1 , 0 to p-1)
'
'
'
'  .................... Generate 3d point values .................         
'
for k = 0 to p-1
      z = -1+2*k/(p-1)
 for j = 0 to m-1
      y = -1 + 2*j/(m-1)
  for i = 0 to n-1
       x = -1 + 2*i/(n-1)
       u =  f3(x,y,z)
       u = abs(u)
   a1(i,j,k) = u
  next i
 next j
next  k
'
      return 0
'
end function
'
' ----------------------------------------------------------------
'
function efilter(a1() as single , a2() as point , n as integer , m as integer , p as integer,l1 as single,l2 as single) as integer
'
'    Filter data  from array a1()  and selectively store in array a2()
'
'
static as integer i,j,k,count,cn
static as single x,y,z,u
static as point p1
'
'
u=l1
if l1 > l2 then
   l1=l2
   l2=u
end if
'
'
'  .................... Select 3d point values .................         
'
count = 0
for k = 0 to p-1
      z = -1+2*k/(p-1)
 for j = 0 to m-1
      y = -1 + 2*j/(m-1)
  for i = 0 to n-1
       x = -1 + 2*i/(n-1)
       u =  a1(i,j,k)
       
'
'   Set level where point becomes visible .
'
if (u >= l1 and u<=l2) then 
     count = count + 1
end if
'
  next i
 next j
next  k
'
if (count > 0 ) then
'
     redim as point  a2(count)
'
    cn =  0
for k = 0 to p-1
      z = -1+2*k/(p-1)
 for j = 0 to m-1
      y = -1 + 2*j/(m-1)
  for i = 0 to n-1
       x = -1 + 2*i/(n-1)
       u = a1(i,j,k)
'
'   Set level where point becomes visible .
'
if (u >= l1 and u<=l2) then   
      p1.x = x
      p1.y = y
      p1.z = z
      p1.u = u
  a2(cn) = p1
        cn = cn + 1
end if
'
  next i
 next j
next  k
'
end if
'
      return count
'
end function
'
' ----------------------------------------------------------------
'
function f3(x as single,y as single,z as single) as single
'
'        3d  sinc  function  
'
'
static as single u,v,w
'
    u=1
    if (abs(x)>0) then u= sin(3*x*Pi)/(x*3*Pi)
    v=1
    if (abs(y)>0) then v= sin(4*y*Pi)/(4*y*Pi)
    w=1
    if (abs(z)>0) then w= sin(5*z*Pi)/(5*z*Pi)
'
    return u*v*w
'
end function
'
' -----------------------------------------------------------
'



Luxan
Posts: 222
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: 3D Geometry , basics

Post by Luxan »

A very minor update , I'm calling this fn3d5a.bas .

I'm now using the log of the magnitude to represent the colors and have set the
levels higher ; this results in a more familiar 3d graph.

Just using points to represent the data appears to produce a sensible object ;
the use of perspective contributes to this.

Code: Select all


'
' --------------------------------------------------------------------
'
'               fn3d5a.bas
'
'    (c) Copyright 2016  sciwise@ihug.co.nz
'
'        Luxan
'
'        Edward.Q.Montague  [ alias ]
'
'
'                 3d function 
'
'
'         Screen pages are used  .
'
'
' -------------------------------------------------------------------
'
type point
         x as single
         y as single
         z as single
         u as single '  possible extension for special coord system
end type
'
' -------------------------------------------------------------------
'
declare function roty(q as point,angy as single) as point
declare function rotx(q as point,angx as single) as point
declare function rotz(q as point,angz as single) as point
declare function persp(q as point,d as single) as point
declare  function pally(pal() as long) as integer
declare function cube( p1a() as point, ByRef np as integer, edge() as integer,ByRef ne as integer) as integer
'
declare function f3(x as single,y as single,z as single) as single
declare function  Generate(a1() as single,n as integer,m as integer,p as integer) as integer
declare function efilter(a1() as single , a2() as point , n as integer , m as integer , p as integer,l1 as single,l2 as single) as integer
'
declare function dproj(a2() as point , count as integer , p1a() as point ,np as integer,edge() as integer,ne as integer, pal() as long , ng as integer , div as integer) as integer
'
' =====================================
'
const Pi = 4*atn(1)
dim as integer n,m,p,fg
'
 n = 128
m = 128
 p = 128
'
redim as single a1(0 to 4,0 to 4 , 0 to 4)
' ----------------------------------------------------------
'
dim pal(15) as long
fg = pally(pal() ) 
'
' -----------------------------------------------------------
'
dim as integer np,ne
dim as point p1a()
dim as integer edge()
dim as single l1,l2   '   levels
fg = cube(p1a() ,np ,edge() ,ne ) 
'
' -----------------------------------------------------------
'
screen 12,8,2
' image for working page #1 (visible page #0)
ScreenSet 1, 0 
Cls
window (-1.5,-1.5)-(1.5,1.5)
line (-1.4,-1.4)-(1.4,1.4),11,b
ScreenCopy 1, 0
'
redim a2(0 to 4) as point 
'
' --------------------------------------------------------------
'
fg = Generate(a1(),n ,m ,p ) 
'
l1=0.85
l2=1.0
fg= efilter(a1()  , a2()  , n  , m  , p ,l1 ,l2 ) 
'
fg= dproj(a2()  , fg , p1a()  ,np ,edge() ,ne , pal()  , 4  , 32 ) 
'
redim a2(4) as point
redim as single a1(0 to 4,0 to 4 , 0 to 4)
'
sleep
'
end
'
' ====================================
'
'
'    number of vertices
'
storeA:
data 8  
'
'    vertex data , easier to keep track of
'  data when we use multiple data statements.
'
store1: 
data  1,1,1
data  -1,1,1
data -1,-1,1
data 1,-1,1
data 1,1,-1
data  -1,1,-1
data  -1,-1,-1
data  1,-1,-1
'
'   Number of edges.
'
storeB:
data 12 
'
'  edge data 
'
store2:
data 1,2
data 1,4
data 1,5
data 2,3
data 2,6
data 3,4
data 3,7
data 4,8
data 5,6
data 5,8
data 6,7
data 7,8
'
'   Number of edges.
'
storeC:
data 9 
'
'  edge data  , back faces .
'
store3:
data   1 , 2
 data   1 , 4
 data   2 , 3
 data   3 , 4
 data   2 , 6
 data   3  , 7
data    6 , 7
 data   7 , 8
 data   4 , 8
'
' -------------------------------------------------------------------------------
'
function rotx(q as point,angx as single) as point
'
'              Rotate around x axis .
'
static as point p
'
             p.x = q.x
             p.y= q.y*cos(angx)-sin(angx)*q.z
             p.z= q.z*cos(angx)+sin(angx)*q.y
'
              return p
'
end function 
'
' -----------------------------------------------------------------------------
'
function roty(q as point,angy as single) as point
'
'              Rotate , a single point , around y axis .
'
static as point p
'
            p.x = sin(angy)*q.z + cos(angy)*q.x
            p.y = q.y
            p.z = cos(angy)*q.z -sin(angy)*q.x
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function rotz(q as point,angz as single) as point
'
'              Rotate around z axis .
'
static as point p
'
            p.x = sin(angz)*q.y + cos(angz)*q.x
            p.y = cos(angz)*q.y-sin(angz)*q.x
            p.z = q.z
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function persp(q as point,d as single) as point
'
'              3d  perspective .  
'
'     Applied to a single point .
'
'    Add 2 to the numerator when using any negative z value.
'
'     In this instance    -1 <= z  <= 1  ,  unit cube .
'
'    Therefore 2 is appropriate .
'
static as point p
'
     p.x = d*q.x/(q.z+3.5)
     p.y = d*q.y/(q.z+3.5)
     p.z = d
'
     return p
'
end function
'
' ---------------------------------------------------------------
'
function pally(pal() as long) as integer
'
'          Assign  palette information .
'

'
'     re defined 15 colour mode .
'
     pal(0)=5
     pal(1)=1
     pal(2)=9
     pal(3)=11
     pal(4)=2
     pal(5)=10
     pal(6)=14
     pal(7)=6
     pal(8)=4
     pal(9)=12
     pal(10)=13
     pal(11)=15
     pal(12)=0
     pal(13)=0
     pal(14)=0
     pal(15)=0
'
'
        return 0
'
end function
'
' --------------------------------------------------------------
'
function cube( p1a() as point, ByRef np as integer, edge() as integer,ByRef ne as integer) as integer
'
'         Read Cube coordinate points and edge sequence data . 
'
static as integer i
'
restore  storeA
read np
'
redim  p1a(1 to np) 
'
restore store1
for i =1 to np
   read p1a(i).x
   read p1a(i).y
   read p1a(i).z
next i
'
'
restore  storeC
read ne
'
redim  edge(1 to ne,0 to 1)
'
restore store3
for i = 1 to ne
   read edge(i,0)
   read edge(i,1)
next i
'
     return  0
'
end function
'
' ----------------------------------------------------------------
'
function dproj(a2() as point , count as integer , p1a() as point ,np as integer,edge() as integer,ne as integer, pal() as long , ng as integer , div as integer) as integer
'
' -----------------------------------------------------------
'
'   Draw cube frame and 
'
'   count  ,  number of points  in array a2()
'  
'
'   p1a()    point array for cube 
'   np number of points for cube array , this is usually a constant value .
'   edge()   edge array for cube
'   ne number of edge pairs for edge() array , this is usually a constant value.
' 
'  pal() , palette array ,  this is a constant length .
'  
'  ng , the number of cycles for the complete 360 deg rotation. 
'
' div , sets the  size of each  discrete rotation  ; 2*Pi/div .    
'
'
'   Generate values from 3d function for each rotation .
'
' -----------------------------------------------------------
'
static as single c1,d1,c2,d2,theta
dim as single x,y,z,u,maxa,thi
dim as single x1,y1,z1,u1,x2,y2,z2,u2
static as integer i,j,k,ni,jt,i1,j1,k1
static as point p1,p4,p2(1 to np)
'
c1=1.2
d1=0.4
c2=1.35
d2=1.3
'
theta = Pi/div
'
for ni=0 to ng
'
for jt = 0 to div
  cls
'
'  ....   Draw color palette ,,,,,,,
'
for i = 0 to 11
     y1=(d2-d1)*i/15+d1
     y2=(d2-d1)*(i+1)/15+d1
      line(c1,y1)-(c2,y2),pal(i),bf    
     line(c1,y1)-(c2,y2),11,b
next i
'  
'  .......... Draw outer cube .............
'
       thi = jt*theta
'       
   for k = 1 to np
     p2(k) = roty(p1a(k),thi)
     p2(k) = rotx(p2(k),-0.75) 
     p2(k)=persp(p2(k),2)
   next k     
'
for i1 = 1 to ne
      j1 = edge(i1,0)
     k1 = edge(i1,1)
   x1 = p2(j1).x
   y1 = p2(j1).y
'   z1 = p2(j1).z    
   x2 = p2(k1).x
   y2 = p2(k1).y
'   z2 = p2(k1).z    
line(x1,y1)-(x2,y2),14 
next i1 
'
'  ....................  3d point values from array a2()  .................         
'
for k=0 to count
      p1 = a2(k)
      u = p1.u
     p4 = roty(p1,thi)
     p4 = rotx(p4,-0.75) 
     p4 = persp(p4,2)      
     pset(p4.x,p4.y),pal(u*11)
next k
'
ScreenCopy 1, 0
'
sleep 100
 next jt
 '
 '    repeat cycle (ng - ni)  imes
 '
 next ni
'
      return 0
'
end function
'
' ----------------------------------------------------------------
'
function  Generate(a1() as single,n as integer,m as integer,p as integer) as integer
'
'   Generate data for array a1() from function f3  .
'
static as integer i,j,k
static as single x,y,z,u
'
redim as single a1(0 to n-1,0 to m-1 , 0 to p-1)
'
'
'
'  .................... Generate 3d point values .................         
'
for k = 0 to p-1
      z = -1+2*k/(p-1)
 for j = 0 to m-1
      y = -1 + 2*j/(m-1)
  for i = 0 to n-1
       x = -1 + 2*i/(n-1)
       u =  f3(x,y,z)
       u = abs(u)
   a1(i,j,k) = u
  next i
 next j
next  k
'
      return 0
'
end function
'
' ----------------------------------------------------------------
'
function efilter(a1() as single , a2() as point , n as integer , m as integer , p as integer,l1 as single,l2 as single) as integer
'
'    Filter data  from array a1()  and selectively store in array a2()
'
'
static as integer i,j,k,count,cn
static as single x,y,z,u
static as point p1
'
'
u=l1
if l1 > l2 then
   l1=l2
   l2=u
end if
'
'
'  .................... Select 3d point values .................         
'
count = 0
for k = 0 to p-1
      z = -1+2*k/(p-1)
 for j = 0 to m-1
      y = -1 + 2*j/(m-1)
  for i = 0 to n-1
       x = -1 + 2*i/(n-1)
       u =  a1(i,j,k)
       
'
'   Set level where point becomes visible .
'
if (u >= l1 and u<=l2) then 
     count = count + 1
end if
'
  next i
 next j
next  k
'
if (count > 0 ) then
'
     redim as point  a2(count)
'
    cn =  0
for k = 0 to p-1
      z = -1+2*k/(p-1)
 for j = 0 to m-1
      y = -1 + 2*j/(m-1)
  for i = 0 to n-1
       x = -1 + 2*i/(n-1)
       u = a1(i,j,k)
'
'   Set level where point becomes visible .
'
if (u >= l1 and u<=l2) then   
      p1.x = x
      p1.y = y
      p1.z = z
      p1.u = u
  a2(cn) = p1
        cn = cn + 1
end if
'
  next i
 next j
next  k
'
end if
'
      return count
'
end function
'
' ----------------------------------------------------------------
'
function f3(x as single,y as single,z as single) as single
'
'        3d  sinc  function  
'
'
static as single u,v,w
'
    u=1
    if (abs(x)>0) then u= sin(4*x*Pi)/(4*x*Pi)
    v=1
    if (abs(y)>0) then v= sin(2*y*Pi)/(2*y*Pi)
    w=1
    if (abs(z)>0) then w= sin(3*z*Pi)/(3*z*Pi)
'
    u=u*v*w
    u=log(abs(u)*12700000 + 1)/log(12700001)

    return u
'
end function
'
' -----------------------------------------------------------
'



Luxan
Posts: 222
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: 3D Geometry , basics

Post by Luxan »

Okay , now looking at the data points (x,y,z) & u , selected
via the efilter function , I notice that these are almost in the correct
format to use as a csv file ; within Paraview. However I don't have
enough experience with 3d viewers and their file formats to be able
to implement this .
Does anyone have any suggestions , or better yet sample
code to show how this might be done.
Luxan
Posts: 222
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: 3D Geometry , basics

Post by Luxan »

I may need to start a separate topic to cover everything that I'm sending
in this post .

Okay , continuing with the theme of 3d data points , I have generated a
3d array from a 3d fft of a 3d rectangular object.
Then I exchanged data values , within this array , to produce a centered
spectrum.
The selection levels were chosen to filter out the smaller magnitude values.
This I then graphed using the 3d functions developed previously , the whole
selection was rotated to highlight the features.

The 3d fft hasn't been thoroughly check , however the results appear to
be correct.
Any comment regards this is appreciated.
All of this code to be placed in the same directory.

There should be 6 portions of code .

All code :

(c) Copyright 2016 Luxan , sciwise@ihug.co.nz
1

Code: Select all


'
' --------------------------------------------------------------------
'
'               fn3d5ald.bi    ,   code as  an include file . 
'
'    (c) Copyright 2016  sciwise@ihug.co.nz
'
'        Luxan
'
'        Edward.Q.Montague  [ alias ]
'
'
'                 3d function
'
'
'         Screen pages are used  .
'
'
' -------------------------------------------------------------------
'
type point
         x as single
         y as single
         z as single
         u as single '  possible extension for special coord system
end type
'
' -------------------------------------------------------------------
'
declare function roty(q as point,angy as single) as point
declare function rotx(q as point,angx as single) as point
declare function rotz(q as point,angz as single) as point
declare function persp(q as point,d as single) as point
declare  function pally(pal() as long) as integer
declare function cube( p1a() as point, ByRef np as integer, edge() as integer,ByRef ne as integer) as integer
'
declare function f3(x as single,y as single,z as single) as single
declare function  Generate(b4() as single,n as integer,m as integer,p as integer) as integer
declare function normb(b4() as single ,  n as integer , m as integer , p as integer) as integer
declare function logb(b4() as single,n as integer,m as integer,p as integer) as integer
declare function efilter(b4() as single , b4c() as point , n as integer , m as integer , p as integer,l1 as single,l2 as single) as integer
'
declare function dproj(b4c() as point , count as integer , p1a() as point ,np as integer,edge() as integer,ne as integer, pal() as long , ng as integer , div as integer) as integer
'
' =====================================
'
const Pi = 4*atn(1)
'dim as integer fg
'
 'n = 8
'm = 8
' p = 8
'
redim as single b4(0 to 4,0 to 4 , 0 to 4)
' ----------------------------------------------------------
'
dim pal(15) as long
'
' -----------------------------------------------------------
'
dim as integer np,ne
dim as point p1a()
dim as integer edge()
dim as single l1,l2   '   levels
'
' ----------------------------------------------------------
'





2

Code: Select all

'
' ---------------------------------------------------
'
'     fn3d5l.bas
'
'
'               Load file for 3d graph
'
'
' ---------------------------------------------------
'
'    number of vertices
'
storeA:
data 8 
'
'    vertex data , easier to keep track of
'  data when we use multiple data statements.
'
store1:
data  1,1,1
data  -1,1,1
data -1,-1,1
data 1,-1,1
data 1,1,-1
data  -1,1,-1
data  -1,-1,-1
data  1,-1,-1
'
'   Number of edges.
'
storeB:
data 12
'
'  edge data
'
store2:
data 1,2
data 1,4
data 1,5
data 2,3
data 2,6
data 3,4
data 3,7
data 4,8
data 5,6
data 5,8
data 6,7
data 7,8
'
'   Number of edges.
'
storeC:
data 9
'
'  edge data  , back faces .
'
store3:
data   1 , 2
 data   1 , 4
 data   2 , 3
 data   3 , 4
 data   2 , 6
 data   3  , 7
data    6 , 7
 data   7 , 8
 data   4 , 8
'
' -------------------------------------------------------------------------------
'
function rotx(q as point,angx as single) as point
'
'              Rotate around x axis .
'
static as point p
'
             p.x = q.x
             p.y= q.y*cos(angx)-sin(angx)*q.z
             p.z= q.z*cos(angx)+sin(angx)*q.y
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function roty(q as point,angy as single) as point
'
'              Rotate , a single point , around y axis .
'
static as point p
'
            p.x = sin(angy)*q.z + cos(angy)*q.x
            p.y = q.y
            p.z = cos(angy)*q.z -sin(angy)*q.x
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function rotz(q as point,angz as single) as point
'
'              Rotate around z axis .
'
static as point p
'
            p.x = sin(angz)*q.y + cos(angz)*q.x
            p.y = cos(angz)*q.y-sin(angz)*q.x
            p.z = q.z
'
              return p
'
end function
'
' -----------------------------------------------------------------------------
'
function persp(q as point,d as single) as point
'
'              3d  perspective . 
'
'     Applied to a single point .
'
'    Add 2 to the numerator when using any negative z value.
'
'     In this instance    -1 <= z  <= 1  ,  unit cube .
'
'    Therefore 2 is appropriate .
'
static as point p
'
     p.x = d*q.x/(q.z+3.5)
     p.y = d*q.y/(q.z+3.5)
     p.z = d
'
     return p
'
end function
'
' ---------------------------------------------------------------
'
function pally(pal() as long) as integer
'
'          Assign  palette information .
'

'
'     re defined 15 colour mode .
'
     pal(0)=5
     pal(1)=1
     pal(2)=9
     pal(3)=11
     pal(4)=2
     pal(5)=10
     pal(6)=14
     pal(7)=6
     pal(8)=4
     pal(9)=12
     pal(10)=13
     pal(11)=15
     pal(12)=0
     pal(13)=0
     pal(14)=0
     pal(15)=0
'
'
        return 0
'
end function
'
' --------------------------------------------------------------
'
function cube( p1a() as point, ByRef np as integer, edge() as integer,ByRef ne as integer) as integer
'
'         Read Cube coordinate points and edge sequence data .
'
static as integer i
'
restore  storeA
read np
'
redim  p1a(1 to np)
'
restore store1
for i =1 to np
   read p1a(i).x
   read p1a(i).y
   read p1a(i).z
next i
'
'
restore  storeC
read ne
'
redim  edge(1 to ne,0 to 1)
'
restore store3
for i = 1 to ne
   read edge(i,0)
   read edge(i,1)
next i
'
     return  0
'
end function
'
' ----------------------------------------------------------------
'
function dproj(b4c() as point , count as integer , p1a() as point ,np as integer,edge() as integer,ne as integer, pal() as long , ng as integer , div as integer) as integer
'
' -----------------------------------------------------------
'
'   Draw cube frame and
'
'   count  ,  number of points  in array b3c()
' 
'
'   p1a()    point array for cube
'   np number of points for cube array , this is usually a constant value .
'   edge()   edge array for cube
'   ne number of edge pairs for edge() array , this is usually a constant value.
'
'  pal() , palette array ,  this is a constant length .
' 
'  ng , the number of cycles for the complete 360 deg rotation.
'
' div , sets the  size of each  discrete rotation  ; 2*Pi/div .   
'
'
'   Generate values from 3d function for each rotation .
'
' -----------------------------------------------------------
'
static as single c1,d1,c2,d2,theta
dim as single x,y,z,u,maxa,thi
dim as single x1,y1,z1,u1,x2,y2,z2,u2
static as integer i,j,k,ni,jt,i1,j1,k1
static as point p1,p4,p2(1 to np)
'
c1=1.2
d1=0.4
c2=1.35
d2=1.3
'
theta = Pi/div
'
for ni=0 to ng
'
for jt = 0 to div
  cls
'
'  ....   Draw color palette ,,,,,,,
'
for i = 0 to 11
     y1=(d2-d1)*i/15+d1
     y2=(d2-d1)*(i+1)/15+d1
      line(c1,y1)-(c2,y2),pal(i),bf   
     line(c1,y1)-(c2,y2),11,b
next i
' 
'  .......... Draw outer cube .............
'
       thi = jt*theta
'       
   for k = 1 to np
     p2(k) = roty(p1a(k),thi)
     p2(k) = rotx(p2(k),-0.75)
     p2(k)=persp(p2(k),2)
   next k     
'
for i1 = 1 to ne
      j1 = edge(i1,0)
     k1 = edge(i1,1)
   x1 = p2(j1).x
   y1 = p2(j1).y
'   z1 = p2(j1).z   
   x2 = p2(k1).x
   y2 = p2(k1).y
'   z2 = p2(k1).z   
line(x1,y1)-(x2,y2),14
next i1
'
'  ....................  3d point values from array b4c()  .................         
'
for k=0 to count
      p1 = b4c(k)
      u = p1.u
     p4 = roty(p1,thi)
     p4 = rotx(p4,-0.75)
     p4 = persp(p4,2)     
     pset(p4.x,p4.y),pal(u*11)
next k
'
ScreenCopy 1, 0
'
sleep 100
 next jt
 '
 '    repeat cycle (ng - ni)  imes
 '
 next ni
'
      return 0
'
end function
'
' ----------------------------------------------------------------
'
function  Generate(b4() as single,n as integer,m as integer,p as integer) as integer
'f
'   Generate data for array b4() from function f3  .
'
static as integer i,j,k
static as single x,y,z,u
'
redim as single b4(0 to n-1,0 to m-1 , 0 to p-1)
'
'
'
'  .................... Generate 3d point values .................         
'
for k = 0 to p-1
      z = -1+2*k/(p-1)
 for j = 0 to m-1
      y = -1 + 2*j/(m-1)
  for i = 0 to n-1
       x = -1 + 2*i/(n-1)
       u =  f3(x,y,z)
       u = abs(u)
   b4(i,j,k) = u
  next i
 next j
next  k
'
      return 0
'
end function
'
' ----------------------------------------------------------------
'
function efilter(b4() as single , b4c() as point , n as integer , m as integer , p as integer,l1 as single,l2 as single) as integer
'
'    Filter data  from array b4()  and selectively store in array b4c()
'
'
static as integer i,j,k,count,cn
static as single x,y,z,u
static as point p1
'
'
u=l1
if l1 > l2 then
   l1=l2
   l2=u
end if
'
'
'  .................... Select 3d point values .................         
'
count = 0
for k = 0 to p-1
      z = -1+2*k/(p-1)
 for j = 0 to m-1
      y = -1 + 2*j/(m-1)
  for i = 0 to n-1
       x = -1 + 2*i/(n-1)
       u =  b4(i,j,k)
       
'
'   Set level where point becomes visible .
'
if (u >= l1 and u<=l2) then
     count = count + 1
end if
'
  next i
 next j
next  k
'
if (count > 0 ) then
'
     redim as point  b4c(count)
'
    cn =  0
for k = 0 to p-1
      z = -1+2*k/(p-1)
 for j = 0 to m-1
      y = -1 + 2*j/(m-1)
  for i = 0 to n-1
       x = -1 + 2*i/(n-1)
       u = b4(i,j,k)
'
'   Set level where point becomes visible .
'
if (u >= l1 and u<=l2) then   
      p1.x = x
      p1.y = y
      p1.z = z
      p1.u = u
  b4c(cn) = p1
        cn = cn + 1
end if
'
  next i
 next j
next  k
'
end if
'
      return count
'
end function
'
' ----------------------------------------------------------------
'
function f3(x as single,y as single,z as single) as single
'
'        3d  sinc  function 
'
'
static as single u,v,w
'
    u=1
    if (abs(x)>0) then u= sin(4*x*Pi)/(4*x*Pi)
    v=1
    if (abs(y)>0) then v= sin(2*y*Pi)/(2*y*Pi)
    w=1
    if (abs(z)>0) then w= sin(3*z*Pi)/(3*z*Pi)
'
    u=u*v*w
    u=log(abs(u)*12700000 + 1)/log(12700001)

    return u
'
end function
'
' -----------------------------------------------------------
'
function normb(b4() as single ,  n as integer , m as integer , p as integer) as integer
'
'       Normalize all data values .
'
'
static as single max,fp
static as integer i,j,k
'
max=0
for k=0 to p-1
   for j=0 to m-1
     for i=0 to n-1
        fp=b4(i,j,k)
        fp=abs(fp)
    if ( fp > max ) then max = fp    
     next i
   next j
next k
'        
if (max =0 ) then max = 1
'
for k=0 to p-1
   for j=0 to m-1
     for i=0 to n-1
        fp=b4(i,j,k)
        fp=abs(fp)/max
        b4(i,j,k) = fp
     next i
   next j
next k
'  
      return 0
'
end function
'
' ----------------------------------------------------
'
function logb(b4() as single,n as integer,m as integer,p as integer) as integer
 '
 '        Natural Logarithm of all data in array b4()
 '
 '
 static as integer i,j,k
 static as single u
 '
 for k=0 to p-1
   for j=0 to m-1
     for i=0 to n-1
          u = b4(i,j,k)
           u=log(abs(u)*12700000 + 1)/log(12700001)   
           b4(i,j,k)=u
     next i
    next j
  next k   
 '
     return 0
 '
 end function

3

Code: Select all

'
' -------------------------------------------------------------------
'
'       3d  graph  of centered  spectrum from 3d fft  .
'
' -------------------------------------------------------------------
'
 #include "fourier.bas"
 #include "fourtst123.bi"
#include  "fn3d5ald.bi"

declare function  uxchg(bg3() as single,n as integer,m as integer,p as integer) as integer

'
' ---------------------------------------------------
'


fg = pally(pal() )
fg = cube(p1a() ,np ,edge() ,ne )
'
' -----------------------------------------------------------
'

n=8
m=8
p=8
'fg= ft1test(n )
'fg = ft2test(n , m ) 
'fg= ft3test(n ,m, p )
redim as single   b4(0 to n-1,0 to m-1,0 to p-1)
redim as point   b4c(4)
'fg=ft3gendat(b4()  ,n ,m , p ) 
'fg=normb(b4(),n,m,p)
'fg=ft3Prtdat(b4()  , n ,m , p ) 

'sleep
'end

screen 12,8,2
' image for working page #1 (visible page #0)
ScreenSet 1, 0
Cls
window (-1.5,-1.5)-(1.5,1.5)
line (-1.4,-1.4)-(1.4,1.4),11,b
ScreenCopy 1, 0
'
'redim b4c(0 to 4) as point
'
' --------------------------------------------------------------
'
n=128
m=128
p=128
redim as single   b4(0 to n-1,0 to m-1,0 to p-1)

'fg =  Generate(b4() ,n ,m ,p ) 

fg=ft3gendat(b4()  ,n ,m , p ) 
 fg=normb(b4(),n,m,p)
 fg=logb(b4() ,n ,m ,p ) 

 '  l1=0.85
 '  l2=1
'
l1=0.78
l2=1
'
fg=uxchg(b4()  , n ,m ,p ) 
fg= efilter(b4()  , b4c()  , n  , m  , p ,l1 ,l2 )
fg= dproj(b4c()  , fg , p1a()  ,np ,edge() ,ne , pal()  , 4  , 32 )
'
redim as single  b4(0 to 4,0 to 4 , 0 to 4)
redim b4c(4) as point

'
 sleep
 end
 '
 ' ================================
 '
#include "fourtst123x.bas"
#include "fn3d5l.bas"

function  uxchg(bg3() as single , n as integer,m as integer,p as integer) as integer
'
'
'                                         Exchange  bg3() values  into bg3()
'
'
static as integer i,j,k,i1,j1
'
'              [k,[i,j]]
'
                          redim b1(0 to m-1) as single
                          redim xb(0 to m-1) as single

                            for k=0 to p-1                               
                               for i=0 to n-1 
                                  for j=0 to m-1 
                                     b1(j)=bg3(i,j,k)
                                  next j
                                 
                                  for j=0 to (m/2)-1
                                      i1=j+(m/2)
                                    xb(j)=b1(i1)
                                    xb(i1)=b1(j) 
                               next j

                               for j=0 to m-1 
                                bg3(i,j,k)=xb(j)
                               next j
                           next i
                         next k

'exit function

'
'                                       [k,[j,i]]

                          redim b1(0 to n-1) as single
                          redim xb(0 to n-1) as single
                             
                           for k=0 to p-1  
                             for j=0 to m-1 

                               for i=0 to n-1 
                                  b1(i)=bg3(i,j,k)
                              next i
'
                               for i=0 to (n/2)-1
                                   j1=i+(n/2)
                               xb(i)=b1(j1)
                               xb(j1)=b1(i) 
                               next i

                                for i=0 to n-1 
                                    bg3(i,j,k)=xb(i)
                                 next i
                             next j
                           next k  
 
 '  exit function
 
 '
 '                                                 [j,[i,k]]
 '
                          redim b1(0 to p-1) as single
                          redim xb(0 to p-1) as single
                             
                        for j=0 to m-1 
                           for i = 0 to n-1 
'
                            for k = 0 to p-1  
                                  b1(k)=bg3(i,j,k)
                            next k
'
                            for k = 0 to (p/2)-1
                                    j1=k+(p/2)
                              xb(k)=b1(j1)
                             xb(j1)=b1(k) 
                           next  k

                                for k = 0 to p-1 
                                    bg3(i,j,k)=xb(k)
                                 next k
                             next i
                           next j  
 '
                   return 0                           
 '
end function                
'
' ---------------------------------------------------
'
4

Code: Select all

' ******************************************************************
' Fast Fourier Transform
' Modified from a Pascal program by Don Cross:
' http://groovit.disjunkt.com/analog/time-domain/fft.html
' ******************************************************************

' ------------------------------------------------------------------
' Error codes
' ------------------------------------------------------------------

#DEFINE FFTOk 0  
' No error

#DEFINE FFTPowErr -1 
' Number of samples is not a power of 2

#DEFINE FFTBoundErr -2 
' Arrays do not start at index 0

' ------------------------------------------------------------------
' Mathematical constant
' ------------------------------------------------------------------

#DEFINE TwoPi 6.283185307179586#
' 2*Pi

' ------------------------------------------------------------------
' Type definition
' ------------------------------------------------------------------

TYPE Complex
  X AS DOUBLE
  Y AS DOUBLE
END TYPE

' ------------------------------------------------------------------
' Global variable
' ------------------------------------------------------------------

COMMON SHARED ErrCode AS INTEGER
' Error code from the last function evaluation

' ******************************************************************

#DEFINE MaxPower 32

FUNCTION IsPowerOfTwo(BYVAL X AS INTEGER) AS INTEGER

  DIM AS INTEGER I, Y

  Y = 2
  FOR I = 1 TO MaxPower - 1
    IF X = Y THEN RETURN -1
    Y = Y SHL 1
  NEXT I
  RETURN 0
END FUNCTION

FUNCTION NumberOfBitsNeeded(BYVAL PowerOfTwo AS INTEGER) AS INTEGER
  
  DIM AS INTEGER I

  FOR I = 0 TO MaxPower
    IF (PowerOfTwo AND (1 SHL I)) <> 0 THEN RETURN I
  NEXT I
  RETURN 0
END FUNCTION

FUNCTION ReverseBits(BYVAL Index AS INTEGER, BYVAL NumBits AS INTEGER) AS INTEGER
  
  DIM AS INTEGER I, K, Rev

  Rev = 0
  I = Index
  FOR K = 0 TO NumBits - 1
    Rev = (Rev SHL 1) OR (I AND 1)
    I = I SHR 1
  NEXT K
  RETURN Rev
END FUNCTION

SUB FourierTransform(BYVAL AngleNumerator AS DOUBLE,  _
                           InArray()      AS Complex, _
                           OutArray()     AS Complex) 

  IF LBOUND(InArray) <> 0 OR LBOUND(OutArray) <> 0 THEN
    ErrCode = FFTBoundErr
    EXIT SUB
  END IF

  DIM AS INTEGER NumSamples
  
  NumSamples = UBOUND(InArray) + 1
  IF (NOT IsPowerOfTwo(NumSamples)) OR (NumSamples < 2) THEN
    ErrCode = FFTPowErr
    EXIT SUB
  END IF

  DIM AS INTEGER NumBits, I, J, K, N, BlockSize, BlockEnd
  DIM AS DOUBLE  Delta_angle, Delta_ar, Alpha, Beta, Tr, Ti, Ar, Ai

  ErrCode = FFTOk
  
  NumBits = NumberOfBitsNeeded(NumSamples)
  FOR I = 0 TO NumSamples - 1
    J = ReverseBits(I, NumBits)
    OutArray(J).X = InArray(I).X
    OutArray(J).Y = InArray(I).Y
  NEXT I

  BlockEnd = 1
  BlockSize = 2

  DO WHILE BlockSize <= NumSamples
    Delta_angle = AngleNumerator / BlockSize
    Alpha = SIN(0.5 * Delta_angle)
    Alpha = 2 * Alpha * Alpha
    Beta = SIN(Delta_angle)

    I = 0
    DO WHILE I < NumSamples
      Ar = 1  ' cos(0)
      Ai = 0  ' sin(0)

      J = I
      FOR N = 0 TO BlockEnd - 1
        K = J + BlockEnd
        Tr = Ar * OutArray(K).X - Ai * OutArray(K).Y
        Ti = Ar * OutArray(K).Y + Ai * OutArray(K).X
        OutArray(K).X = OutArray(J).X - Tr
        OutArray(K).Y = OutArray(J).Y - Ti
        OutArray(J).X = OutArray(J).X + Tr
        OutArray(J).Y = OutArray(J).Y + Ti
        Delta_ar = Alpha * Ar + Beta * Ai
        Ai = Ai - (Alpha * Ai - Beta * Ar)
        Ar = Ar - Delta_ar
        J = J + 1
      NEXT N

      I = I + BlockSize
    LOOP

    BlockEnd = BlockSize
    BlockSize = BlockSize SHL 1
  LOOP
END SUB

SUB FFT(InArray() AS Complex, OutArray() AS Complex)
' ------------------------------------------------------------------
' Calculates the Fast Fourier Transform of the array of
' complex numbers 'InArray' to produce the output complex
' numbers in 'OutArray'
' ------------------------------------------------------------------

  FourierTransform TwoPi, InArray(), OutArray()
END SUB

SUB IFFT(InArray() AS Complex, OutArray() AS Complex)
' ------------------------------------------------------------------
' Calculates the Inverse Fast Fourier Transform of the array of
' complex numbers represented by 'InArray' to produce the output
' complex numbers in 'OutArray'
' ------------------------------------------------------------------
  
  FourierTransform -TwoPi, InArray(), OutArray()

  IF ErrCode <> FFTOk THEN EXIT SUB
  
  DIM AS INTEGER NumSamples, I
  
  NumSamples = UBOUND(InArray) + 1
  
  ' Normalize the resulting time samples 
  FOR I = 0 TO NumSamples - 1
    OutArray(I).X = OutArray(I).X / NumSamples
    OutArray(I).Y = OutArray(I).Y / NumSamples
  NEXT I
END SUB

FUNCTION CalcFrequency(BYVAL FrequencyIndex AS INTEGER, _
                             InArray()      AS Complex) AS Complex
' ------------------------------------------------------------------
' This function calculates the complex frequency sample at a given
' index directly. Use this instead of 'FFT' when you only need one
' or two frequency samples, not the whole spectrum.
'
' It is also useful for calculating the Discrete Fourier Transform
' (DFT) of a number of data which is not an integer power of 2.
' For example, you could calculate the DFT of 100 points instead
' of rounding up to 128 and padding the extra 28 array slots with
' zeroes.
' ------------------------------------------------------------------

  IF LBOUND(InArray) <> 0 THEN
    ErrCode = FFTBoundErr
    EXIT FUNCTION
  END IF

  ErrCode = FFTOk
  
  DIM AS INTEGER NumSamples, K
  DIM AS DOUBLE  Cos1, Cos2, Cos3, Theta, Beta, Sin1, Sin2, Sin3
  DIM AS Complex Res

  NumSamples = UBOUND(InArray) + 1
  
  Res.X = 0
  Res.Y = 0

  Theta = TwoPi * FrequencyIndex / NumSamples
  Sin1  = SIN(- 2 * Theta)
  Sin2  = SIN(- Theta)
  Cos1  = COS(- 2 * Theta)
  Cos2  = COS(- Theta)
  Beta  = 2 * Cos2

  FOR K = 0 TO NumSamples - 1
    ' Update trig values
    Sin3 = Beta * Sin2 - Sin1
    Sin1 = Sin2
    Sin2 = Sin3

    Cos3 = Beta * Cos2 - Cos1
    Cos1 = Cos2
    Cos2 = Cos3

    Res.X = Res.X + InArray(K).X * Cos3 - InArray(K).Y * Sin3
    Res.Y = Res.Y + InArray(K).Y * Cos3 + InArray(K).X * Sin3
  NEXT K
  
  RETURN Res
END FUNCTION

5

Code: Select all

declare function fft1(a1() as Complex,b1() as Complex) as integer
declare function fft2(a2() as Complex  , b2() as Complex) as integer
declare function fft3(a3() as Complex,b3() as Complex) as integer
'
declare function ft1test(n as integer) as integer
declare function ft2test(n as integer,m as integer) as integer
declare function ft3test(n as integer,m as integer, p as integer) as integer
declare function ft3gendat(b4() as single ,n as integer,m as integer, p as integer) as integer
declare function ft3Prtdat(b4() as single , n as integer,m as integer, p as integer) as integer
'
'
'   Subroutine to call .
'
'
' FFT(InArray() AS Complex, OutArray() AS Complex)
'
dim as integer n,m,p,fg

6

Code: Select all

function fft1(a1() as Complex,b1() as Complex) as integer
 '
 '    One dimensional , complex fft
 '
 '     Input      array a1()
 '     Output   array b1()
 '
 '
           FFT(a1() , b1())
 '
        return 0
 '
 end function
'
' ------------------------------------------------------------
'
function fft2(a2() as Complex , b2() as Complex) as integer
'
'      Two dimensional , complex fft
'
'     Input      array a2()
'     Output   array b2()
'
'
static as integer i,j,n,m
n = ubound(a2)
m= ubound(a2,2)
'
'print " n , m  ";n;" , ";m
n=n+1
m=m+1
'
redim as complex a1(0 to n-1),b1(0 to n-1)
'
for j=0 to m-1
 for i=0 to n-1
     b2(i,j).X = a2(i,j).X
     b2(i,j).Y = a2(i,j).Y 
 next i
next j
'
'
for j=0 to m-1
'
 for i=0 to n-1
     a1(i).X = b2(i,j).X
     a1(i).Y = b2(i,j).Y
 next i
'
 FFT(a1() , b1())  
'
 for i=0 to n-1
    b2(i,j).X = b1(i).X
    b2(i,j).Y = b1(i).Y
 next i
'
next j
'
redim as complex a1(0 to m-1),b1(0 to m-1)
'
 for i = 0 to n - 1
 '
    for j = 0 to m - 1
      a1(j).X = b2(i,j).X
      a1(j).Y = b2(i,j).Y
    next j
'
    FFT(a1() , b1())  
' 
    for j=0 to m-1
       b2(i,j).X = b1(j).X
       b2(i,j).Y = b1(j).Y 
    next j
'
 next i
'
        return 0
'
end function
'
' ------------------------------------------------------------
'
function fft3(a3() as Complex,b3() as Complex) as integer
'
'      Three dimensional , complex fft
'
'     Input      array a3()
'     Output   array b3() '  alternately , just return a3()
'
'
'
static as integer i, j , k
static as integer n , m , p
n = ubound(a3)
m = ubound(a3,2)
p = ubound(a3,3)
'
redim as Complex a1(0 to 4) ,b1(0 to 4)
'
n  = n  + 1
m = m + 1
p  =  p + 1
'
'        Copy a3() to b3() 
'
for i=0 to n-1
  for j=0 to m-1
    for k=0 to p-1
       b3(i,j,k).X = a3(i,j,k).X
        b3(i,j,k).Y = a3(i,j,k).Y      
    next k
  next j
next i    
'
'     ---------- [ i , [ j , k ] ] ---------------
'
redim as Complex a1(0 to n-1) ,b1(0 to n-1)
'
 for k=0 to p-1
     for j=0 to m-1
   '
   for i=0 to n-1
     a1(i)=b3(i,j,k)
   next i
   '
   FFT(a1(),b1())
   '
  for i=0 to n-1
         b3(i,j,k)=b1(i)
  next i
  '
    next j
next k
'
' ...................            .[j,[i,k]]                                                  {k,i,j}
'
redim as Complex a1(0 to m-1) ,b1(0 to m-1)
'
for k=0 to p-1
  for i=0 to n-1
'
for j = 0 to m-1
    a1(j) = b3(i,j,k)
next j   
'
    FFT(a1() , b1())
'
for j = 0 to m-1
   b3(i,j,k) =  b1(j)
next  j
'
  next i
next k
'
' ...........
'
redim as Complex a1(0 to p-1) ,b1(0 to p-1)
'
for j=0 to m-1
for i=0 to n-1
'
for k = 0 to p-1
   a1(k) = b3(i,j,k)
next k
'
 FFT(a1() , b1())
'
for k = 0 to p-1
        b3(i,j,k)  =  b1(k)
next k
'
next i
next j
'
redim as Complex a1(0 to 4) ,b1(0 to 4)
'
        return 0
'
end function
'
' -----------------------------------------------------------
'
function ft1test(n as integer) as integer
'
'    Test fft1
'
'
static as integer i,fg
redim as Complex a1(0 to n-1),  b1(0 to n-1)
'
for i=0 to (n/2)-1
    a1(i).X = 1
next i
'
print
for i=0 to n - 1
  print " ";  a1(i).X;" , ";a1(i).Y
next i
'
print
'
   FFT(a1() , b1())
'fg= fft1(a1() ,b1() ) 
 '
print
for i=0 to n - 1
  print " ";  b1(i).X;" , ";b1(i).Y
next i
'
print
'
redim as Complex a1(0 to 3),  b1(0 to 3)
'
      return 0
'
end function
'
' ----------------------------------------------------------
'
function ft2test(n as integer,m as integer) as integer
'
'    Test  fft2
'
'
static as integer i,j,fg
static as single fp
redim as Complex a2(0 to n-1,0 to m-1),  b2(0 to n-1,0 to m-1)
'
for j = 0 to (m/2) - 1
 for i = 0 to (n/2) - 1
    a2(i,j).X = 1
 next i
next j
'
print
print " .........................................................."
print "   2d data , real "
print
for j = 0 to m - 1
 for i = 0 to n - 1
     print  "  ";a2(i,j).X ;
 next i
print
next j
print
'
fg = fft2(a2() ,b2() )
'
print " .........................................................."
print "  fft2d  , real   "
print
for j = 0 to m - 1
 for i = 0 to n - 1
      fp = b2(i,j).X
      fp=int(fp*1000)/1000
    print  "  ";fp ;
 next i
 print
next j
print
'
print " .........................................................."
print "     fft2d  ,  imag   "
print
for j = 0 to m - 1
print
 for i = 0 to n - 1
     fp = b2(i,j).Y
      fp=int(fp*1000)/1000
    print  "  ";fp ;
 next i
next j
print
'
redim as Complex a2(0 to 3,0 to 3),  b2(0 to 3,0 to 3)
'
     return 0
'
end function


function ft3test(n as integer,m as integer, p as integer) as integer
'
'    Test  fft3
'
'
static as integer i,j,k,fg
static as single fp,ep
redim as Complex a3(0 to n-1,0 to m-1,0 to p-1),  b3(0 to n-1,0 to m-1,0 to p-1)
'
for k=0 to (p/2)-1
for j = 0 to (m/2) - 1
 for i = 0 to (n/2) - 1
    a3(i,j,k).X = 1
 next i
next j
next k
'
'
'
print
print " ------------------------------------------"
print
print "   3d data , real "
print
for k=0 to p-  1
for j = 0 to m - 1
 for i = 0 to n - 1
     print  "  ";a3(i,j,k).X ;
 next i
print
next j
print
next k
print
'
fg = fft3(a3() ,b3() )
'
print
print" ..........................................................."
print "  fft3d  , real   "
print
for k=0 to p-1
for j = 0 to m - 1
 for i = 0 to n - 1
      fp = b3(i,j,k).X
      fp=int(fp*1000)/1000
    print  "  ";fp ;
 next i
  print
next j
 print
next k
print
'
print
print " .........................................................."
print "     fft3d  ,  imag   "
print
for k=0 to p-1
for j = 0 to m - 1
 for i = 0 to n - 1
      fp = b3(i,j,k).Y
      fp=int(fp*1000)/1000
    print  "  ";fp ;
 next i
  print
next j
 print
next k
print
'
'
print
print " ..........................................................."
print "     fft3d  ,  mag   "
print
for k=0 to p-1
for j = 0 to m - 1
 for i = 0 to n - 1
      ep = b3(i,j,k).X
      fp = b3(i,j,k).Y
      fp=ep*ep+fp*fp
      if (fp>0) then fp=sqr(fp)      
      fp=int(fp*1000)/1000
    print  "  ";fp ;
 next i
  print
next j
 print
next k
print
'
redim as Complex a3(0 to 3,0 to 3,0 to 3),  b3(0 to 3,0 to 3,0 to 3)
'
     return 0
'
end function
'
' ------------------------------------------------------------
'
function ft3gendat(b4() as single , n as integer,m as integer, p as integer) as integer
'
'    Test  fft3
'
'
static as integer i,j,k,fg
static as single fp,ep
redim as Complex a3(0 to n-1,0 to m-1,0 to p-1),  b3(0 to n-1,0 to m-1,0 to p-1)
redim as single b4(0 to n-1,0 to m-1,0 to p-1) 
'
for k=0 to (p/8)-1
for j = 0 to (m/3) - 1
 for i = 0 to (n/4) - 1
    a3(i,j,k).X = 1
 next i
next j
next k
'
fg = fft3(a3() ,b3() )
'
for k=0 to p-1
for j = 0 to m - 1
 for i = 0 to n - 1
      ep = b3(i,j,k).X
      fp = b3(i,j,k).Y
      fp=ep*ep+fp*fp
      if (fp>0) then fp=sqr(fp)      
      fp=int(fp*1000)/1000
      b4(i,j,k) = fp
 next i
next j
next k
'
redim as Complex a3(0 to 3,0 to 3,0 to 3)
'    ,  b3(0 to 3,0 to 3,0 to 3)
'
     return 0
'
end function
'
' ------------------------------------------------------------
'
function ft3Prtdat(b4() as single , n as integer,m as integer, p as integer) as integer
'
'       Print  data from array b3()
'
static as single fp
static as integer i,j,k
'
print
print " ..........................................................."
print "   ft3Prtdat ,   fft3d  ,  mag   "
print
for k=0 to p-1
for j = 0 to m - 1
 for i = 0 to n - 1
      fp = b4(i,j,k)
    print  "  ";fp ;
 next i
  print
next j
 print
next k
print
'
       return 0
'
end function
'
' --------------------------------------------------------------
'

Luxan
Posts: 222
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: 3D Geometry , basics

Post by Luxan »

The code fourier.bas , sent in previous post is part of fbmath module
and hence can't be copyrighted by me ; everything else however is
my own work.
Post Reply