Radon

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

Radon

Post by Luxan »

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 RadonTransform(image() as single, radon() as single)
Declare sub InverseRadonTransform(radon() as single, image1() 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 Radon transform
RadonTransform(image(), radon())
' Perform the inverse Radon transform
InverseRadonTransform(radon(), reconstructed())
' 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
    rd2=radius^2
 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
sub RadonTransform(image() as single, radon() as single)
    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)
                    radon(t, r_index) += image(x, y)
                end if
            next
        next
    next
'    
end sub
'
' ---------------------------------------------------------------
'
sub InverseRadonTransform(radon() as single, image1() as single)
    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)
                image1(x, y) += radon(t, r_index)
            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

Re: Radon

Post by dafhi »

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 RadonTransform(image() as single, radon() as single)
Declare sub InverseRadonTransform(radon() as single, image1() 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 Radon transform
RadonTransform(image(), radon())
' Perform the inverse Radon transform
InverseRadonTransform(radon(), reconstructed())
' 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 
    rd2=radius^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
sub RadonTransform(image() as single, radon() as single)
    
    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
                    radon(t, r_index) += image(x, y)
                end if
            next
        next
    next
'    
end sub
'
' ---------------------------------------------------------------
'
sub InverseRadonTransform(radon() as single, image1() as single)
    
'    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()
                image1(x, y) += radon(t, r_index)
            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
var types: const pi, radius
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 RadonTransform(image() as single, radon() as single)
Declare sub InverseRadonTransform(radon() as single, image1() 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 Radon transform
RadonTransform(image(), radon())
' Perform the inverse Radon transform
InverseRadonTransform(radon(), reconstructed())
' 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
    rd2=radius^2
 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
sub RadonTransform(image() as single, radon() as single)
    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)
                    radon(t, r_index) += image(x, y)
                end if
            next
        next
    next
'    
end sub
'
' ---------------------------------------------------------------
'
sub InverseRadonTransform(radon() as single, image1() as single)
    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)
                image1(x, y) += radon(t, r_index)
            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

Re: Radon

Post by Luxan »

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 RadonTransform(image() as single, radon() as single)
Declare sub InverseRadonTransform(radon() as single, image1() as single)
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 Radon transform
RadonTransform(image(), radon())
' Perform the inverse Radon transform
InverseRadonTransform(radon(), reconstructed())
' 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
    rd2=radius^2
 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
sub RadonTransform(image() as single, radon() as single)
    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)
                    radon(t, r_index) += image(x, y)
                end if
            next
        next
    next
'    
end sub
'
' ---------------------------------------------------------------
'
sub InverseRadonTransform(radon() as single, image1() as single)
    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)
                image1(x, y) += radon(t, r_index)
            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
    rd2=radius^2
              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

Re: Radon

Post by dafhi »

'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

Re: Radon

Post by Luxan »

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

Re: Radon

Post by dafhi »

got that to work. i feel so accomplished.

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

for my benefit https://www.youtube.com/watch?v=pZ7JlXagT0w

[edit] - 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

Re: Radon

Post by dafhi »

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

Re: Radon

Post by Luxan »

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

Re: Radon

Post by Luxan »

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 RadonTransform(image() as single, radon() as single)
Declare sub InverseRadonTransform(radon() as single, image1() as single)
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)
' Perform the Radon transform
RadonTransform(image(), radon())
' Filter Rows , via fft , Ram_Lak , ifft
filterRows( radon() )
' Perform the inverse Radon transform
InverseRadonTransform(radon(), reconstructed())
' 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
    rd2=radius^2
 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
sub RadonTransform(image() as single, radon() as single)
    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)
                    radon(t, r_index) += image(x, y)
                end if
            next
        next
    next
'    
end sub
'
' ---------------------------------------------------------------
'
sub InverseRadonTransform(radon() as single, image1() as single)
    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)
                image1(x, y) += radon(t, r_index)
            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
    rd2=radius^2
              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
'
' ----------------------------------------------------------------------
'
sub filterRows( radon() as single)
'
'
'
dim as integer  n,m,i,j,q
'
n=ubound(radon,1)
m=ubound(radon,2)
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
     rrow(i).x=radon(i,j)
 next i
 '
 FFT(rrow(), frrow())  
 '
Ram_Lak(frrow() )
 '
 IFFT(frrow(), rrow())     
'
  for i=0 to n
     radon(i,j)=rrow(i).x
 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

Re: Radon

Post by dafhi »

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 RadonTransform(image() as single, radon() as single)
Declare sub InverseRadonTransform(radon() as single, image1() as single)
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)
' Perform the Radon transform
RadonTransform(image(), radon())
' Filter Rows , via fft , Ram_Lak , ifft
filterRows( radon() )
' Perform the inverse Radon transform
InverseRadonTransform(radon(), reconstructed())
' 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
sub RadonTransform(image() as single, radon() as single)
    
    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)
                    radon(t, r_index) += image(x, y)
                end if
            next
        next
    next
'    
end sub
'
' ---------------------------------------------------------------
'
sub InverseRadonTransform(radon() as single, image1() as single)
    
    '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)
                image1(x, y) += radon(t, r_index)
            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
    rd2=radius^2
              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
'
' ----------------------------------------------------------------------
'
sub filterRows( radon() as single)
'
'
'
dim as integer  n,m,i,j,q
'
n=ubound(radon,1)
m=ubound(radon,2)
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
     rrow(i).x=radon(i,j)
 next i
 '
 FFT(rrow(), frrow())  
 '
Ram_Lak(frrow() )
 '
 IFFT(frrow(), rrow())     
'
  for i=0 to n
     radon(i,j)=rrow(i).x
 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

Re: Radon

Post by Luxan »

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 RadonTransform(image() as single, radon() as single)
Declare sub InverseRadonTransform(radon() as single, image1() as single)
'Declare sub filterRows(radon() as single)
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)

' Perform the Radon transform
Dim as single radon()
RadonTransform(image(), radon())

' Filter Rows, via fft, Ram_Lak, ifft
'filterRows(radon())

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

' 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
sub RadonTransform(image() as single, radon() as single)
    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)
                    radon(t, r_index) += image(x, y)
                end if
            next
        next
    next
end sub
'
' -------------------------------------------------------------------------------
'
' Function to perform Inverse Radon Transform
sub InverseRadonTransform(radon() as single, image1() as single)
    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)
                image1(x, y) += radon(t, r_index)
            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
    rd2=radius^2
              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
'    
'    rd2=radius^2
 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

Re: Radon

Post by dafhi »

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

Re: Radon

Post by Luxan »

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

Re: Radon

Post by Luxan »

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) = CAdd(even(k), temp)
        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

Re: Radon

Post by dafhi »

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
Post Reply