## 3D Geometry , basics

User projects written in or related to FreeBASIC.
Luxan
Posts: 75
Joined: Feb 18, 2009 12:47
Location: New Zealand

### 3D Geometry , basics

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 systemend type'const Pi=4*atn(1)'dim as single x1,y1,z1,x2,y2,z2dim 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 pointdeclare function roty(q as point,angy as single) as pointdeclare function rotz(q as point,angz as single) as pointdeclare function tranx(q as point,movx as single) as pointdeclare function trany(q as point,movy as single) as pointdeclare function tranz(q as point,movz as single) as pointdeclare 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 store1for i=1 to 8   read p1(i).x   read p1(i).y   read p1(i).znext i'restore store2for i=1 to 12   read edge(i,0)   read edge(i,1)next i'' -----------------------------------------------------------------------------'screen 12window (-1.5,-1.5)-(1.5,1.5)line (-1.4,-1.4)-(1.4,1.4),11,b''clsTrall( p1() ,8,edge(), 32 )sleepend'' ===================================''     vertex data , easier to keep track of'  data when we use multiple data statements.'store1: data  1,1,1data -1,1,1data-1,-1,1data 1,-1,1data 1,1,-1data -1,1,-1data -1,-1,-1data 1,-1,-1''  edge data 'store2:data 1,2data 1,4data 1,5data 2,3data 2,6data 3,4data 3,7data 4,8data 5,6data 5,8data 6,7data 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,z2static as integer i,j,kstatic 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 jnext i'     return 0'end function'' --------------------------------------------------------------------------------'`
BasicCoder2
Posts: 3705
Joined: Jan 01, 2009 7:03
Location: Australia

### Re: 3D Geometry , basics

Have played around a bit myself although Dodicat has it down to a tee.
viewtopic.php?f=8&t=22842&hilit=3d+cube
viewtopic.php?f=7&t=23053&hilit=3d+cube
.
Luxan
Posts: 75
Joined: Feb 18, 2009 12:47
Location: New Zealand

### Re: 3D Geometry , basics

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

At the moment my programming momentum is in a slightly different direction ;

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: 75
Joined: Feb 18, 2009 12:47
Location: New Zealand

### Re: 3D Geometry , basics

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 systemend type'const Pi=4*atn(1)'dim as single x1,y1,z1,x2,y2,z2dim 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 pointdeclare function roty(q as point,angy as single) as pointdeclare function rotz(q as point,angz as single) as pointdeclare function tranx(q as point,movx as single) as pointdeclare function trany(q as point,movy as single) as pointdeclare function tranz(q as point,movz as single) as pointdeclare 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  storeAread np'dim as point p1(1 to np)'restore store1for i =1 to np   read p1(i).x   read p1(i).y   read p1(i).znext i''restore  storeBread ne'dim as integer edge(1 to ne,0 to 1)'restore store2for i = 1 to ne   read edge(i,0)   read edge(i,1)next i'' -----------------------------------------------------------------------------'screen 12window (-1.5,-1.5)-(1.5,1.5)line (-1.4,-1.4)-(1.4,1.4),11,b''clsk=Trall( p1() ,3,edge() , 32 ,np ,ne ) sleepend'' ===================================''    number of vertices'storeA:data 8  ''    vertex data , easier to keep track of'  data when we use multiple data statements.'store1: data  1,1,1data  -1,1,1data  -1,-1,1data 1,-1,1data  1,1,-1data  -1,1,-1data  -1,-1,-1data  1,-1,-1''   Number of edges.'storeB:data 12 ''  edge data 'store2:data 1,2data 1,4data 1,5data 2,3data 2,6data 3,4data 3,7data 4,8data 5,6data 5,8data 6,7data 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,z2static as integer i,j,kstatic 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 jnext i'     return 0'end function'' --------------------------------------------------------------------------------'`
Last edited by Luxan on Jan 26, 2016 6:37, edited 1 time in total.
BasicCoder2
Posts: 3705
Joined: Jan 01, 2009 7:03
Location: Australia

### Re: 3D Geometry , basics

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: 75
Joined: Feb 18, 2009 12:47
Location: New Zealand

### Re: 3D Geometry , basics

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: 75
Joined: Feb 18, 2009 12:47
Location: New Zealand

### Re: 3D Geometry , basics

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 systemend type'' -------------------------------------------------------------------'declare function roty(q as point,angy as single) as pointdeclare 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,p4dim as integer n,m,pdim as single x,y,z,u,maxa,thi,theta,div' n = 128m = 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,k1dim as single x1,y1,x2,y2'restore  storeAread np'dim as point p1a(1 to np),p2(1 to np)'restore store1for i =1 to np   read p1a(i).x   read p1a(i).y   read p1a(i).znext i''restore  storeBread ne'dim as integer edge(1 to ne,0 to 1)'restore store2for 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 Clswindow (-1.5,-1.5)-(1.5,1.5)line (-1.4,-1.4)-(1.4,1.4),11,bScreenCopy 1, 0'' -----------------------------------------------------------''   Draw cube frame and ''   Generate from 3d function for each rotation .'' -----------------------------------------------------------'dim as single c1,d1,c2,d2'c1=1.2d1=0.4c2=1.35d2=1.3'   div = 32theta = 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,bnext 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 jnext  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,1data  -1,1,1data -1,-1,1data 1,-1,1data 1,1,-1data  -1,1,-1data  -1,-1,-1data  1,-1,-1''   Number of edges.'storeB:data 12 ''  edge data 'store2:data 1,2data 1,4data 1,5data 2,3data 2,6data 3,4data 3,7data 4,8data 5,6data 5,8data 6,7data 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: 75
Joined: Feb 18, 2009 12:47
Location: New Zealand

### Re: 3D Geometry , basics

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 systemend type'' -------------------------------------------------------------------'declare function roty(q as point,angy as single) as pointdeclare function rotx(q as point,angx as single) as pointdeclare function rotz(q as point,angz as single) as pointdeclare function persp(q as point,d as single) as pointdeclare 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 integerdeclare  function pally(pal() as long) as integerdeclare 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 = 128m = 128 p = 128''dim as single a1(0 to n-1,0 to m-1 , 0 to p-1)' ----------------------------------------------------------'dim pal(15) as longfg = pally(pal() ) '' -----------------------------------------------------------'dim as integer np,nedim 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 Clswindow (-1.5,-1.5)-(1.5,1.5)line (-1.4,-1.4)-(1.4,1.4),11,bScreenCopy 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,1data  -1,1,1data -1,-1,1data 1,-1,1data 1,1,-1data  -1,1,-1data  -1,-1,-1data  1,-1,-1''   Number of edges.'storeB:data 12 ''  edge data 'store2:data 1,2data 1,4data 1,5data 2,3data 2,6data 3,4data 3,7data 4,8data 5,6data 5,8data 6,7data 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  , 7data    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  storeAread np'redim  p1a(1 to np) 'restore store1for i =1 to np   read p1a(i).x   read p1a(i).y   read p1a(i).znext i''restore  storeCread ne'redim  edge(1 to ne,0 to 1)'restore store3for 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,thetadim as single x,y,z,u,maxa,thidim as single x1,y1,z1,u1,x2,y2,z2,u2static as integer i,j,k,ni,jt,i1,j1,k1static as point p1,p4,p2(1 to np)'c1=1.2d1=0.4c2=1.35d2=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,bnext 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 jnext  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: 75
Joined: Feb 18, 2009 12:47
Location: New Zealand

### Re: 3D Geometry , basics

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 systemend type'' -------------------------------------------------------------------'declare function roty(q as point,angy as single) as pointdeclare function rotx(q as point,angx as single) as pointdeclare function rotz(q as point,angz as single) as pointdeclare function persp(q as point,d as single) as pointdeclare 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 integerdeclare  function pally(pal() as long) as integerdeclare 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 integerdeclare 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 = 128m = 128 p = 128''dim as single a1(0 to n-1,0 to m-1 , 0 to p-1)' ----------------------------------------------------------'dim pal(15) as longfg = pally(pal() ) '' -----------------------------------------------------------'dim as integer np,nedim 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 Clswindow (-1.5,-1.5)-(1.5,1.5)line (-1.4,-1.4)-(1.4,1.4),11,bScreenCopy 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 pointend'' ===================================='''    number of vertices'storeA:data 8  ''    vertex data , easier to keep track of'  data when we use multiple data statements.'store1: data  1,1,1data  -1,1,1data -1,-1,1data 1,-1,1data 1,1,-1data  -1,1,-1data  -1,-1,-1data  1,-1,-1''   Number of edges.'storeB:data 12 ''  edge data 'store2:data 1,2data 1,4data 1,5data 2,3data 2,6data 3,4data 3,7data 4,8data 5,6data 5,8data 6,7data 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  , 7data    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  storeAread np'redim  p1a(1 to np) 'restore store1for i =1 to np   read p1a(i).x   read p1a(i).y   read p1a(i).znext i''restore  storeCread ne'redim  edge(1 to ne,0 to 1)'restore store3for 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,thetadim as single x,y,z,u,maxa,thidim as single x1,y1,z1,u1,x2,y2,z2,u2static as integer i,j,k,ni,jt,i1,j1,k1static as point p1,p4,p2(1 to np)'c1=1.2d1=0.4c2=1.35d2=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,bnext 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 jnext  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,cnstatic as single x,y,z,ustatic as point p1'''  .................... Generate 3d point values .................         'count = 0for 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 = zend if'  next i next jnext  k'if (count > 0 ) then'     redim as point  a2(count)'    cn =  0for 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 + 1end if'  next i next jnext  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,thetadim as single x,y,z,u,maxa,thidim as single x1,y1,z1,u1,x2,y2,z2,u2static as integer i,j,k,ni,jt,i1,j1,k1static as point p1,p4,p2(1 to np)'c1=1.2d1=0.4c2=1.35d2=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,bnext 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: 75
Joined: Feb 18, 2009 12:47
Location: New Zealand

### Re: 3D Geometry , basics

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 systemend type'' -------------------------------------------------------------------'declare function roty(q as point,angy as single) as pointdeclare function rotx(q as point,angx as single) as pointdeclare function rotz(q as point,angz as single) as pointdeclare function persp(q as point,d as single) as pointdeclare  function pally(pal() as long) as integerdeclare 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 integerdeclare 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 = 128m = 128 p = 128''dim as single a1(0 to n-1,0 to m-1 , 0 to p-1)' ----------------------------------------------------------'dim pal(15) as longfg = pally(pal() ) '' -----------------------------------------------------------'dim as integer np,nedim as point p1a()dim as integer edge()dim as single l1,l2   '   levelsfg = cube(p1a() ,np ,edge() ,ne ) '' -----------------------------------------------------------'screen 12,8,2' image for working page #1 (visible page #0)ScreenSet 1, 0 Clswindow (-1.5,-1.5)-(1.5,1.5)line (-1.4,-1.4)-(1.4,1.4),11,bScreenCopy 1, 0'redim a2(0 to 4) as point l1=0.1l2=0.15fg = dfilter(a2()  , n  , m  , p , l1 , l2 ) fg= dproj(a2()  , fg , p1a()  ,np ,edge() ,ne , pal()  , 4  , 32 ) redim a2(4) as pointend'' ===================================='''    number of vertices'storeA:data 8  ''    vertex data , easier to keep track of'  data when we use multiple data statements.'store1: data  1,1,1data  -1,1,1data -1,-1,1data 1,-1,1data 1,1,-1data  -1,1,-1data  -1,-1,-1data  1,-1,-1''   Number of edges.'storeB:data 12 ''  edge data 'store2:data 1,2data 1,4data 1,5data 2,3data 2,6data 3,4data 3,7data 4,8data 5,6data 5,8data 6,7data 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  , 7data    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  storeAread np'redim  p1a(1 to np) 'restore store1for i =1 to np   read p1a(i).x   read p1a(i).y   read p1a(i).znext i''restore  storeCread ne'redim  edge(1 to ne,0 to 1)'restore store3for 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,cnstatic as single x,y,z,ustatic as point p1''u=l1if l1 > l2 then   l1=l2   l2=uend if'''  .................... Generate 3d point values .................         'count = 0for 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 = zend if'  next i next jnext  k'if (count > 0 ) then'     redim as point  a2(count)'    cn =  0for 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 + 1end if'  next i next jnext  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,thetadim as single x,y,z,u,maxa,thidim as single x1,y1,z1,u1,x2,y2,z2,u2static as integer i,j,k,ni,jt,i1,j1,k1static as point p1,p4,p2(1 to np)'c1=1.2d1=0.4c2=1.35d2=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,bnext 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: 75
Joined: Feb 18, 2009 12:47
Location: New Zealand

### Re: 3D Geometry , basics

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 systemend type'' -------------------------------------------------------------------'declare function roty(q as point,angy as single) as pointdeclare function rotx(q as point,angx as single) as pointdeclare function rotz(q as point,angz as single) as pointdeclare function persp(q as point,d as single) as pointdeclare  function pally(pal() as long) as integerdeclare 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 singledeclare function  Generate(a1() as single,n as integer,m as integer,p as integer) as integerdeclare 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 = 128m = 128 p = 128'redim as single a1(0 to 4,0 to 4 , 0 to 4)' ----------------------------------------------------------'dim pal(15) as longfg = pally(pal() ) '' -----------------------------------------------------------'dim as integer np,nedim as point p1a()dim as integer edge()dim as single l1,l2   '   levelsfg = cube(p1a() ,np ,edge() ,ne ) '' -----------------------------------------------------------'screen 12,8,2' image for working page #1 (visible page #0)ScreenSet 1, 0 Clswindow (-1.5,-1.5)-(1.5,1.5)line (-1.4,-1.4)-(1.4,1.4),11,bScreenCopy 1, 0'redim a2(0 to 4) as point '' --------------------------------------------------------------'fg = Generate(a1(),n ,m ,p ) 'l1=0.1l2=0.15fg= efilter(a1()  , a2()  , n  , m  , p ,l1 ,l2 ) 'fg= dproj(a2()  , fg , p1a()  ,np ,edge() ,ne , pal()  , 4  , 32 ) 'redim a2(4) as pointredim 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,1data  -1,1,1data -1,-1,1data 1,-1,1data 1,1,-1data  -1,1,-1data  -1,-1,-1data  1,-1,-1''   Number of edges.'storeB:data 12 ''  edge data 'store2:data 1,2data 1,4data 1,5data 2,3data 2,6data 3,4data 3,7data 4,8data 5,6data 5,8data 6,7data 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  , 7data    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  storeAread np'redim  p1a(1 to np) 'restore store1for i =1 to np   read p1a(i).x   read p1a(i).y   read p1a(i).znext i''restore  storeCread ne'redim  edge(1 to ne,0 to 1)'restore store3for 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,thetadim as single x,y,z,u,maxa,thidim as single x1,y1,z1,u1,x2,y2,z2,u2static as integer i,j,k,ni,jt,i1,j1,k1static as point p1,p4,p2(1 to np)'c1=1.2d1=0.4c2=1.35d2=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,bnext 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,kstatic 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 jnext  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,cnstatic as single x,y,z,ustatic as point p1''u=l1if l1 > l2 then   l1=l2   l2=uend if'''  .................... Select 3d point values .................         'count = 0for 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 + 1end if'  next i next jnext  k'if (count > 0 ) then'     redim as point  a2(count)'    cn =  0for 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 + 1end if'  next i next jnext  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: 75
Joined: Feb 18, 2009 12:47
Location: New Zealand

### Re: 3D Geometry , basics

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 systemend type'' -------------------------------------------------------------------'declare function roty(q as point,angy as single) as pointdeclare function rotx(q as point,angx as single) as pointdeclare function rotz(q as point,angz as single) as pointdeclare function persp(q as point,d as single) as pointdeclare  function pally(pal() as long) as integerdeclare 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 singledeclare function  Generate(a1() as single,n as integer,m as integer,p as integer) as integerdeclare 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 = 128m = 128 p = 128'redim as single a1(0 to 4,0 to 4 , 0 to 4)' ----------------------------------------------------------'dim pal(15) as longfg = pally(pal() ) '' -----------------------------------------------------------'dim as integer np,nedim as point p1a()dim as integer edge()dim as single l1,l2   '   levelsfg = cube(p1a() ,np ,edge() ,ne ) '' -----------------------------------------------------------'screen 12,8,2' image for working page #1 (visible page #0)ScreenSet 1, 0 Clswindow (-1.5,-1.5)-(1.5,1.5)line (-1.4,-1.4)-(1.4,1.4),11,bScreenCopy 1, 0'redim a2(0 to 4) as point '' --------------------------------------------------------------'fg = Generate(a1(),n ,m ,p ) 'l1=0.85l2=1.0fg= efilter(a1()  , a2()  , n  , m  , p ,l1 ,l2 ) 'fg= dproj(a2()  , fg , p1a()  ,np ,edge() ,ne , pal()  , 4  , 32 ) 'redim a2(4) as pointredim 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,1data  -1,1,1data -1,-1,1data 1,-1,1data 1,1,-1data  -1,1,-1data  -1,-1,-1data  1,-1,-1''   Number of edges.'storeB:data 12 ''  edge data 'store2:data 1,2data 1,4data 1,5data 2,3data 2,6data 3,4data 3,7data 4,8data 5,6data 5,8data 6,7data 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  , 7data    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  storeAread np'redim  p1a(1 to np) 'restore store1for i =1 to np   read p1a(i).x   read p1a(i).y   read p1a(i).znext i''restore  storeCread ne'redim  edge(1 to ne,0 to 1)'restore store3for 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,thetadim as single x,y,z,u,maxa,thidim as single x1,y1,z1,u1,x2,y2,z2,u2static as integer i,j,k,ni,jt,i1,j1,k1static as point p1,p4,p2(1 to np)'c1=1.2d1=0.4c2=1.35d2=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,bnext 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,kstatic 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 jnext  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,cnstatic as single x,y,z,ustatic as point p1''u=l1if l1 > l2 then   l1=l2   l2=uend if'''  .................... Select 3d point values .................         'count = 0for 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 + 1end if'  next i next jnext  k'if (count > 0 ) then'     redim as point  a2(count)'    cn =  0for 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 + 1end if'  next i next jnext  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: 75
Joined: Feb 18, 2009 12:47
Location: New Zealand

### Re: 3D Geometry , basics

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: 75
Joined: Feb 18, 2009 12:47
Location: New Zealand

### Re: 3D Geometry , basics

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 systemend type'' -------------------------------------------------------------------'declare function roty(q as point,angy as single) as pointdeclare function rotx(q as point,angx as single) as pointdeclare function rotz(q as point,angz as single) as pointdeclare function persp(q as point,d as single) as pointdeclare  function pally(pal() as long) as integerdeclare 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 singledeclare function  Generate(b4() as single,n as integer,m as integer,p as integer) as integerdeclare function normb(b4() as single ,  n as integer , m as integer , p as integer) as integerdeclare function logb(b4() as single,n as integer,m as integer,p as integer) as integerdeclare 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,nedim 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,1data  -1,1,1data -1,-1,1data 1,-1,1data 1,1,-1data  -1,1,-1data  -1,-1,-1data  1,-1,-1''   Number of edges.'storeB:data 12''  edge data'store2:data 1,2data 1,4data 1,5data 2,3data 2,6data 3,4data 3,7data 4,8data 5,6data 5,8data 6,7data 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  , 7data    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  storeAread np'redim  p1a(1 to np)'restore store1for i =1 to np   read p1a(i).x   read p1a(i).y   read p1a(i).znext i''restore  storeCread ne'redim  edge(1 to ne,0 to 1)'restore store3for 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,thetadim as single x,y,z,u,maxa,thidim as single x1,y1,z1,u1,x2,y2,z2,u2static as integer i,j,k,ni,jt,i1,j1,k1static as point p1,p4,p2(1 to np)'c1=1.2d1=0.4c2=1.35d2=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,bnext 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),14next 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,kstatic 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 jnext  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,cnstatic as single x,y,z,ustatic as point p1''u=l1if l1 > l2 then   l1=l2   l2=uend if'''  .................... Select 3d point values .................         'count = 0for 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 + 1end if'  next i next jnext  k'if (count > 0 ) then'     redim as point  b4c(count)'    cn =  0for 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 + 1end if'  next i next jnext  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,fpstatic as integer i,j,k'max=0for 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 jnext 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 jnext 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=8m=8p=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'endscreen 12,8,2' image for working page #1 (visible page #0)ScreenSet 1, 0Clswindow (-1.5,-1.5)-(1.5,1.5)line (-1.4,-1.4)-(1.4,1.4),11,bScreenCopy 1, 0''redim b4c(0 to 4) as point'' --------------------------------------------------------------'n=128m=128p=128redim 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.78l2=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 DOUBLEEND TYPE' ------------------------------------------------------------------' Global variable' ------------------------------------------------------------------COMMON SHARED ErrCode AS INTEGER' Error code from the last function evaluation' ******************************************************************#DEFINE MaxPower 32FUNCTION 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 0END FUNCTIONFUNCTION 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 0END FUNCTIONFUNCTION 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 RevEND FUNCTIONSUB 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  LOOPEND SUBSUB 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 SUBSUB 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 IEND SUBFUNCTION 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 ResEND FUNCTION`

5

Code: Select all

`declare function fft1(a1() as Complex,b1() as Complex) as integerdeclare function fft2(a2() as Complex  , b2() as Complex) as integerdeclare function fft3(a3() as Complex,b3() as Complex) as integer'declare function ft1test(n as integer) as integerdeclare function ft2test(n as integer,m as integer) as integerdeclare function ft3test(n as integer,m as integer, p as integer) as integerdeclare function ft3gendat(b4() as single ,n as integer,m as integer, p as integer) as integerdeclare 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,mn = ubound(a2)m= ubound(a2,2)''print " n , m  ";n;" , ";mn=n+1m=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 inext 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 , kstatic as integer n , m , pn = ubound(a3)m = ubound(a3,2)p = ubound(a3,3)'redim as Complex a1(0 to 4) ,b1(0 to 4)'n  = n  + 1m = m + 1p  =  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 jnext 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 jnext 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 inext k'' ...........'redim as Complex a1(0 to p-1) ,b1(0 to p-1)'for j=0 to m-1for 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 inext 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,fgredim as Complex a1(0 to n-1),  b1(0 to n-1)'for i=0 to (n/2)-1    a1(i).X = 1next i'printfor i=0 to n - 1  print " ";  a1(i).X;" , ";a1(i).Ynext i'print'   FFT(a1() , b1())'fg= fft1(a1() ,b1() )  'printfor i=0 to n - 1  print " ";  b1(i).X;" , ";b1(i).Ynext 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,fgstatic as single fpredim 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 inext j'printprint " .........................................................."print "   2d data , real "printfor j = 0 to m - 1 for i = 0 to n - 1     print  "  ";a2(i,j).X ; next iprintnext jprint'fg = fft2(a2() ,b2() )'print " .........................................................."print "  fft2d  , real   "printfor 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 printnext jprint'print " .........................................................."print "     fft2d  ,  imag   "printfor j = 0 to m - 1print for i = 0 to n - 1     fp = b2(i,j).Y      fp=int(fp*1000)/1000    print  "  ";fp ; next inext jprint'redim as Complex a2(0 to 3,0 to 3),  b2(0 to 3,0 to 3)'     return 0'end functionfunction ft3test(n as integer,m as integer, p as integer) as integer''    Test  fft3''static as integer i,j,k,fgstatic as single fp,epredim 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)-1for j = 0 to (m/2) - 1 for i = 0 to (n/2) - 1    a3(i,j,k).X = 1 next inext jnext k'''printprint " ------------------------------------------"printprint "   3d data , real "printfor k=0 to p-  1for j = 0 to m - 1 for i = 0 to n - 1     print  "  ";a3(i,j,k).X ; next iprintnext jprintnext kprint'fg = fft3(a3() ,b3() )'printprint" ..........................................................."print "  fft3d  , real   "printfor k=0 to p-1for 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  printnext j printnext kprint'printprint " .........................................................."print "     fft3d  ,  imag   "printfor k=0 to p-1for 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  printnext j printnext kprint''printprint " ..........................................................."print "     fft3d  ,  mag   "printfor k=0 to p-1for 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  printnext j printnext kprint'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,fgstatic as single fp,epredim 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)-1for j = 0 to (m/3) - 1 for i = 0 to (n/4) - 1    a3(i,j,k).X = 1 next inext jnext k'fg = fft3(a3() ,b3() )'for k=0 to p-1for 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 inext jnext 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 fpstatic as integer i,j,k'printprint " ..........................................................."print "   ft3Prtdat ,   fft3d  ,  mag   "printfor k=0 to p-1for j = 0 to m - 1 for i = 0 to n - 1      fp = b4(i,j,k)    print  "  ";fp ; next i  printnext j printnext kprint'       return 0'end function'' --------------------------------------------------------------'`
Luxan
Posts: 75
Joined: Feb 18, 2009 12:47
Location: New Zealand

### Re: 3D Geometry , basics

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.