General FreeBASIC programming questions.
Luxan
Posts: 236
Joined: Feb 18, 2009 12:47
Location: New Zealand

This topic is very messy, there's a routine for this in the skimage package of python,
however it's not easy to run or determine what he's doing.

I've written similar code, without the assistance of ChatGPT , previously in FreeBasic.

The inverse radon transform is [ was ] used for computerised axial tomography .

Therefore ChatGPT, and I, have produced a more basic version that many Basic programmers can
comprehend; documentation is on the agenda .

As typical, I'm using my own coding style .
I did have to wrestle with the code from ChatGPT , it didn't run as expected first time and manual
debugging was necessary .

The left side image is the original, the right hand side image has been reconstructed.

Anyway here's the code :

Code: Select all

``````
'
'             sciwiseg@gmail.com , 2024 ; lgpl , if applicable
'
'             74% help from ChatGPT
'
' Define the size of the image
const IMG_WIDTH =  256 '128
const IMG_HEIGHT = 256 '128
const N_ANGLES = 180
'
' -------------------------- declarations -----------------
'
Declare sub GenerateTestImage(image() as single)
Declare sub plot2img(image1() as single,image2() as single)
'
' ---------------------------------------------------------------
'
'       Main program
'
redim as single  image( IMG_WIDTH - 1,IMG_HEIGHT - 1)
redim as single radon(0 to 0, 0 to 0)
redim as single reconstructed( IMG_WIDTH - 1,IMG_HEIGHT - 1)
'
' ---------------------------------------------------------------
'
' Generate a test image
GenerateTestImage(image())
' Perform the inverse Radon transform
' Plot the original and reconstructed images .
plot2img(image(),  reconstructed())

sleep
end
'
' ===========================
'
'
' Function to generate a simple test image (a filled circle)
sub GenerateTestImage(image() as single)
dim as integer x, y
dim as single radius = 30
dim as single cx = IMG_WIDTH / 2
dim as single cy =IMG_HEIGHT / 2
dim as single dist, rd2
redim image( IMG_WIDTH - 1,IMG_HEIGHT - 1)
' Initialize the image with zeros
For y = 0 To IMG_HEIGHT - 1
For x = 0 To IMG_WIDTH - 1
image(x, y) = 0
Next
Next

for y = 0 to IMG_HEIGHT - 1
for x = 0 to IMG_WIDTH - 1
dist = ((x - cx) * (x - cx) + (y - cy) * (y - cy))
if dist <= rd2  then  image(x, y) = 1 end if
next x
next y
'
end sub
'
' ---------------------------------------------------------------
'
' Function to perform Radon Transform
dim as integer x, y, t
dim as single theta, rho
dim as integer r_index
dim as single pi = atn(1) * 4

dim as integer max_rho = sqr(IMG_WIDTH * IMG_WIDTH + IMG_HEIGHT * IMG_HEIGHT)
redim radon(0 to N_ANGLES - 1, 0 to 2 * max_rho - 1)

for t = 0 to N_ANGLES - 1
theta = t * pi / N_ANGLES
for y = 0 to IMG_HEIGHT - 1
for x = 0 to IMG_WIDTH - 1
if image(x, y) > 0 then
rho = x * cos(theta) + y * sin(theta)
r_index = max_rho + int(rho)
end if
next
next
next
'
end sub
'
' ---------------------------------------------------------------
'
dim as integer x, y, t
dim as single theta, rho
dim as integer r_index
dim as single pi = atn(1) * 4

dim as integer max_rho = sqr(ubound(image1, 1) * ubound(image1, 1) + ubound(image1, 2) * ubound(image1, 2))

redim image1( IMG_WIDTH - 1,IMG_HEIGHT - 1)
' Clear the image
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
image1(x, y) = 0
next
next

for t = 0 to N_ANGLES - 1
theta = t * pi / N_ANGLES
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
rho = x * cos(theta) + y * sin(theta)
r_index = max_rho + int(rho)
next
next
next
'
end sub
'
' ---------------------------------------------------------------
'
sub plot2img(image1() as single,image2() as single)
'
'       Plot 2 images, alongside each other .
'
dim as integer x, y
dim as single max_val1, min_val1
dim as single max_val2, min_val2
'
max_val1 = 0
min_val1=100
' Find the maximum  and minimum value in the  array .
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
if image1(x, y) > max_val1 then
max_val1 = image1(x, y)
end if
if image1(x, y) < min_val1 then
min_val1 = image1(x, y)
end if
next
next
'
'min_val1=0
max_val1=abs(max_val1-min_val1)
if max_val1=0 then max_val1=1 end if
'
' ......................................................................
'
'
max_val2 = 0
min_val2=100
' Find the maximum  and minimum value in the  array .
for y = 0 to ubound(image2, 2)
for x = 0 to ubound(image2, 1)
if image2(x, y) > max_val2 then
max_val2 = image2(x, y)
end if
if image2(x, y) < min_val2 then
min_val2 = image2(x, y)
end if
next
next
'
'min_val1=0
max_val2=abs(max_val2-min_val2)
if max_val2=0 then max_val2=1 end if
'
' ---------------------- plot images -----------------------
'
x=2*IMG_WIDTH
y=IMG_HEIGHT
screenres x, y, 32

'screenres 256, 128, 32
'
' Plot the radon array values as grayscale intensities
dim as integer intensity
'
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
intensity = int(255 * (image1(x, y)-min_val1) / max_val1)
pset (x, y), rgb(intensity, intensity, intensity)
next
next
'
' ......................................................................
'
for y = 0 to ubound(image2, 2)
for x = 0 to ubound(image2, 1)
intensity = int(255 * (image2(x, y)-min_val2) / max_val2)
pset (x+ IMG_WIDTH, y), rgb(intensity, intensity, intensity)
next
next
'
'
end sub
'
' ---------------------------------------------------------------
'

``````
dafhi
Posts: 1673
Joined: Jun 04, 2005 9:51

examination: N_ANGLES = 2

update 3

Code: Select all

``````dim as single  image() ' IMG_WIDTH - 1,IMG_HEIGHT - 1)
dim as single radon() '0 to 0, 0 to 0)
dim as single reconstructed() ' IMG_WIDTH - 1,IMG_HEIGHT - 1)
``````

Code: Select all

``````    redim radon(0 to N_ANGLES - 1, 0 to 1 * max_rho - 1)
'redim radon(0 to N_ANGLES - 1, 0 to 2 * max_rho - 1)
..
' compensated r_index calc
``````
full

Code: Select all

``````'
'             sciwiseg@gmail.com , 2024 ; lgpl , if applicable
'
'             74% help from ChatGPT
'
' Define the size of the image
const IMG_WIDTH =  256 '128
const IMG_HEIGHT = 256 '128
const N_ANGLES = 2'180
'
' -------------------------- declarations -----------------
'
Declare sub GenerateTestImage(image() as single)
Declare sub plotBoth(image1() as single,image2() as single)
'
' ---------------------------------------------------------------
'
'       Main program
'
dim as single  image() ' IMG_WIDTH - 1,IMG_HEIGHT - 1)
dim as single radon() '0 to 0, 0 to 0)
dim as single reconstructed() ' IMG_WIDTH - 1,IMG_HEIGHT - 1)
'
' ---------------------------------------------------------------
'
' Generate a test image
GenerateTestImage(image())
' Perform the inverse Radon transform
' Plot the original and reconstructed images .
plotBoth(image(),  reconstructed())

sleep
end
'
' ===========================
dim shared as integer x, y, t, r_index
dim shared as single _max, _min
dim shared as single dist, rd2
dim shared as single rho,  cost, sint, ysin
'    dim shared as single theta

const as single pi = atn(1) * 4
const as single radius = 30
'
'
' Function to generate a simple test image
sub GenerateTestImage(image() as single)
dim as single cx = IMG_WIDTH / 2
dim as single cy =IMG_HEIGHT / 2
redim image( IMG_WIDTH - 1,IMG_HEIGHT - 1)

For y = 0 To IMG_HEIGHT - 1
For x = 0 To IMG_WIDTH - 1
'' if statements are slow
'            if ((x - cx)^2 + (y - cy)^2) <= rd2  then  image(x, y) = 1
image(x, y) = -((x - cx)^2 + (y - cy)^2 <= rd2)
if rnd < .02 then image(x, y) = rnd '' test pattern by dafhi
Next
Next
'
end sub
'
' ---------------------------------------------------------------
'
' Function to perform Radon Transform

dim as integer max_rho = sqr(IMG_WIDTH ^ 2 + IMG_HEIGHT ^ 2)
redim radon(0 to N_ANGLES - 1, 0 to 1 * max_rho - 1)

for t = 0 to N_ANGLES - 1
cost = cos(t * pi / N_ANGLES)
sint = sin(t * pi / N_ANGLES)
for y = 0 to IMG_HEIGHT - 1
ysin = y * sint'(theta)
for x = 0 to IMG_WIDTH - 1
if image(x, y) > 0 then
rho = x * cost + ysin
r_index = int(rho + .5) '' + .5 eliminates artifact
'                    r_index = max_rho + int(rho + .5) '' + .5 eliminates artifact
end if
next
next
next
'
end sub
'
' ---------------------------------------------------------------
'

'    dim as integer max_rho = sqr( ubound(image1, 1) ^2 + ubound(image1, 2)^2 )
'    dim as integer max_rho = sqr(ubound(image1, 1) * ubound(image1, 1) + ubound(image1, 2) * ubound(image1, 2))

redim image1( IMG_WIDTH - 1,IMG_HEIGHT - 1)

for t = 0 to N_ANGLES - 1
cost = cos(t * pi / N_ANGLES)
sint = sin(t * pi / N_ANGLES)
for y = 0 to ubound(image1, 2)
ysin = y * sint'(theta)
for x = 0 to ubound(image1, 1)
rho = x * cost + ysin
r_index = int(rho)
'                r_index = max_rho + int(rho) '' + .5 doesn't seem necessary here, like it does in RadonTransform()
next
next
next
'
end sub
'
' ---------------------------------------------------------------
'
sub _plot_common( im() as single, x_offset as integer = 0)
_max = 0
_min=100
' Find the maximum  and minimum value in the  array .
for y = 0 to ubound(im, 2)
for x = 0 to ubound(im, 1)
if im(x, y) > _max then
_max = im(x, y)
elseif im(x, y) < _min then
_min = im(x, y)
end if
next
next
'
'_min=0
_max=abs(_max-_min)
if _max=0 then _max=1' end if
dim as integer intensity
_max = 255.999 / _max '' update 2

for y = 0 to ubound(im, 2)
for x = 0 to ubound(im, 1)
intensity = int( _max * ( im(x, y) - _min ))
pset (x + x_offset, y), rgb(intensity, intensity, intensity)
next
next
end sub

sub plotBoth(image1() as single,image2() as single)
'
'       Plot 2 images, alongside each other .
'
x=2*IMG_WIDTH
y=IMG_HEIGHT
screenres x, y, 32

_plot_common image1()
_plot_common image2(), IMG_WIDTH
end sub
'
' ---------------------------------------------------------------
'
``````

- update 2

speedup: moved sin() cos() to outermost loop, stored result in variables
-- RadonTransform() .. r_index = max_rho + int(rho + .5)
smaller code: dim shared same-name variables
rename: plot2img -> plotBoth (plot2image could mean several things)
speedup: _plot_common() .. intensity = int( _max * ( im(x, y) - _min ))

- update 1
smaller code & different pattern

Code: Select all

``````'
'             sciwiseg@gmail.com , 2024 ; lgpl , if applicable
'
'             74% help from ChatGPT
'
' Define the size of the image
const IMG_WIDTH =  256 '128
const IMG_HEIGHT = 256 '128
const N_ANGLES = 180
'
' -------------------------- declarations -----------------
'
Declare sub GenerateTestImage(image() as single)
Declare sub plot2img(image1() as single,image2() as single)
'
' ---------------------------------------------------------------
'
'       Main program
'
redim as single  image( IMG_WIDTH - 1,IMG_HEIGHT - 1)
redim as single radon(0 to 0, 0 to 0)
redim as single reconstructed( IMG_WIDTH - 1,IMG_HEIGHT - 1)
'
' ---------------------------------------------------------------
'
' Generate a test image
GenerateTestImage(image())
' Perform the inverse Radon transform
' Plot the original and reconstructed images .
plot2img(image(),  reconstructed())

sleep
end
'
' ===========================
'
'
' Function to generate a simple test image
sub GenerateTestImage(image() as single)
dim as integer x, y
dim as single radius = 30
dim as single cx = IMG_WIDTH / 2
dim as single cy =IMG_HEIGHT / 2
dim as single dist, rd2
redim image( IMG_WIDTH - 1,IMG_HEIGHT - 1)

For y = 0 To IMG_HEIGHT - 1
For x = 0 To IMG_WIDTH - 1
image(x, y) = 0'rnd - .9
if rnd < .02 then image(x, y) = rnd
Next
Next
/'
for y = 0 to IMG_HEIGHT - 1
for x = 0 to IMG_WIDTH - 1
dist = ((x - cx) * (x - cx) + (y - cy) * (y - cy))
if dist <= rd2  then  image(x, y) = 1 end if
next x
next y
'/
'
end sub
'
' ---------------------------------------------------------------
'
' Function to perform Radon Transform
dim as integer x, y, t
dim as single theta, rho
dim as integer r_index
dim as single pi = atn(1) * 4

dim as integer max_rho = sqr(IMG_WIDTH * IMG_WIDTH + IMG_HEIGHT * IMG_HEIGHT)
redim radon(0 to N_ANGLES - 1, 0 to 2 * max_rho - 1)

for t = 0 to N_ANGLES - 1
theta = t * pi / N_ANGLES
for y = 0 to IMG_HEIGHT - 1
for x = 0 to IMG_WIDTH - 1
if image(x, y) > 0 then
rho = x * cos(theta) + y * sin(theta)
r_index = max_rho + int(rho)
end if
next
next
next
'
end sub
'
' ---------------------------------------------------------------
'
dim as integer x, y, t
dim as single theta, rho
dim as integer r_index
dim as single pi = atn(1) * 4

dim as integer max_rho = sqr(ubound(image1, 1) * ubound(image1, 1) + ubound(image1, 2) * ubound(image1, 2))

redim image1( IMG_WIDTH - 1,IMG_HEIGHT - 1)
' Clear the image
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
image1(x, y) = 0
next
next

for t = 0 to N_ANGLES - 1
theta = t * pi / N_ANGLES
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
rho = x * cos(theta) + y * sin(theta)
r_index = max_rho + int(rho)
next
next
next
'
end sub
'
' ---------------------------------------------------------------
'
sub _plot_common( im() as single, x_offset as integer = 0)
dim as integer x, y
dim as single _max, _min
_max = 0
_min=100
' Find the maximum  and minimum value in the  array .
for y = 0 to ubound(im, 2)
for x = 0 to ubound(im, 1)
if im(x, y) > _max then
_max = im(x, y)
elseif im(x, y) < _min then
_min = im(x, y)
end if
next
next
'
'_min=0
_max=abs(_max-_min)
if _max=0 then _max=1' end if
dim as integer intensity

for y = 0 to ubound(im, 2)
for x = 0 to ubound(im, 1)
intensity = int( 255.999 * (im(x, y) - _min) / _max)
pset (x + x_offset, y), rgb(intensity, intensity, intensity)
next
next
end sub

sub plot2img(image1() as single,image2() as single)
'
'       Plot 2 images, alongside each other .
'
dim as integer x, y
x=2*IMG_WIDTH
y=IMG_HEIGHT
screenres x, y, 32

_plot_common image1()
_plot_common image2(), IMG_WIDTH
end sub
'
' ---------------------------------------------------------------
'
``````
Luxan
Posts: 236
Joined: Feb 18, 2009 12:47
Location: New Zealand

Interesting tweaks and improvements, I should time the improved code.

To generate a selection of test image data, from a function, I started to construct something like this :

Code: Select all

``````
'
'      rd0004.bas
'
'             sciwiseg@gmail.com , 2024 ; lgpl , if applicable
'
'             74% help from ChatGPT
'

' Define the size of the image
const IMG_WIDTH = 256 '128
const IMG_HEIGHT = 256 '128
const N_ANGLES = 180
'
' -------------------------- declarations -----------------
'
Declare sub GenerateTestImage(image() as single,choice as integer)
Declare sub plot2img(image1() as single,image2() as single)
'
Declare function f2(x as single,y as single,choice as integer) as single
'
' ---------------------------------------------------------------
'
'       Main program
'
redim as single  image( IMG_WIDTH - 1,IMG_HEIGHT - 1)
redim as single radon(0 to 0, 0 to 0)
redim as single reconstructed( IMG_WIDTH - 1,IMG_HEIGHT - 1)
'
' ---------------------------------------------------------------
'
Dim as integer choice
choice=1 ' filled circle
choice=2 ' rectangles

' Generate a test image
GenerateTestImage(image(),choice)
' Perform the inverse Radon transform
' Plot the original and reconstructed images .
plot2img(image(),  reconstructed())

sleep
end
'
' ===========================
'
'
' Function to generate a simple test image (a filled circle)
sub GenerateTestImage(image() as single,choice as integer)
dim as integer x, y
dim as single radius = 30
dim as single cx = IMG_WIDTH / 2
dim as single cy =IMG_HEIGHT / 2
dim as single dist, rd2
redim image( IMG_WIDTH - 1,IMG_HEIGHT - 1)
' Initialize the image with zeros
For y = 0 To IMG_HEIGHT - 1
For x = 0 To IMG_WIDTH - 1
image(x, y) = 0
Next
Next

for y = 0 to IMG_HEIGHT - 1
for x = 0 to IMG_WIDTH - 1

image(x, y) = f2(x,y,choice)
next x
next y

exit sub

for y = 0 to IMG_HEIGHT - 1
for x = 0 to IMG_WIDTH - 1
dist = ((x - cx) * (x - cx) + (y - cy) * (y - cy))
if dist <= rd2  then  image(x, y) = 1 end if
next x
next y
'
end sub
'
' ---------------------------------------------------------------
'
' Function to perform Radon Transform
dim as integer x, y, t
dim as single theta, rho
dim as integer r_index
dim as single pi = atn(1) * 4

dim as integer max_rho = sqr(IMG_WIDTH * IMG_WIDTH + IMG_HEIGHT * IMG_HEIGHT)
redim radon(0 to N_ANGLES - 1, 0 to 2 * max_rho - 1)

for t = 0 to N_ANGLES - 1
theta = t * pi / N_ANGLES
for y = 0 to IMG_HEIGHT - 1
for x = 0 to IMG_WIDTH - 1
if image(x, y) > 0 then
rho = x * cos(theta) + y * sin(theta)
r_index = max_rho + int(rho)
end if
next
next
next
'
end sub
'
' ---------------------------------------------------------------
'
dim as integer x, y, t
dim as single theta, rho
dim as integer r_index
dim as single pi = atn(1) * 4

dim as integer max_rho = sqr(ubound(image1, 1) * ubound(image1, 1) + ubound(image1, 2) * ubound(image1, 2))

redim image1( IMG_WIDTH - 1,IMG_HEIGHT - 1)
' Clear the image
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
image1(x, y) = 0
next
next

for t = 0 to N_ANGLES - 1
theta = t * pi / N_ANGLES
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
rho = x * cos(theta) + y * sin(theta)
r_index = max_rho + int(rho)
next
next
next
'
end sub
'
' ---------------------------------------------------------------
'
sub plot2img(image1() as single,image2() as single)
'
'       Plot 2 images, alongside each other .
'
dim as integer x, y
dim as single max_val1, min_val1
dim as single max_val2, min_val2
'
max_val1 = 0
min_val1=100
' Find the maximum  and minimum value in the  array .
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
if image1(x, y) > max_val1 then
max_val1 = image1(x, y)
end if
if image1(x, y) < min_val1 then
min_val1 = image1(x, y)
end if
next
next
'
'min_val1=0
max_val1=abs(max_val1-min_val1)
if max_val1=0 then max_val1=1 end if
'
' ......................................................................
'
'
max_val2 = 0
min_val2=100
' Find the maximum  and minimum value in the  array .
for y = 0 to ubound(image2, 2)
for x = 0 to ubound(image2, 1)
if image2(x, y) > max_val2 then
max_val2 = image2(x, y)
end if
if image2(x, y) < min_val2 then
min_val2 = image2(x, y)
end if
next
next
'
'min_val1=0
max_val2=abs(max_val2-min_val2)
if max_val2=0 then max_val2=1 end if
'
' ---------------------- plot images -----------------------
'
x=2*IMG_WIDTH
y=IMG_HEIGHT
screenres x, y, 32

'screenres 256, 128, 32
'
' Plot the radon array values as grayscale intensities
dim as integer intensity
'
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
intensity = int(255 * (image1(x, y)-min_val1) / max_val1)
pset (x, y), rgb(intensity, intensity, intensity)
next
next
'
' ......................................................................
'
for y = 0 to ubound(image2, 2)
for x = 0 to ubound(image2, 1)
intensity = int(255 * (image2(x, y)-min_val2) / max_val2)
pset (x+ IMG_WIDTH, y), rgb(intensity, intensity, intensity)
next
next
'
'
end sub
'
' ---------------------------------------------------------------
'
function f2(x as single,y as single,choice as integer) as single
'
'     Various , to produce test images .
'
dim as single z
'
z=0
'
if choice = 2 then
if  x>IMG_WIDTH/8 and x<IMG_WIDTH/2 and y>IMG_HEIGHT/4 and y<IMG_HEIGHT*3/4 then z=1 end if
if  x>=IMG_WIDTH/2 and x<5*IMG_WIDTH/6 and y>=IMG_HEIGHT*3/4 and y<IMG_HEIGHT-8 then z=1 end if
else

static as single radius = 30
static as single cx = IMG_WIDTH / 2
static as single cy =IMG_HEIGHT / 2
static as single dist, rd2
dist = ((x - cx) * (x - cx) + (y - cy) * (y - cy))
if dist <= rd2  then  z = 1 end if
end if

return z
'
end function
'
' ---------------------------------------------------------------
'

``````
dafhi
Posts: 1673
Joined: Jun 04, 2005 9:51

'radon' gives a pleasing bloom effect
..
want to draw attention to rendering artifact exposed when n_angles = 2
fix is to use int(rho + .5) in RadonTransform()
Luxan
Posts: 236
Joined: Feb 18, 2009 12:47
Location: New Zealand

Yes I noticed that also, I thought about the GLSL code when I observed this.

I believe this configuration involves a single source and a single detector .

Before we investigate the next stage; which is clearing up the bloom effect, here's
some code to test the not so obvious fft routines in the library fbmath .

Code: Select all

``````
'
'    fft_test.bas
'
'     sciwiseg@gmail.com
'     LGPL   2024
'
#INCLUDE "fbmath.bi"

declare sub center(u() as Complex, v() as Complex)

dim as integer n,i
n=8

dim as Complex s1(0 to n-1), p1(0 to n-1), q1(0 to n-1)
'
screen 11
'
for i=0 to (n/2) - 1
s1(i).X=1.0
next i
'
' FFT(InArray() AS Complex, OutArray() AS Complex)

FFT(s1(), q1())
'

for i=0 to n-1
' s1(i).X=i +1
next i

'center(q1() , p1() )
'
for i=0 to n - 1
print i;" , ";q1(i).x
next i
'
print
'
for i=0 to n - 1
print i;" , ";q1(i).y
next i
'
print "   "
print " fft "
for i=0 to n - 1
'   print i;" , ";p1(i).X;" , ";p1(i).Y
next i
'
sleep
end
'
' ======================================================================
'
sub center(u() as Complex, v() as Complex)
'
'   Center values, before and after fft .
'
static as integer i,j,n
static as Complex a,b
n = ubound(u,1) + 1
'
print " center, n = ",n
for i=0 to (n/2)-1
j=i+(n/2)
a=u(i)
b=u(j)
v(j)=a
v(i)=b
next i
'
end sub
'
' ----------------------------------------------------------------------
'

``````
dafhi
Posts: 1673
Joined: Jun 04, 2005 9:51

got that to work. i feel so accomplished.

watching youtu.be/f0sxjhGHRPo?t=213 .. too long

 - this is really cool. i think one could do a plain jane reconstruction with arc sampling. i'll attempt in next few weeks hopefully
dafhi
Posts: 1673
Joined: Jun 04, 2005 9:51

i've wondered how one might perform super resolution from multiple images

tomography may be the ticket. i'll save that for when i have more experience

atm i'm working out details for a parallax transform. will it work? i have no clue
Luxan
Posts: 236
Joined: Feb 18, 2009 12:47
Location: New Zealand

Good to hear you managed to setup fbmath and successfully test the fft example.

I noticed that video, before I began the post, it's fairly well presented.
Still, there's nothing like having code that actually does something; even better, fast understandable
code.

I'm part way through coding for what's known as a Ram_Lak filter, this is a basic high pass filter.

I'm curious about super resolution with this reconstruction method too; nowadays though super resolution might just be called up scaling, whereas previously this involved manipulation in the
Fourier domain.
Luxan
Posts: 236
Joined: Feb 18, 2009 12:47
Location: New Zealand

I used the NextPowerOfTwo routine, as fbmath appears to use radix 2 fft and ifft.
For Ram_Lak routine, I set the first value to 1, others set this to 0, this multiplies the dc term from the fft;
I also used Abs(freq)/n , the division is to avoid numerical overflow.

The image is flatter than previously, though there are now some overly large values at some edges; a consequence
of Ram Lak filter accentuating even the highest frequencies, some roll of is required.

You might find a better filter to go with the reconstruction .

Code: Select all

``````'
'      rd0004a.bas
'
'             sciwiseg@gmail.com , 2024 ; lgpl , if applicable
'
'             74% help from ChatGPT
'
'  Resize radon() row to closests power of 2, prior to fft .
'
#INCLUDE "fbmath.bi"
'
' Define the size of the image
const IMG_WIDTH = 256 '128
const IMG_HEIGHT = 256 '128
const N_ANGLES = 180
'
' -------------------------- declarations -----------------
'
Declare sub GenerateTestImage(image() as single,choice as integer)
Declare sub plot2img(image1() as single,image2() as single)
'
Declare function f2(x as single,y as single,choice as integer) as single
Declare sub Ram_Lak(FreqProjected() as complex)
Declare sub filterRows( radon() as single)
Declare Function NextPowerOfTwo(n As Integer) As Integer

'
' ---------------------------------------------------------------
'
'       Main program
'
redim as single  image( IMG_WIDTH - 1,IMG_HEIGHT - 1)
redim as single radon(0 to 0, 0 to 0)
redim as single reconstructed( IMG_WIDTH - 1,IMG_HEIGHT - 1)
'
' ---------------------------------------------------------------
'
Dim as integer choice
choice=1 ' filled circle
choice=2 ' rectangles

' Generate a test image
GenerateTestImage(image(),choice)
' Filter Rows , via fft , Ram_Lak , ifft
' Perform the inverse Radon transform
' Plot the original and reconstructed images .
plot2img(image(),  reconstructed())

sleep
end
'
' ===========================
'
'
' Function to generate a simple test image (a filled circle)
sub GenerateTestImage(image() as single,choice as integer)
dim as integer x, y
dim as single radius = 30
dim as single cx = IMG_WIDTH / 2
dim as single cy =IMG_HEIGHT / 2
dim as single dist, rd2
redim image( IMG_WIDTH - 1,IMG_HEIGHT - 1)
' Initialize the image with zeros
For y = 0 To IMG_HEIGHT - 1
For x = 0 To IMG_WIDTH - 1
image(x, y) = 0
Next
Next

for y = 0 to IMG_HEIGHT - 1
for x = 0 to IMG_WIDTH - 1

image(x, y) = f2(x,y,choice)
next x
next y

exit sub

for y = 0 to IMG_HEIGHT - 1
for x = 0 to IMG_WIDTH - 1
dist = ((x - cx) * (x - cx) + (y - cy) * (y - cy))
if dist <= rd2  then  image(x, y) = 1 end if
next x
next y
'
end sub
'
' ---------------------------------------------------------------
'
' Function to perform Radon Transform
dim as integer x, y, t
dim as single theta, rho
dim as integer r_index
'  dim as single pi = atn(1) * 4

dim as integer max_rho = sqr(IMG_WIDTH * IMG_WIDTH + IMG_HEIGHT * IMG_HEIGHT)
redim radon(0 to N_ANGLES - 1, 0 to 2 * max_rho - 1)

for t = 0 to N_ANGLES - 1
theta = t * pi / N_ANGLES
for y = 0 to IMG_HEIGHT - 1
for x = 0 to IMG_WIDTH - 1
if image(x, y) > 0 then
rho = x * cos(theta) + y * sin(theta)
r_index = max_rho + int(rho)
end if
next
next
next
'
end sub
'
' ---------------------------------------------------------------
'
dim as integer x, y, t
dim as single theta, rho
dim as integer r_index
'    dim as single pi = atn(1) * 4  ' predefined in fbmath .

dim as integer max_rho = sqr(ubound(image1, 1) * ubound(image1, 1) + ubound(image1, 2) * ubound(image1, 2))

redim image1( IMG_WIDTH - 1,IMG_HEIGHT - 1)
' Clear the image
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
image1(x, y) = 0
next
next

for t = 0 to N_ANGLES - 1
theta = t * pi / N_ANGLES
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
rho = x * cos(theta) + y * sin(theta)
r_index = max_rho + int(rho)
next
next
next
'
end sub
'
' ---------------------------------------------------------------
'
sub plot2img(image1() as single,image2() as single)
'
'       Plot 2 images, alongside each other .
'
dim as integer x, y
dim as single max_val1, min_val1
dim as single max_val2, min_val2
'
max_val1 = 0
min_val1=100
' Find the maximum  and minimum value in the  array .
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
if image1(x, y) > max_val1 then
max_val1 = image1(x, y)
end if
if image1(x, y) < min_val1 then
min_val1 = image1(x, y)
end if
next
next
'
'min_val1=0
max_val1=abs(max_val1-min_val1)
if max_val1=0 then max_val1=1 end if
'
' ......................................................................
'
'
max_val2 = 0
min_val2=100
' Find the maximum  and minimum value in the  array .
for y = 0 to ubound(image2, 2)
for x = 0 to ubound(image2, 1)
if image2(x, y) > max_val2 then
max_val2 = image2(x, y)
end if
if image2(x, y) < min_val2 then
min_val2 = image2(x, y)
end if
next
next
'
'min_val2=0
max_val2=abs(max_val2-min_val2)
if max_val2=0 then max_val2=1 end if
'
' ---------------------- plot images -----------------------
'
x=2*IMG_WIDTH
y=IMG_HEIGHT
screenres x, y, 32

'screenres 256, 128, 32
'
' Plot the radon array values as grayscale intensities
dim as integer intensity
'
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
intensity = int(255 * (image1(x, y)-min_val1) / max_val1)
pset (x, y), rgb(intensity, intensity, intensity)
next
next
'
' ......................................................................
'
for y = 0 to ubound(image2, 2)
for x = 0 to ubound(image2, 1)
intensity = int(255 * (image2(x, y)-min_val2) / max_val2)
pset (x+ IMG_WIDTH, y), rgb(intensity, intensity, intensity)
next
next
'
'
end sub
'
' ---------------------------------------------------------------
'
function f2(x as single,y as single,choice as integer) as single
'
'     Various , to produce test images .
'
dim as single z
'
z=0
'
if choice = 2 then
if  x>IMG_WIDTH/8 and x<IMG_WIDTH/2 and y>IMG_HEIGHT/4 and y<IMG_HEIGHT*3/4 then z=1 end if
if  x>=IMG_WIDTH/2 and x<5*IMG_WIDTH/6 and y>=IMG_HEIGHT*3/4 and y<IMG_HEIGHT-8 then z=1 end if
else

static as single radius = 30
static as single cx = IMG_WIDTH / 2
static as single cy =IMG_HEIGHT / 2
static as single dist, rd2
dist = ((x - cx) * (x - cx) + (y - cy) * (y - cy))
if dist <= rd2  then  z = 1 end if
end if

return z
'
end function
'
' ---------------------------------------------------------------
'
sub Ram_Lak(FreqProjected() as complex)
'
' Apply Ram-Lak filter in frequency domain
'  FreqProjected = fft(projected)
'
dim as integer i,n,freq
'
n=ubound(FreqProjected)+1

For i  = 0 To n - 1
freq  = i
If i=0 Then freq=1 end if
If i > n / 2 Then freq = freq - n end if
FreqProjected(i).x =FreqProjected(i).x*Abs(freq)/n
FreqProjected(i).y =FreqProjected(i).y*Abs(freq)/n
Next
'
end sub
'
' ----------------------------------------------------------------------
'
'
'
'
dim as integer  n,m,i,j,q
'
q=NextPowerOfTwo(n )
'
dim as complex rrow(0 to q-1),frrow(0 to q-1)
for i=0 to q-1
rrow(i).x=0
rrow(i).y=0
frrow(i).x=0
frrow(i).y=0
next i
'
for j=0 to m
for i=0 to n
next i
'
FFT(rrow(), frrow())
'
Ram_Lak(frrow() )
'
IFFT(frrow(), rrow())
'
for i=0 to n
next i
'
next j
'
redim as complex rrow(0 to 1),frrow(0 to 1)
'
end sub
'
' ---------------------------------------------------------------
'
Function NextPowerOfTwo(n As Integer) As Integer
Dim p As Integer = 1
While p < n
p *= 2
Wend
Return p
End Function
'
' ---------------------------------------------------------------
'

``````
dafhi
Posts: 1673
Joined: Jun 04, 2005 9:51

same as yours with fewer lines. you can choose not to use it - i have my own copy

[updated]
new GenerateTestImage, plot2img

Code: Select all

``````'
'      rd0004a.bas
'
'             sciwiseg@gmail.com , 2024 ; lgpl , if applicable
'
'             74% help from ChatGPT
'
'  Resize radon() row to closests power of 2, prior to fft .
'
#INCLUDE "../inc/fbmath.bi"
'
' Define the size of the image
const IMG_WIDTH = 256 '128
const IMG_HEIGHT = 256 '128
const N_ANGLES = 180
'
' -------------------------- declarations -----------------
'
Declare sub GenerateTestImage(image() as single,choice as integer)
Declare sub plot2img(image1() as single,image2() as single)
'
Declare function f2(x as single,y as single,choice as integer) as single
Declare sub Ram_Lak(FreqProjected() as complex)
Declare sub filterRows( radon() as single)
Declare Function NextPowerOfTwo(n As Integer) As Integer

'
' ---------------------------------------------------------------
'
'       Main program
'
redim as single  image( IMG_WIDTH - 1,IMG_HEIGHT - 1)
redim as single radon(0 to 0, 0 to 0)
redim as single reconstructed( IMG_WIDTH - 1,IMG_HEIGHT - 1)
'
' ---------------------------------------------------------------
'
Dim as integer choice
choice=1 ' filled circle
choice=2 ' rectangles

' Generate a test image
GenerateTestImage(image(),choice)
' Filter Rows , via fft , Ram_Lak , ifft
' Perform the inverse Radon transform
' Plot the original and reconstructed images .
plot2img(image(),  reconstructed())

sleep
end
'
' ===========================
dim shared as integer x, y, t, r_index, intensity '' _plot_common() intensity
'
'
' Function to generate a simple test image
sub GenerateTestImage(image() as single,choice as integer)

screenres 2*IMG_WIDTH, IMG_HEIGHT, 32 '' 2024 June 1

redim image( IMG_WIDTH - 1,IMG_HEIGHT - 1)

dim as any ptr im = imagecreate( img_width, img_height )

randomize 2 ' seed

for t = 1 to 10
circle im, (rnd * img_width, rnd * img_height), 5 + rnd*50, rgb(255,255,255)
next

for y = 0 to IMG_HEIGHT - 1
for x = 0 to IMG_WIDTH - 1
dim as ulong col = point(x,y,im)
image(x,y) = ((col shr 8) and 255)/255
'image(x, y) = f2(x,y,choice)
next
next

if imageinfo(im) = 0 then imagedestroy im

end sub
'
' ---------------------------------------------------------------
dim shared as single _max, _min
dim shared as single dist, rd2
dim shared as single rho,  cost, sint, ysin

' Function to perform Radon Transform

dim as integer max_rho = sqr(IMG_WIDTH ^ 2 + IMG_HEIGHT ^ 2)
redim radon(0 to N_ANGLES - 1, 0 to 1 * max_rho - 1)

for t = 0 to N_ANGLES - 1
cost = cos(t * pi / N_ANGLES)
sint = sin(t * pi / N_ANGLES)
for y = 0 to IMG_HEIGHT - 1
ysin = y * sint'(theta)
for x = 0 to IMG_WIDTH - 1
if image(x, y) > 0 then
rho = x * cost + ysin
r_index = int(rho)
'                    r_index = max_rho + int(rho)
end if
next
next
next
'
end sub
'
' ---------------------------------------------------------------
'

'max_rho = sqr( ubound(image1, 1) ^2 + ubound(image1, 2)^2 )
'    dim as integer max_rho = sqr(ubound(image1, 1) * ubound(image1, 1) + ubound(image1, 2) * ubound(image1, 2))

redim image1( IMG_WIDTH - 1, IMG_HEIGHT - 1)

for t = 0 to N_ANGLES - 1
cost = cos(t * pi / N_ANGLES)
sint = sin(t * pi / N_ANGLES)
for y = 0 to ubound(image1, 2)
ysin = y * sint'(theta)
for x = 0 to ubound(image1, 1)
rho = x * cost + ysin
r_index = int(rho)
'                r_index = max_rho + int(rho)
next
next
next
'
end sub
'
' ---------------------------------------------------------------
'
sub _plot_common( im() as single, x_offset as integer = 0)

_max = im( lbound(im,1), lbound(im,2) ) '' 2024 June 1
_min = _max

' Find the maximum  and minimum value in the  array .
for y = 0 to ubound(im, 2)
for x = 0 to ubound(im, 1)
if im(x, y) > _max then
_max = im(x, y)
elseif im(x, y) < _min then
_min = im(x, y)
end if
next
next
'
_max=255.999 / abs(_max-_min)

for y = 0 to ubound(im, 2)
for x = 0 to ubound(im, 1)
intensity = int( _max * ( im(x, y) - _min ))
pset (x + x_offset, y), rgb(intensity, intensity, intensity)
next
next

end sub

sub plot2img(image1() as single,image2() as single)
'
'       Plot 2 images, alongside each other .
'
_plot_common image1()
_plot_common image2(), IMG_WIDTH
end sub
'
' ---------------------------------------------------------------
'
function f2(x as single,y as single,choice as integer) as single
'
'     Various , to produce test images .
'
dim as single z
'
z=0
'
if choice = 2 then
if  x>IMG_WIDTH/8 and x<IMG_WIDTH/2 and y>IMG_HEIGHT/4 and y<IMG_HEIGHT*3/4 then z=1 end if
if  x>=IMG_WIDTH/2 and x<5*IMG_WIDTH/6 and y>=IMG_HEIGHT*3/4 and y<IMG_HEIGHT-8 then z=1 end if
else

static as single radius = 30
static as single cx = IMG_WIDTH / 2
static as single cy =IMG_HEIGHT / 2
static as single dist, rd2
dist = ((x - cx) * (x - cx) + (y - cy) * (y - cy))
if dist <= rd2  then  z = 1 end if
end if

return z
'
end function
'
' ---------------------------------------------------------------
'
sub Ram_Lak(FreqProjected() as complex)
'
' Apply Ram-Lak filter in frequency domain
'  FreqProjected = fft(projected)
'
dim as integer i,n,freq
'
n=ubound(FreqProjected)+1

For i  = 0 To n - 1
freq  = i
If i=0 Then freq=1 end if
If i > n / 2 Then freq = freq - n end if
FreqProjected(i).x =FreqProjected(i).x*Abs(freq)/n
FreqProjected(i).y =FreqProjected(i).y*Abs(freq)/n
Next
'
end sub
'
' ----------------------------------------------------------------------
'
'
'
'
dim as integer  n,m,i,j,q
'
q=NextPowerOfTwo(n )
'
dim as complex rrow(0 to q-1),frrow(0 to q-1)
for i=0 to q-1
rrow(i).x=0
rrow(i).y=0
frrow(i).x=0
frrow(i).y=0
next i
'
for j=0 to m
for i=0 to n
next i
'
FFT(rrow(), frrow())
'
Ram_Lak(frrow() )
'
IFFT(frrow(), rrow())
'
for i=0 to n
next i
'
next j
'
redim as complex rrow(0 to 1),frrow(0 to 1)
'
end sub
'
' ---------------------------------------------------------------
'
Function NextPowerOfTwo(n As Integer) As Integer
Dim p As Integer = 1
While p < n
p *= 2
Wend
Return p
End Function
'
' ---------------------------------------------------------------
'
``````
Luxan
Posts: 236
Joined: Feb 18, 2009 12:47
Location: New Zealand

So, we're puzzled by the artefacts in the reconstructed images .

If we consider the radon then inverse radon routines, upon the 2d image; as being some type of convolution,
then the 2d fft might be useful in giving us an indication of what this ' transfer function ' might be .
You'll need to look closely to observe a few dots in the log magnitude image of the transfer function .

This is just experimental and might be a departure from what's typically considered when finding a filter .

Here's a little extra code to do that.

Code: Select all

``````
'
'      fft2d_test.bas
'
'             sciwiseg@gmail.com , 2024 ; lgpl , if applicable
'
'             71% help from ChatGPT
'
#include "fbmath.bi"

' Define the size of the image
const IMG_WIDTH = 256 '256
const IMG_HEIGHT =256 ' 256
const N_ANGLES = 180

Declare sub GenerateTestImage(image() as single, choice as integer)
Declare sub Ram_Lak(FreqProjected() as complex)
Declare Function NextPowerOfTwo(n As Integer) As Integer
Declare sub Calculate2DFFT(image() as single, fft_result() as complex)
Declare sub Calculate2DTransferFunction(original_fft() as complex, reconstructed_fft() as complex, transfer_function() as complex)
Declare sub plot2img(image1() as single, image2() as single)

Declare sub plot2dSpect(image1() as complex)
'
' -----------------------------------------------------------------------
'
Dim as integer choice
choice = 1 ' Use 1 for filled circle, 2 for rectangles
'choice = 2 ' filled rectangles

' Generate a test image
Redim as single image(IMG_WIDTH - 1, IMG_HEIGHT - 1)
GenerateTestImage(image(), choice)

' Filter Rows, via fft, Ram_Lak, ifft

' Perform the inverse Radon transform
Redim as single reconstructed(IMG_WIDTH - 1, IMG_HEIGHT - 1)

' Calculate 2D FFTs
Dim as complex original_fft(IMG_WIDTH - 1, IMG_HEIGHT - 1)
Dim as complex reconstructed_fft(IMG_WIDTH - 1, IMG_HEIGHT - 1)
Calculate2DFFT(image(), original_fft())
Calculate2DFFT(reconstructed(), reconstructed_fft())

' Determine the 2D transfer function
Dim as complex transfer_function(IMG_WIDTH - 1, IMG_HEIGHT - 1)
Calculate2DTransferFunction(original_fft(), reconstructed_fft(), transfer_function())

' Plot the original and reconstructed images
plot2img(image(), reconstructed())
sleep
plot2dSpect(transfer_function())

Sleep
end
'
' ===========================
'
' Function to perform Radon Transform
dim as integer x, y, t
dim as single theta, rho
dim as integer r_index
dim as integer max_rho = sqr(IMG_WIDTH * IMG_WIDTH + IMG_HEIGHT * IMG_HEIGHT)
redim radon(0 to N_ANGLES - 1, 0 to 2 * max_rho - 1)

for t = 0 to N_ANGLES - 1
theta = t * pi / N_ANGLES
for y = 0 to IMG_HEIGHT - 1
for x = 0 to IMG_WIDTH - 1
if image(x, y) > 0 then
rho = x * cos(theta) + y * sin(theta)
r_index = max_rho + int(rho)
end if
next
next
next
end sub
'
' -------------------------------------------------------------------------------
'
' Function to perform Inverse Radon Transform
dim as integer x, y, t
dim as single theta, rho
dim as integer r_index
dim as integer max_rho = sqr(ubound(image1, 1) * ubound(image1, 1) + ubound(image1, 2) * ubound(image1, 2))

redim image1(IMG_WIDTH - 1, IMG_HEIGHT - 1)
' Clear the image
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
image1(x, y) = 0
next
next

for t = 0 to N_ANGLES - 1
theta = t * pi / N_ANGLES
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
rho = x * cos(theta) + y * sin(theta)
r_index = max_rho + int(rho)
next
next
next
end sub
'
' -------------------------------------------------------------------------------------
'
' Function to calculate 2D FFT
sub Calculate2DFFT(image() as single, fft_result() as complex)
dim as integer u, v, x, y
dim as integer N = IMG_WIDTH
dim as integer M = IMG_HEIGHT
redim fft_result(N - 1, M - 1)
dim as complex row_fft(N - 1)

' Calculate FFT for each row
for y = 0 to M - 1
for x = 0 to N - 1
row_fft(x).x = image(x, y)
row_fft(x).y = 0
next
FFT(row_fft(),row_fft() )
for x = 0 to N - 1
fft_result(x, y) = row_fft(x)
next
next

' Calculate FFT for each column of the row FFT result
for x = 0 to N - 1
for y = 0 to M - 1
row_fft(y) = fft_result(x, y)
next
FFT(row_fft(), row_fft())
for y = 0 to M - 1
fft_result(x, y) = row_fft(y)
next
next
end sub
'
' -------------------------------------------------------------
'
' Function to calculate 2D transfer function
sub Calculate2DTransferFunction(original_fft() as complex, reconstructed_fft() as complex, transfer_function() as complex)
dim as integer x, y
dim as complex w,z
dim as double u

for x = 0 to IMG_WIDTH - 1
for y = 0 to IMG_HEIGHT - 1
w= original_fft(x, y)
z = reconstructed_fft(x, y)
u=w.x^2 + w.y^2
if u <> 0 then
transfer_function(x, y).x = (z.x * w.x + z.y *w.y)/u
transfer_function(x, y).y = (z.y * w.x - z.x * w.y)/u
else
transfer_function(x, y).x = 0
transfer_function(x, y).y = 0
end if
next
next
end sub
'
' ---------------------------------------------------------------
'
sub plot2img(image1() as single,image2() as single)
'
'       Plot 2 images, alongside each other .
'
dim as integer x, y
dim as single max_val1, min_val1,z
dim as single max_val2, min_val2
'
max_val1 = 0
min_val1=100
' Find the maximum  and minimum value in the  array .
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
if image1(x, y) > max_val1 then
max_val1 = image1(x, y)
end if
if image1(x, y) < min_val1 then
min_val1 = image1(x, y)
end if
next
next
'
'min_val1=0
max_val1=abs(max_val1-min_val1)
if max_val1=0 then max_val1=1 end if
'
' ......................................................................
'
'
'    max_val2 = 0
min_val2=100
' Find the maximum  and minimum value in the  array .
for y = 0 to ubound(image2, 2)
for x = 0 to ubound(image2, 1)
if image2(x, y) > max_val2 then
max_val2 = image2(x, y)
end if
if image2(x, y) < min_val2 then
min_val2 = image2(x, y)
end if
next
next
'
min_val2=0
max_val2=abs(max_val2-min_val2)
if max_val2=0 then max_val2=1 end if
'
' ---------------------- plot images -----------------------
'
x=2*IMG_WIDTH
y=IMG_HEIGHT
screenres x, y, 32

'screenres 256, 128, 32
'
' Plot the radon array values as grayscale intensities
dim as integer intensity
'
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
intensity = int(255 * (image1(x, y)-min_val1) / max_val1)
pset (x, y), rgb(intensity, intensity, intensity)
next
next
'
' ......................................................................
'
for y = 0 to ubound(image2, 2)
for x = 0 to ubound(image2, 1)
z= (image2(x, y)-min_val2) / max_val2
if z<0.05 then z=0 end if
intensity = int(255 *z)
pset (x+ IMG_WIDTH, y), rgb(intensity, intensity, intensity)
next
next
'
'
end sub
'
' ---------------------------------------------------------------
'
sub plot2dSpect(image1() as complex)
'
'       Plot 2 images, alongside each other .
'
dim as integer x, y,w
dim as complex z
dim as single max_val1, min_val1
'    dim as single max_val2, min_val2
'
max_val1 = 0
min_val1=100
' Find the maximum  and minimum value in the  array .
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
z=image1(x,y)
w=z.x^2 + z.y^2
if w>0 then w=sqr(w) end if
if  w > max_val1 then
max_val1 = w
end if
if w < min_val1 then
min_val1 = w
end if
next
next
'
'min_val1=0
max_val1=abs(max_val1-min_val1)
if max_val1=0 then max_val1=1 end if
'
' ---------------------- plot image -----------------------
'
x=IMG_WIDTH
y=IMG_HEIGHT
screenres x, y, 32

'screenres 256, 128, 32
'
' Plot the magnitude of  array values as grayscale intensities
dim as integer intensity
dim as single Ag, Pd

Ag=10^16
Pd=log(Ag+1)
'
for y = 0 to ubound(image1, 2)
for x = 0 to ubound(image1, 1)
z=image1(x, y)
w=z.x*z.x + z.y*z.y
if w > 0 then w = sqr(w) end if
w=(w-min_val1) / max_val1
w=log(Ag*w+1)/Pd
intensity = int(255 *w )
pset (x, y), rgb(intensity, intensity, intensity)
next
next
'
' ......................................................................
'
end sub
'
' ---------------------------------------------------------------
'
function f2(x as single,y as single,choice as integer) as single
'
'     Various , to produce test images .
'
dim as single z
'
z=0
'
if choice = 2 then
if  x>IMG_WIDTH/8 and x<IMG_WIDTH/2 and y>IMG_HEIGHT/4 and y<IMG_HEIGHT*3/4 then z=1 end if
if  x>=IMG_WIDTH/2 and x<5*IMG_WIDTH/6 and y>=IMG_HEIGHT*3/4 and y<IMG_HEIGHT-8 then z=1 end if
else

static as single radius = 30
static as single cx = IMG_WIDTH / 2
static as single cy =IMG_HEIGHT / 2
static as single dist, rd2
dist = ((x - cx) * (x - cx) + (y - cy) * (y - cy))
if dist <= rd2  then  z = 1 end if
end if

return z
'
end function
'
' ---------------------------------------------------------------
'
' Function to generate a simple test image (a filled circle)
sub GenerateTestImage(image() as single,choice as integer)
dim as integer x, y
'   dim as single radius = 30
dim as single cx = IMG_WIDTH / 2
dim as single cy =IMG_HEIGHT / 2
dim as single dist, rd2
'
redim image( IMG_WIDTH - 1,IMG_HEIGHT - 1)
' Initialize the image with zeros
For y = 0 To IMG_HEIGHT - 1
For x = 0 To IMG_WIDTH - 1
image(x, y) = 0
Next
Next
'
for y = 0 to IMG_HEIGHT - 1
for x = 0 to IMG_WIDTH - 1
image(x, y) = f2(x,y,choice)
next x
next y
'
end sub
'
' ---------------------------------------------------------------
'

``````
dafhi
Posts: 1673
Joined: Jun 04, 2005 9:51

i'm a noob to FFTs and filters, plus i'm running windows now and can't get fbmath to work.

i like this whole tomography thing and have a new idea, so will proceed
Luxan
Posts: 236
Joined: Feb 18, 2009 12:47
Location: New Zealand

Don't worry about the fbmath library for freebasic, as I also have pure code for the fft and ifft routines.

I was intending to upload this , if there was an interest or need.
It'll take awhile for me to put it together .
Luxan
Posts: 236
Joined: Feb 18, 2009 12:47
Location: New Zealand

So many files, so disorganised ; I actually have multiple incarnations of the fft and ifft routines,
in different programming languages, just didn't want to sort through them all.

I asked copilot about copyright ownership of generated material, the U.S legal system, at the moment, is mostly
siding with the AI companies. If that isn't always going to be the instance, for generated code, then the onus
should be upon the AI companies to get their house in order. As a customer, or user, wouldn't you expect them
to do so ?

There's an online code search utility that you, and the AI firms, can use :
https://searchcode.com

Well, after all of that here's the code, that should be compatible with what I've already uploaded to this thread.

Code: Select all

``````
'
'    Microsoft copilot generated
'    sciwiseg@gmail.com  (c) 2024 LGPL , probably .
'
'
const M_PI =4*atn(1)

Type Complex
x As Double
y As Double
End Type
'
' ----------------------- Declarations --------------------
'
Declare Function CAdd(a As Complex, b As Complex) As Complex
Declare Function CExp(theta As Double) As Complex
Declare Function CMul(a As Complex, b As Complex) As Complex
Declare Function CSub(a As Complex, b As Complex) As Complex
Declare Sub FFT(data1() As Complex, ByVal n As Integer)
Declare Sub IFFT(data1() As Complex, ByVal n As Integer)
'
' --------------------------- main 1 -------------------------
'
' Example usage:
Dim As Integer size = 8 ' Must be a power of two
Dim As Complex signal(size - 1)

' Initialize the signal with some complex numbers
signal(0).x = 1.0 : signal(0).y = 0.0
signal(1).x = 1.0 : signal(1).y = 0.0
signal(2).x = 1.0 : signal(2).y = 0.0
signal(3).x = 1.0 : signal(3).y = 0.0

' Output the original array
For i As Integer = 0 To size - 1
Print "Signal["; i; "] = ("; signal(i).x; ", "; signal(i).y; ")"
Next
print

' Perform FFT
FFT(signal(), size)

' Output the transformed array
For i As Integer = 0 To size - 1
Print "Signal["; i; "] = ("; signal(i).x; ", "; signal(i).y; ")"
Next
print

' Perform IFFT
IFFT(signal(), size)

' Output the inverse transformed array
For i As Integer = 0 To size - 1
Print "Signal["; i; "] = ("; signal(i).x; ", "; signal(i).y; ")"
Next
Sleep

end
'
' ===========================
'

Function CAdd(a As Complex, b As Complex) As Complex
Dim As Complex result
result.x = a.x + b.x
result.y = a.y + b.y
Return result
End Function

Function CSub(a As Complex, b As Complex) As Complex
Dim As Complex result
result.x = a.x - b.x
result.y = a.y - b.y
Return result
End Function

Function CMul(a As Complex, b As Complex) As Complex
Dim As Complex result
result.x = a.x * b.x - a.y * b.y
result.y = a.x * b.y + a.y * b.x
Return result
End Function

Function CExp(theta As Double) As Complex
Dim As Complex result
result.x = Cos(theta)
result.y = Sin(theta)
Return result
End Function

Sub FFT(data1() As Complex, ByVal n As Integer)
If n <= 1 Then Return

Dim As Integer halfSize = n / 2
Dim As Integer k
Dim As Complex temp, e, o, twiddle

Dim As Complex even(halfSize - 1)
Dim As Complex odd(halfSize - 1)

For k = 0 To halfSize - 1
even(k) = data1(2 * k)
odd(k) = data1(2 * k + 1)
Next

FFT(even(), halfSize)
FFT(odd(), halfSize)

For k = 0 To halfSize - 1
twiddle = CExp(-2.0 * M_PI * k / n)
temp = CMul(twiddle, odd(k))

data1(k + halfSize) = CSub(even(k), temp)
Next
End Sub

Sub IFFT(data1() As Complex, ByVal n As Integer)
If n <= 1 Then Return

' Conjugate the complex numbers
For i As Integer = 0 To n - 1
data1(i).y = -data1(i).y
Next

' Call the FFT with conjugated complex numbers
FFT(data1(), n)

' Conjugate the complex numbers again
For i As Integer = 0 To n - 1
data1(i).y = -data1(i).y
Next

' Scale the numbers
For i As Integer = 0 To n - 1
data1(i).x = data1(i).x / n
data1(i).y = data1(i).y/ n
Next
End Sub

``````
dafhi
Posts: 1673
Joined: Jun 04, 2005 9:51

i copied fbmath's fourier.bas to a local folder.

from your previous post, this sub crashes unless i comment out certain lines

Code: Select all

``````sub Calculate2DFFT(image() as single, fft_result() as complex)
dim as integer u, v, x, y
dim as integer N = IMG_WIDTH
dim as integer M = IMG_HEIGHT

'    redim fft_result(N - 1, M - 1)
dim as complex row_fft(N - 1)

' Calculate FFT for each row
for y = 0 to M - 1
for x = 0 to N - 1
row_fft(x).x = image(x, y)
row_fft(x).y = 0
next
FFT(row_fft(),row_fft() )
for x = 0 to N - 1
'            fft_result(x, y) = row_fft(x)
next
next

' Calculate FFT for each column of the row FFT result
for x = 0 to N - 1
for y = 0 to M - 1
'            row_fft(y) = fft_result(x, y)
next
FFT(row_fft(), row_fft())
for y = 0 to M - 1
'            fft_result(x, y) = row_fft(y)
next
next
end sub
``````