Radon

General FreeBASIC programming questions.
dafhi
Posts: 1682
Joined: Jun 04, 2005 9:51

Re: Radon

Post by dafhi »

my newest post is hidden, like a previous post

fxm and i commented on it. i'll reiterate.

my latest post is hidden also, it's on page 2 and you'll only see it if you preview reply
fxm wrote: Jun 04, 2024 11:50
dafhi wrote: Jun 04, 2024 11:30 unable to view your post until i reply Preview

Last hidden post from Luxan !
Luxan (04 Jun 2024, 03:50) wrote: This might be because of how arrays are being re dimension ed ; perhaps the array that's being passed to the
fft_result() needs to be declared as a redim array , before calling the routine.

Also check how the type complex is defined , often re , im are used ; instead of x , y .

I considered the fft routines from the fbmath source, I think I already have it on disk, probably worthwhile figuring
that out .
fxm
Moderator
Posts: 12219
Joined: Apr 22, 2009 12:46
Location: Paris suburbs, FRANCE

Re: Radon

Post by fxm »

dafhi wrote: Jun 15, 2024 7:55 my newest post is hidden, like a previous post

Indeed, second hidden post from other user in this same topic (visible only from preview reply):
dafhi (15 Jun 2024, 09:48) wrote: gen2 accumulator. forget ushort .. i'm using a ubyte

might be a possibility for a gen3

Code: Select all

#define dbl  as double '' reduced text
#define sng  as single

#define min( a, b)    iif( (a)<(b), (a), (b) )
#define max( a, b)    iif( (a)>(b), (a), (b) )
  
#macro sw( a, b, tmp )
  tmp = a: a = b: b = tmp
#endmacro

function modu( in as double, m as double = 1 ) as double
  return in - m * int( in / m )
end function

function round(in as double, places as ubyte = 2) as string
  dim as integer mul = 10 ^ places
  return str(csng( int(in * mul + .5) / mul) )
End Function

Union uARGB
  As Ulong        col
  Type: As UByte  B,G,R,A
  End Type
  declare constructor( as ulong = 0 )
  declare operator cast as string
  declare operator cast as ulong
  declare operator let( as ulong)
End Union

constructor uARGB( in as ulong )
  col = in
end constructor

operator uARGB.let( in as ulong )
  r = (in shr 16) and 255
  g = (in shr 8) and 255
  b = in and 255
end operator

operator uARGB.cast as string
  return str(r) + " " + str(g) + " " + str(b)
end operator

operator uARGB.cast as ulong
  return col
end operator

  dim shared as ulong ptr p32
  dim shared as long      wm, hm, pitchBy

sub _gfx_release( byref im as any ptr )
  if imageinfo(im) = 0 then imagedestroy im
  im = 0
end sub

  const tau = 8*atn(1)


'' image class
''
type t_image_info
  dim as long     w,h, bypp, pitch,rate
  dim as any ptr  pixels, im
  dim as string   driver_name
  declare sub     get_info( as any ptr = 0 )
end type

  sub _get_screen( byref i as t_image_info )
    _gfx_release i.im
    ScreenInfo i.w, i.h, , i.bypp, i.pitch, i.rate, i.driver_name
    i.pixels = screenptr
  end sub

  sub _get_image( byref i as t_image_info, im as any ptr )
    if im<>i.im then _gfx_release i.im '' 2024 June 3
    ImageInfo im, i.w, i.h, i.bypp, i.pitch, i.pixels
    i.im = im
  end sub
  
sub t_image_info.get_info( im as any ptr )
  if im = 0 then _get_screen this: exit sub
  _get_image this, im
end sub


  namespace line2d '' rasterizer  2024 June 15  by dafhi
  
/'

  example:  screenres 800,600,32
  
  myline.render_target 0
  myline.draw rnd*800,rnd*600,rnd*800,rnd*600, rgba(255,255,255,50)

'/

dim as t_image_info   im

sub render_target( _i as any ptr )
  im.get_info _i
  pitchBy = im.pitch \ 4 '' integer divide
  p32 = im.pixels
  wm = im.w - 1
  hm = im.h - 1
end sub
  
  function _no_intersect( byref x sng, byref y sng, byref x1 sng, byref y1 sng, slope sng, _max sng, _min sng = 0 ) as long
    static sng ptr xmin, xmax, ymin, ymax, xmin_y, xmax_y, ymin_x, ymax_x
    if x1 < x then
      xmin = @x1: xmax = @x
      xmin_y = @y1: xmax_y = @y
    else
      xmin = @x: xmax = @x1
      xmin_y = @y: xmax_y = @y1
    endif
    dim as long retval = *xmin >= _max orelse *xmax < _min
    if *xmin < _min then
      *xmin_y += (_min - *xmin) * slope
      *xmin = _min
    endif
    _max -= .0001
    if *xmax > _max then
      *xmax_y -= (*xmax - _max) * slope
      *xmax = _max
    endif
    return retval
  end function
  
  sub _make_unrenderable( byref x sng, byref x1 sng )
    x = 1
    x1 = 0
  end sub
  
    type t_returns
      sng ptr px,py
      sng     x,y, _x, _y, x1, slope
  end type
  
  dim as t_returns  ret

  dim sng   slope, s_temp, dx, dy
  
sub _absdxdy_sorted( byref x sng, byref y sng, byref x1 sng, byref y1 sng, byref px sng, byref py sng, _im_w sng, _im_h sng, border_test sng = 10 )
    if x > x1 then
  sw( x, x1, s_temp )
  sw( y, y1, s_temp )
  endif
  dx = x1 - x
  dy = y1 - y:  slope = dy / dx
  if _no_intersect( x,y, x1,y1, slope, _im_w-border_test, border_test ) then _make_unrenderable ret._x, ret.x1: exit sub
  if _no_intersect( y,x, y1,x1, dx/dy, _im_h-border_test, border_test ) then _make_unrenderable ret._x, ret.x1: exit sub
  ret._x = x - .4999
  ret._y = y - .4999
  ret.x1 = x1 - .4999
  ret.slope = slope
end sub

sub _common( byref x sng, byref y sng, byref x1 sng, byref y1 sng )
  const sng border_test = 0.0'113 '' .022 crashes .. 023 might :P
  if abs(y1 - y) > abs(x1 - x) then
    ret.px = @ret.y:  ret.py = @ret.x
    _absdxdy_sorted y, x, y1, x1, x,y, im.h, im.w, border_test
  else
    ret.px = @ret.x:  ret.py = @ret.y
    _absdxdy_sorted x, y, x1, y1, x,y, im.w, im.h, border_test
  endif
end sub

  #define line2d_loopS  _common x, y, x1, y1: ret.x=ret._x: ret.y=ret._y:  while ret.x < ret.x1: _
  if *ret.px > -0.5 andalso *ret.py > -0.5 andalso _
   *ret.px < (wm+.5) andalso *ret.py < (hm+.5) then
  
  #define line2d_loopE  endif: ret.y += ret.slope: ret.x += 1:  wend

sub draw( x sng, y sng, x1 sng, y1 sng, col as ulong)
    line2d_loopS
'  ppset *ret.px, *ret.py, col
  line2d_loopE
end sub

end namespace


  namespace roll_accum '' 2024 June 15  by dafhi
  
dim as single   _max = 255
dim as single   _min = 0

  type  accum_int as ubyte
  
const lenx8 = len(accum_int) * 8
const adj_ = max(1, lenx8 / 16)
  
const           cb_sum = lenx8 \ 2 + adj_

const           cb_c = lenx8 - cb_sum
const as single cb_v = 2^(cb_sum - adj_) \ 2^(cb_c) 'iif(0, 1, lenx8 - 2 * cb_c)

  const c_mag = 2 ^ cb_c, v_mag = 2 ^ cb_v, sum_mag = 2 ^ (lenx8 - cb_c)
  dim as single _scale_in, _scale_out
 
sub set_range( ma sng = 255, mi sng = 0 )
  _max = ma:  _min = mi
  _scale_in = (v_mag - 0) / ((ma - mi))
  _scale_out = 1 / _scale_in
end sub
  const       c_shl = lenx8 - cb_c
  
function fc( i as accum_int) as accum_int
  return i shr c_shl
end function
  
function readv( i as accum_int ) as single
  return _min + _scale_out * (i and (sum_mag - 1)) / fc(i)
end function


sub addv( byref des as accum_int, f as single )
  
    if des > 0 then:  static as ushort c, i
  
  f = (f - _min) * _scale_in + (des and (sum_mag - 1))
  c = fc(des)
    if c = (c_mag-1) then
  f /= (c + 1)
  i = c_mag - 1
  f *= i
    else
  i = c + 1
  endif
  
  des = (i shl c_shl) or (f and (sum_mag-1))
  '? tab(1); iif(rollover,"*", ""); round(v);" ";c+1;" "; round(s_new_val); " ";
  
    else '' des = 0  [uninitialized]
    
  set_range _max, _min
  des = (1 shl c_shl) or (f - _min) * _scale_in
  endif

end sub
  
end namespace ' -- roll_accum

  
  namespace radon_test '' 2024 June 15 - by dafhi

dim as t_image_info source

using roll_accum

type cRGB
  as accum_int        r,g,b
  declare sub         vals( as ubyte, as ubyte, as ubyte )
  declare sub         reset
  declare operator    cast as ulong
end type

sub cRGB.reset
  r = 0: g = 0: b = 0
end sub

sub cRGB.vals( rr as ubyte, gg as ubyte, bb as ubyte )
  reset
  addv r, rr
  addv g, gg
  addv b, bb
end sub

operator cRGB.cast as ulong
  return rgb( roll_accum.readv(r), roll_accum.readv(g), roll_accum.readv(b))
end operator

dim as cRGB           accRGB(), raycol
'dim as cRGB           backup_scan()

  sub _encode_initialize( im as any ptr )
source.get_info im
wm = source.w-1
hm = source.h-1
redim accRGB( wm, hm )
line2d.render_target imagecreate( source.w, source.h, rgb(0,0,0) )
'redim backup_scan( max( wm, hm ) )
end sub

  sub _cRGB_to_im
    for y as long = 0 to hm
  for x as long = 0 to wm
dim byref as cRGB q = accRGB( x,y )
p32[y*pitchBy + x] = rgb( readv(q.r), readv(q.g), readv(q.b) )
next
next
end sub

  sub shoot_xray( x sng, y sng, x1 sng, y1 sng ): using line2d
raycol.reset
dim as uargb ptr a = source.pixels
line2d_loopS
dim as long i = *ret.px + clng(*ret.py) * pitchBy
roll_accum.addv raycol.r, a[ i ].r
roll_accum.addv raycol.g, a[ i ].g
roll_accum.addv raycol.b, a[ i ].b
line2d_loopE
end sub
  
  sub _line_accum( x sng, y sng, x1 sng, y1 sng, c as ulong ): using line2d
  #if 0
line line2d.im.im, (x,y)-(x1,y1),c
  #else
  line2d_loopS
dim as uargb uar = c
dim byref as cRGB q = accRGB( *ret.px, *ret.py )
roll_accum.addv q.r, uar.r
roll_accum.addv q.g, uar.g
roll_accum.addv q.b, uar.b
line2d_loopE
#endif
end sub

  sub process( im as any ptr )
_encode_initialize im
dim as long c_pixels = source.w * source.h
dim as single diagonal = sqr( source.w^2 + source.h^2)
for i_line as long = 1 to 4999
  dim as single si = rnd*c_pixels
  dim as single y = si / source.w
  dim as single x = si - int(y) * source.w
  dim as single a = rnd * tau
  dim as single slen = diagonal'sqr(diagonal) + rnd * diagonal / 3
  dim as single x1 = x + slen * cos(a): x -= slen * cos(a)
  dim as single y1 = y + slen * sin(a): y -= slen * sin(a)
  shoot_xray x,y,x1,y1
  _line_accum x,y,x1,y1, raycol
next
_cRGB_to_im
put (source.w,0), line2d.im.im, pset
end sub

end namespace ' -- radon_test


screenres 800,600, 32

dim as t_image_info   buf

buf.get_info

dim as any ptr im = imagecreate( buf.w/2, buf.h, iif(1, rgb(255,255,255), 0) )

randomize 7

for i as long = 0 to 9
  circle im, (rnd * buf.w/2, rnd * buf.h), 4 + rnd * 50, rgb(0,rnd*255,0),,,,f
'  circle im, (rnd * buf.w/2, rnd * buf.h), 4 + rnd * 50, rnd * culng(-1),,,,f
next

put (0,0), im, pset

radon_test.process im

sleep


Note:
- The first post that was hidden now reappears at the top of page 2 (did some admin do something, or consequence of saving forum, or another cause?).
- This second hidden post also corresponds to a page change!
- Let is wait for the next page break to page 4, to see if this hidden post reappears at the top of page 3!
Last edited by fxm on Jun 15, 2024 13:32, edited 5 times in total.
Reason: Added note.
Luxan
Posts: 241
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: Radon

Post by Luxan »

Hidden posts, yes annoying.

I'm following your software development, the use of types is also found in the programming language Lean.
I'm also familiar with there use at the machine code level.

I managed to get CUDA set up upon my computer, without the use of a virtual environment.
Using that, I ran Python neural network code to find the optimum weights for the filter after the Radon
Transform . Using these weights for inference, I was indeed able to reconstruct a more accurate image.
Doesn't seem very intuitive though.

Using different geometries for the sources and detectors, the reconstructed image may not exhibit unwanted
edge effects .
dafhi
Posts: 1682
Joined: Jun 04, 2005 9:51

Re: Radon

Post by dafhi »

i did a search on 1d sharpen, found a post mentioning
origin [0 1 0] + sharp [-1 2 -1]

but it doesn't do radon justice. i might have an idea of how to proceed next:

make a delta of source and inverse radon, then analyze radon of that.

[edit]
by "analyze radon" i mean compare a radon angle sample with corresponding original angle result.

some transformation of original radon will maybe need to equal that. working on it
dafhi
Posts: 1682
Joined: Jun 04, 2005 9:51

Re: Radon

Post by dafhi »

unless i'm missing something, i seem to be reconstructing an image nearly perfectly, just not the original image

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"
'#include "/inc/fourier.bas"
'
' Define the size of the image
const IMG_WIDTH = 256 '128
const IMG_HEIGHT = 256 '128
const N_ANGLES = 38

const as single pi = 4 * atn(1)

#define min( a, b)    iif( (a)<(b), (a), (b) )
#define max( a, b)    iif( (a)>(b), (a), (b) )

#define dbl  as double '' reduced text
#define sng  as single
'
' -------------------------- declarations -----------------
'
Declare sub GenerateTestImage(image() as single,choice as integer=0)
Declare sub Radon_And_Inverse(image() as single, image1() as single, off_x as long = 0, off_y as long = 0)
declare sub _show( im() as single, x_offset as long = 0, y_offset as long = 0)
declare sub gen_delta( src() sng, xform0() sng, del() sng )
'
' ---------------------------------------------------------------
'
'       Main program
'
dim as single image(), constructed1()
dim as single delta(), constructed2()
'
dim as integer show_x = img_width
dim as integer show_y = img_height
' ---------------------------------------------------------------
'
GenerateTestImage(image())
_show image()

redim constructed1( img_width-1, img_height-1 )
Radon_And_Inverse( image(), constructed1(), show_x )
locate 40,40
print "difference"
sleep 1000
gen_delta image(), constructed1(), delta()
_show delta(), IMG_WIDTH, IMG_HEIGHT
line (img_width * 1.5, img_height * .5) - (img_width * .5, img_height * 1.5)
sleep 1500
redim constructed2( img_width-1, img_height-1 )
Radon_And_Inverse( constructed1(), constructed2(), 0, show_y )
line (img_width * 1.5, img_height * 1.5) - (img_width * .5, img_height * 1.5)
sleep 1000
Radon_And_Inverse( delta(), constructed2(), 0, show_y )

sleep
end
'
' ===========================
    dim shared as integer x, y, t, r_index, c, wm, hm, max_rho
'
'
' ---------------------------------------------------------------
    dim shared as single rho,  cost, sint, ysin
  
#macro _at_angle_init()
    cost = cos(a)
    sint = sin(a)
    for y = 0 to hm
        ysin = (y+0/2) * sint
        for x = 0 to wm
            rho = (x+0/2) * cost + ysin
            r_index = rho'int(rho)
'    if r_index < 0 then ? r_index; " ";': sleep
'    if r_index < 0 or r_index > ubound(radon) then ? r_index; " ";': sleep
#endmacro

sub radon_at_angle( image() as single, radon() as single, a sng )
        _at_angle_init()
        radon(r_index) += image(x, y)
    next
  next
end sub

sub inverse_at_angle( des() as single, radon() as single, a sng )
      _at_angle_init()
      des(x, y) += radon(r_index)
    next
  next
end sub


sub Radon_And_Inverse( src() as single, des() as single, off_x as long, off_y as long )
  max_rho = sqr(IMG_WIDTH ^ 2 + IMG_HEIGHT ^ 2)
  
  wm = ubound(src, 1)
  hm = ubound(src, 2)
  
  'redim des( wm, hm )
  
  for t = 0 to n_angles - 1
  
    dim sng radon(-max_rho to max_rho - 1)
    dim sng a = t * pi / n_angles
    
    radon_at_angle src(), radon(), a
    inverse_at_angle des(), radon(), a
    
    #if 1
    _show des(), off_x, off_y
    sleep 1000 / n_angles
    #endif
  
  next
  
'sleep
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
            for x = 0 to IMG_WIDTH - 1
                rho = x * cost + ysin
                r_index = int(rho)
                radon(t, r_index) += image(x, y)
            next
        next
    next
    
  wm = ubound(image, 1)
  hm = ubound(image, 2)
'    
end sub
'
' ---------------------------------------------------------------
    dim shared as single _max, _min

sub get__min_and_max( im() as single )
    _max = im( lbound(im,1), lbound(im,2) ) '' 2024 June 1
    _min = _max
    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
end sub

sub remap( des() sng, src() sng )
    
    get__min_and_max des()
  dim sng des_min = _min
  dim sng des_max = _max
  
    get__min_and_max src()
  dim sng scalar = (_max - _min) / (des_max - des_min)
  
    for y = 0 to hm
  for x = 0 to wm
  des(x,y) -= des_min
  des(x,y) *= scalar
  next
  next
  
end sub

sub gen_delta( src() sng, xform() sng, del() sng )
  remap xform(), src()
  redim del( IMG_WIDTH - 1, IMG_HEIGHT - 1)
    for y as long = 0 to hm
  for x as long = 0 to wm
  del(x,y) = src(x,y) - xform(x,y)
  next
  next
end sub
'
' ---------------------------------------------------------------
'
sub _show( im() as single, x_offset as long, y_offset as long)
  get__min_and_max im()

  _max = 255.999 / abs(_max-_min)

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

'
' Function to generate a simple test image
sub GenerateTestImage(image() as single,choice as integer)
  
  screenres 2*IMG_WIDTH, 2*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 8 ' seed
  
  for t = 1 to 9
    var a = int(rnd*256)
    circle im, (rnd * img_width, rnd * img_height), 5 + rnd*50, rgb(a,a,a),,,,f
  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
Luxan
Posts: 241
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: Radon

Post by Luxan »

There's a fair bit going on there, let's outline
what's happening.

1. Generate test image, a collection of, variable density,
filled circles.

2. Attempt to reconstruct that image via radon and inverse radon
routines.

3. Obtain an image that represents the difference between the original
image, from 1 and the reconstructed image from 2 .

4. Use radon and inverse radon routines to reconstruct the first
reconstructed image .

5. Use radon and inverse radon routines to reconstruct the
difference image from 3 .

Comments

Obviously in a real world situation, where you're scanning a physical
object, you can't obtain an exact image of a cross section; unless you
actually slice through it. Maybe like a knife through Swiss cheese.
Hence you wouldn't be able to obtain the difference.

Applying the radon then inverse radon routines twice seems to improve
the reconstruction of the original image.
This might be equivalent, in the real world, to doing a physical
scan and reconstruction, then simulating a scan and reconstruction
upon that data.
dafhi
Posts: 1682
Joined: Jun 04, 2005 9:51

Re: Radon

Post by dafhi »

sharpen 1d and coefficients optimizer

Code: Select all

/' 
  
  radon transform with sharpen 1d and coefficients optimizer
  
  2024 June 30  by dafhi


    update

    new tool - scaled lerp with expon
  
    lerp "endpoints"
  1. baseline Radon and Inverse
  2. 1d-sharpened radon and Inverse
 
    for my
  1d sharpen i use array "coeffs" with last 2 elements reserved for lerp scale & pow
  (sub named scale_delta)
 
'/

' -------- derived from original --------------------------------
'
'      rd0004a.bas
'
'             sciwiseg@gmail.com , 2024 ; lgpl , if applicable
'
'             74% help from ChatGPT
'
' ----------------------------

const IMG_WIDTH = 192
const IMG_HEIGHT = 192

const N_ANGLES = 99

const as single pi = 4 * atn(1)

#define min( a, b)    iif( (a)<(b), (a), (b) )
#define max( a, b)    iif( (a)>(b), (a), (b) )

#define dbl  as double '' reduced text
#define sng  as single

function round(in as double, places as ubyte = 2) as string
  dim as integer mul = 10 ^ places
  return str(csng( int(in * mul + .5) / mul) )
End Function


  namespace array2d_util

dim as integer x, y, wm, hm

sub _redim( a() sng )
  if ubound(a,1) <> wm or ubound(a,2) <> hm then redim a(wm,hm)
end sub

dim sng ptr     ps, pd, p0
dim sng         _max, _min, _scalar
dim as double   sum

#macro _src_des_loopinit( ades, asrc )
  p0 = @ades(0,0): ps = @asrc(0,0)
  for pd = p0 to p0 + (wm+1)*(hm+1)-1
#endmacro

'' get array min and max
sub get_minmax( im() as single )
  _max = im( 0, 0 ): _min = _max
    _src_des_loopinit( im, im )
  if *pd > _max then: _max = *pd
  elseif *pd < _min then _min = *pd: endif
  next
  _scalar = 1 / (_max - _min)
end sub

'' array copy
sub cpy( des() sng, src() sng )
  _redim des()
    _src_des_loopinit( des, src )
  *pd = ps[pd - p0]
  next
end sub

'' add R to L
sub ary_add( L() sng, R() sng )
    _src_des_loopinit( L, R )
  *pd += ps[pd - p0]
  next
end sub

'' invert:  a(i) = max - a(i)
sub invert( a() sng )
  get_minmax a()
    _src_des_loopinit( a, a )
  *pd = _max - *pd
  next
end sub

'' negate
sub negate( a() sng )
    _src_des_loopinit( a, a )
  *pd = -*pd
  next
end sub

'' like normalize but if base0 = false, min can be non-zero
sub remap( des() sng, base0 as long = true )
  get_minmax des()
  dim sng local_min = -(base0 <> false) * _min
    _src_des_loopinit( des, des )
  *pd -= local_min
  *pd *= _scalar
  next
end sub

function f_err( a() sng, b() sng ) as single
  dim as long c:  sum = 0
  p0 = @a(0,0): ps = @b(0,0)
  for pd = p0 to p0 + (wm+1)*(hm+1)-1 step 189 '' yeah baby
    sum += sqr( ( *pd - ps[pd-p0] ) ^ 2 )
    c += 1
  next
  return sum / c
end function

sub scale_delta( des() sng, src() sng, scale sng = 2, _pow sng = 1 )
    _src_des_loopinit( des, src )
  dim sng delta = *pd - ps[pd - p0]
  *pd = ps[pd - p0] + sgn(delta) * abs( scale * delta ) ^ _pow
  next
end sub

end namespace

    sub _show( im() as single, x_offset as long = 0, y_offset as long = 0)
      using array2d_util
      get_minmax im()
      dim sng scale = 255.999 * _scalar
        for y = 0 to hm
      for x = 0 to wm
      var c = int( scale * (im(x, y) - _min) )
      pset (x + x_offset, y + y_offset), rgb(c,c,c)
      next
      next
    end sub


  namespace ns_RadonTransform
'
dim as integer  t, r_index, max_rho

' ---- hack for reduced text -----------------------------------------
#macro _at_angle_init()
  for y = 0 to hm
      dim sng ysin = y * sint
      for x = 0 to wm
          dim sng rho = x * cost + ysin
          r_index = rho'int(rho)
#endmacro
' --------------------------------------------------------------------

  dim as single cost, sint
  using array2d_util
 
  sub redim_radon2( radon2() as single )
    max_rho = sqr( (wm+1) ^ 2 + (hm+1) ^ 2 )
    redim radon2( -1*max_rho to 1*max_rho - 1 )
  end sub
 
  sub radon_at_angle( src() as single, radon() as single, a as single )
    redim radon(-max_rho to max_rho - 1)
    cost = cos(a)
    sint = sin(a)
      _at_angle_init()
    radon(r_index) += src(x, y)
    next
    next
  end sub
  
  sub inverse_at_angle( des() as single, radon() as single, a as single )
      _at_angle_init()
    des(x, y) += radon(r_index)
    next
    next
  end sub

      sub at_angle( des() as single, src() as single, radon() as single, a as single )
        radon_at_angle src(), radon(), a
        inverse_at_angle des(), radon(), a
      end sub

    sub Radon_And_Inverse( des() as single, src() as single, radon() as single, off_x as long = 0, off_y as long = 0)
      max_rho = sqr( (wm+1) ^ 2 + (hm+1) ^ 2 )
        for t = 0 to n_angles - 1
      dim sng a = t * pi / n_angles
      at_angle des(), src(), radon(), a
      next
      remap des()
    end sub

  dim as single er, er_best
  
'' -- A main sub
'
sub Radon_And_FilteredInverse( des() sng, src() sng, radon() sng, radon2() sng, _
  my_filt as sub( ()sng, ()sng ), a0 sng = 0)

  redim_radon2 radon2()
  
    for t = 0 to n_angles - 1
  dim sng a = 2 * t * pi / n_angles + a0
  radon_at_angle src(), radon(), a
  my_filt( radon2(), radon() )
  inverse_at_angle des(), radon2(), a
  next
end sub

  sub cpy_coeffs( des() as single, src() as single )
    redim des(lbound(src) to ubound(src))
    for t = lbound(src) to ubound(src)
      des(t) = src(t)
    next
  end sub

dim as single coeffs(11), c_best(), shift_read_pos
 
 
sub my_sharpen( radon2() as single, radon() as single )

  dim as long l_src = lbound(radon)
  dim as long u_src = ubound(radon), idx_src
  
  shift_read_pos = -(ubound(coeffs)+1) / 4.5 '' experiment
  
    for i as long = l_src to u_src
  radon2(i) = 0
 
  '' last 2 coefficients i use for non-sharpen things
  for icoeff as long = lbound(coeffs) to ubound(coeffs) - 2
  
  idx_src = max( l_src, min( u_src, icoeff + shift_read_pos + i ) )
  
  radon2(i) += coeffs( (icoeff) ) * radon( idx_src )
  next
  next

end sub

  function f_rnd as single
    return 12 * (rnd - .5) '' wide search space :D
  end function
  
  sub coeff_append( idx as long, v as single )
    if rnd < .55 then exit sub
    coeffs(idx) = v
  end sub

sub randomize_
  for t = lbound(coeffs) to ubound(coeffs)
    coeff_append t, f_rnd
  next
end sub

end namespace

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

sub GenerateTestImage( image() as single )
  
  screenres 2*IMG_WIDTH, 2*IMG_HEIGHT, 32 '' 2024 June 1
  
  using array2d_util
  wm = IMG_WIDTH - 1
  hm = IMG_HEIGHT - 1
  
  dim as any ptr im = imagecreate( img_width, img_height )
  
  randomize 1 ' seed
  
  for x = 1 to 6
    var a = int(rnd*256)
    circle im, (rnd * img_width, rnd * img_height), 25 + rnd*50, rgb(a,a,a),,,,f
    a = int(rnd*256)
    
    #define ix rnd*img_width
    #define iy rnd*img_height
    line im, (ix,iy)-(ix,iy), rgb(a,a,a), bf
    
    #undef ix
    #undef iy
  next
  
  _redim image()
  
      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
  next
  next
  
  remap image()
  
  if imageinfo(im) = 0 then imagedestroy im

end sub
'
' ---------------------------------------------------------------
'
dim shared as integer show_plotx = img_width
dim shared as integer show_ploty = img_height
'
' ---------------------------------------------------------------
'
'       Main program
'
' ---------------------------------------------------------------
'
dim as single   image(), constructed1(), constructed2()
dim as single   delta()
dim as single   radon(), radon2()

GenerateTestImage(image())
_show image()

randomize

using array2d_util
using ns_RadonTransform

_redim constructed1()
Radon_And_Inverse constructed1(), image(), radon(), show_plotx
_show constructed1(), show_plotx

_redim constructed2()
randomize_
Radon_And_FilteredInverse constructed2(), image(), radon(), radon2(), @my_sharpen
er_best = f_err( constructed2(), image() )
cpy_coeffs c_best(), coeffs()

dim as single scalar, _pow', repurposed_as_angle0
dim as string kstr
dim as double t = timer, t_demo_end = t + 7, t_showinfo = t

do

  cpy_coeffs coeffs(), c_best()
  randomize_
  
  _redim constructed2()
    var u = ubound( coeffs )
  '  repurposed_as_angle0 = coeffs(u-2)
  Radon_And_FilteredInverse constructed2(), image(), radon(), radon2(), @my_sharpen', repurposed_as_angle0
 
  ' experimental
    scalar = coeffs(u)
    _pow = .1 + abs(coeffs(u-1)) / 1.5
  scale_delta constructed2(), constructed1(), scalar, _pow
  remap constructed2()
  er = f_err( constructed2(), image() )

  if er < er_best then
    er_best = er
    cpy_coeffs c_best(), coeffs()
    _show constructed2(), show_plotx, show_ploty
    locate 3 + 0, 41: print "coeffs"
    for t = 0 to ubound(coeffs)
      locate 5 + t, 41
      print round(coeffs(t), 3)
    next
    windowtitle "err " + round(er, 4)
  endif
  
  t = timer
  if t > t_showinfo then
    locate 2,2:  t_showinfo += 1.5
    ? "timer "; round(t_demo_end - t)
  endif
  sleep 15
  kstr = inkey
loop until kstr = chr(27) or t > t_demo_end

locate 2,2
? "fin!        "

sleep

- old post
i'm pretty much clueless but i built some tooling

Code: Select all

/' 
  radon transform retooling demo - 2024 June 27  by dafhi

    quick note:  remappings occur at the end of subs
  GenerateTestImage()
  Radon_And_Inverse()

'/

' -------- derived from original --------------------------------
'
'      rd0004a.bas
'
'             sciwiseg@gmail.com , 2024 ; lgpl , if applicable
'
'             74% help from ChatGPT
'
' ----------------------------

const IMG_WIDTH = 256 '128
const IMG_HEIGHT = 256 '128
const N_ANGLES = 38

const as single pi = 4 * atn(1)

#define min( a, b)    iif( (a)<(b), (a), (b) )
#define max( a, b)    iif( (a)>(b), (a), (b) )

#define dbl  as double '' reduced text
#define sng  as single

' ==============================================================

  namespace array2d_util
'
dim as integer x, y, wm, hm
'
sub _redim( a() sng )
  if ubound(a,1) <> wm or ubound(a,2) <> hm then redim a(wm,hm)
end sub
'
dim sng ptr    ps, pd, p0
dim sng        _max, _min, _scalar
'
#macro _src_des_loopinit( ades, asrc )
  p0 = @ades(0,0): ps = @asrc(0,0)
  for pd = p0 to p0 + (wm+1)*(hm+1)-1
#endmacro

'' get array min and max
sub get_minmax( im() as single )
  _max = im( wm, hm ): _min = _max
    _src_des_loopinit( im, im )
  if *pd > _max then: _max = *pd
  elseif *pd < _min then _min = *pd: endif
  next
  _scalar = 1 / (_max - _min)
end sub

'' array copy
sub cpy( des() sng, src() sng )
  _redim des()
    _src_des_loopinit( des, src )
  *pd = ps[pd - p0]
  next
end sub

'' add R to L
sub ary_add( L() sng, R() sng )
    _src_des_loopinit( L, R )
  *pd += ps[pd - p0]
  next
end sub

'' invert:  a(i) = max - a(i)
sub invert( a() sng )
  get_minmax a()
    _src_des_loopinit( a, a )
  *pd = _max - *pd
  next
end sub

'' negate
sub negate( a() sng )
    _src_des_loopinit( a, a )
  *pd = -*pd
  next
end sub

'' like normalize but if base0 = false, min can be non-zero
sub remap( des() sng, base0 as long = false )
  get_minmax des()
  dim sng local_min = -(base0 <> false) * _min
    _src_des_loopinit( des, des )
  *pd -= local_min
  *pd *= _scalar
  next
end sub

end namespace

    sub _show( im() as single, x_offset as long = 0, y_offset as long = 0)
      using array2d_util
      get_minmax im()
      dim sng scale = 255.999 * _scalar
        for y = 0 to hm
      for x = 0 to wm
      var c = int( scale * (im(x, y) - _min) )
      pset (x + x_offset, y + y_offset), rgb(c,c,c)
      next
      next
    end sub


  namespace ns_RadonTransform
'
dim as integer  t, r_index, max_rho

' ---- hack for reduced text --------------------
#macro _at_angle_init()
  for y = 0 to hm
      dim sng ysin = y * sint
      for x = 0 to wm
          dim sng rho = x * cost + ysin
          r_index = rho'int(rho)
#endmacro

  dim sng  radon()

  sub at_angle( des() as single, src() as single, a sng )
    
    redim radon(-max_rho to max_rho - 1)
  
    dim sng cost = cos(a)
    dim sng sint = sin(a)
    
    using array2d_util
    
      _at_angle_init()
    radon(r_index) += src(x, y)
    next
    next
    
      _at_angle_init()
    des(x, y) += radon(r_index)
    next
    next

  end sub

sub Radon_And_Inverse( des() as single, src() as single, off_x as long = 0, off_y as long = 0)
  
  using array2d_util

  max_rho = sqr( (wm+1) ^ 2 + (hm+1) ^ 2 )
  
    for t = 0 to n_angles - 1
  dim sng a = t * pi / n_angles
  at_angle des(), src(), a
  _show des(), off_x, off_y
  sleep 300 / n_angles
  next
  
  remap des()
  
end sub

end namespace

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

sub GenerateTestImage(image() as single)
  
  screenres 2*IMG_WIDTH, 2*IMG_HEIGHT, 32 '' 2024 June 1
  
  using array2d_util
  wm = IMG_WIDTH - 1
  hm = IMG_HEIGHT - 1
  
  dim as any ptr im = imagecreate( img_width, img_height )
  
  randomize 8 ' seed
  
  for x = 1 to 1
    var a = int(rnd*256)
    circle im, (rnd * img_width, rnd * img_height), 25 + rnd*50, rgb(a,a,a),,,,f
  next
  
  _redim image()
  
      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
  next
  next
  
  remap image()
  
  if imageinfo(im) = 0 then imagedestroy im

end sub
'
' ---------------------------------------------------------------
'
dim shared as integer show_plotx = img_width
dim shared as integer show_ploty = img_height
'
' ---------------------------------------------------------------
'
'       Main program
'
' ---------------------------------------------------------------
'
dim as single   image(), constructed1()
dim as single   delta(), constructed2()

GenerateTestImage(image())
_show image()

using array2d_util

_redim constructed1()
ns_RadonTransform.Radon_And_Inverse constructed1(), image(), show_plotx

sleep 1000
_redim delta()
ns_RadonTransform.at_angle delta(), constructed1(), 0
_show delta(), show_plotx, show_ploty

sleep 1000
negate delta()
_show delta(), show_plotx, show_ploty

sleep 1000
line (img_width * .5, img_height * .5) - (img_width * 1.5, img_height * 1.5)
sleep 1000

ns_RadonTransform.at_angle delta(), image(), 0
_show delta(), show_plotx, show_ploty

sleep 500
locate 2,2
? "fin!"

sleep
end

Last edited by dafhi on Jun 30, 2024 11:19, edited 1 time in total.
deltarho[1859]
Posts: 4345
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Radon

Post by deltarho[1859] »

If you cannot see previous posts whilst composing, run another instance of the forum.

Alternatively, compose offline, as I do most of the time, especially with longish posts where I may save them and return a few times to add more before the release post.
dafhi
Posts: 1682
Joined: Jun 04, 2005 9:51

Re: Radon

Post by dafhi »

it's a thing of a post being invisible, aka some kind of forum bug. well it could be anything. a collusion between gov't and aliens.

i have composed offline but when i'm inspired, stuff usually flows. i often find myself sitting for minutes trying to word something so it's either somewhat clear and hopefully thoughtful.

here's my latest update. sometimes my coefficient search writes new data that looks questionable so i thought why not write a human pick-and-choose. it's kinda fun

Code: Select all

/' 
  
  radon transform with sharpen 1d and "pick fave"
  
  2024 June 29  by dafhi

'/

' -------- derived from original --------------------------------
'
'      rd0004a.bas
'
'             sciwiseg@gmail.com , 2024 ; lgpl , if applicable
'
'             74% help from ChatGPT
'
' ----------------------------

const IMG_WIDTH = 192
const IMG_HEIGHT = 192

const N_ANGLES = 99

const as single pi = 4 * atn(1)

#define min( a, b)    iif( (a)<(b), (a), (b) )
#define max( a, b)    iif( (a)>(b), (a), (b) )

#define dbl  as double '' reduced text
#define sng  as single

function round(in as double, places as ubyte = 2) as string
  dim as integer mul = 10 ^ places
  return str(csng( int(in * mul + .5) / mul) )
End Function

  const rect_border = 1


  namespace array2d_util

dim as integer x, y, wm, hm

sub _redim( a() sng )
  if ubound(a,1) <> wm or ubound(a,2) <> hm then redim a(wm,hm)
end sub

dim sng ptr     ps, pd, p0
dim sng         _max, _min, _scalar
dim as double   sum

#macro _src_des_loopinit( ades, asrc )
  p0 = @ades(0,0): ps = @asrc(0,0)
  for pd = p0 to p0 + (wm+1)*(hm+1)-1
#endmacro

'' get array min and max
sub get_minmax( im() as single )
  _max = im( 0, 0 ): _min = _max
    _src_des_loopinit( im, im )
  if *pd > _max then: _max = *pd
  elseif *pd < _min then _min = *pd: endif
  next
  _scalar = 1 / (_max - _min)
end sub

'' array copy
sub cpy( des() sng, src() sng )
  _redim des()
    _src_des_loopinit( des, src )
  *pd = ps[pd - p0]
  next
end sub

'' add R to L
sub ary_add( L() sng, R() sng )
    _src_des_loopinit( L, R )
  *pd += ps[pd - p0]
  next
end sub

'' invert:  a(i) = max - a(i)
sub invert( a() sng )
  get_minmax a()
    _src_des_loopinit( a, a )
  *pd = _max - *pd
  next
end sub

'' negate
sub negate( a() sng )
    _src_des_loopinit( a, a )
  *pd = -*pd
  next
end sub

'' like normalize but if base0 = false, min can be non-zero
sub remap( des() sng, base0 as long = true )
  get_minmax des()
  dim sng local_min = -(base0 <> false) * _min
    _src_des_loopinit( des, des )
  *pd -= local_min
  *pd *= _scalar
  next
end sub

function f_err( a() sng, b() sng ) as single
  dim as long c:  sum = 0
  p0 = @a(0,0): ps = @b(0,0)
  for pd = p0 to p0 + (wm+1)*(hm+1)-1 step 19 '' 
    sum += ( ( *pd - ps[pd-p0] ) ^ 2 )
    c += 1
  next
  return sum / c
end function

end namespace

    sub _show( im() as single, x_offset as long = 0, y_offset as long = 0)
      using array2d_util
      get_minmax im()
      dim sng scale = 255.999 * _scalar
        for y = 0 to hm
      for x = 0 to wm
      var c = int( scale * (im(x, y) - _min) )
      pset (x + x_offset, y + y_offset), rgb(c,c,c)
      next
      next
    end sub


  namespace ns_RadonTransform
'
dim as integer  t, r_index, max_rho

' ---- hack for reduced text -----------------------------------------
#macro _at_angle_init()
  for y = 0 to hm
      dim sng ysin = y * sint
      for x = 0 to wm
          dim sng rho = x * cost + ysin
          r_index = rho'int(rho)
#endmacro
' --------------------------------------------------------------------

  dim as single cost, sint
  using array2d_util
 
  sub redim_radon2( radon2() as single )
    max_rho = sqr( (wm+1) ^ 2 + (hm+1) ^ 2 )
    redim radon2( -1*max_rho to 1*max_rho - 1 )
  end sub
 
  sub radon_at_angle( src() as single, radon() as single, a as single )
    redim radon(-max_rho to max_rho - 1)
    cost = cos(a)
    sint = sin(a)
      _at_angle_init()
    radon(r_index) += src(x, y)
    next
    next
  end sub
  
  sub inverse_at_angle( des() as single, radon() as single, a as single )
      _at_angle_init()
    des(x, y) += radon(r_index)
    next
    next
  end sub

  dim as single er, er_best
  
'' -- A main sub
'
sub Radon_And_FilteredInverse( des() sng, src() sng, radon() sng, radon2() sng, _
  my_filt as sub( ()sng, ()sng ), a0 sng = 0)
  
  redim_radon2 radon2()
  
    for t = 0 to n_angles - 1
  dim sng a = t * pi / n_angles + a0
  radon_at_angle src(), radon(), a
  my_filt( radon2(), radon() )
  inverse_at_angle des(), radon2(), a
  next
  
  remap des()
  er = f_err( des(), src() )
  
end sub

  sub cpy_coeffs( des() as single, src() as single )
    if ubound(des) <> ubound(src) then redim des(lbound(src) to ubound(src))
    for t = lbound(src) to ubound(src)
      des(t) = src(t)
    next
  end sub

'' coefficients
  dim as single coeffs(5), c_best()

sub my_sharpen( radon2() as single, radon() as single )

  dim as long l_src = lbound(radon)
  dim as long u_src = ubound(radon), idx_src
  
    for i as long = l_src to u_src
  radon2(i) = 0
  for ic as long = -ubound(coeffs) to ubound(coeffs)
  idx_src = max( l_src, min( u_src, ic+i ) )
  radon2(i) += coeffs( abs(ic) ) * radon( idx_src )
  next
  next

end sub

  function f_rnd as single
    return -rnd * 1.2
  end function
  
  dim as boolean  something_changed
  
  sub coeff_append( idx as long, v as single )
    if rnd < .65 then exit sub
    something_changed = true
    coeffs(idx) = v
  end sub

sub randomize_
  dim as long u = ubound(coeffs)
  something_changed = false
  coeff_append 1, f_rnd
  coeff_append u, f_rnd
  for t = 2 to u-1
    #if 1
      coeff_append t, f_rnd
    #else
      lerp = (t-1)*f
      coeff_append t, c(1) + lerp * (c(u) - c(1))
    #endif
  next
  sum = 0
  for t = 1 to u
    sum += coeffs(t)
  next
  if something_changed then coeffs(0) = -2 * sum - .1 + rnd
  if (rnd < .25) then
'  if (rnd < .2) or something_changed = false then
  coeffs(0) += rnd * 12
  endif
end sub

end namespace

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

sub helpscreen
  locate 1,1
  ?
  ? "pick fave"
'  ?
'  ? "ctrL - s:  save coeffcients.txt"
end sub
  dim shared as boolean b_showhelp

sub toggle_help
  b_showhelp = not b_showhelp
  if b_showhelp then
    screencopy 0,1
    helpscreen
  else
    screencopy 1,0
  endif
end sub
 
const rcw = IMG_WIDTH + rect_border*2
const rch = IMG_Height + rect_border*2

Type mouse '' https://www.freebasic.net/wiki/KeyPgGetmouse
    As Long res
    As Long x, y, wheel, clip
    Union
        buttons As Long
        Type
            Left:1 As Long
            Right:1 As Long
            middle:1 As Long
        End Type
    End Union
End Type

type t_coeff_rect
  as single     coeffs( ubound(ns_RadonTransform.coeffs) ) = {1}
  as long       keep
end type


  namespace ns_selection

dim as t_coeff_rect   coef_rc(0 to 2, 0 to 1)
dim as mouse          m
dim as long           plotx, ploty, a, ix, iy
dim as boolean        updating = true
  
  using ns_RadonTransform

sub do_mouse( )
  static as long buttons_previous
  
  m.res = GetMouse( m.x, m.y, m.wheel, m.buttons, m.clip )
  
  if m.buttons = 0 then
    if buttons_previous > 0 then
      ix = m.x \ rcw '' integer divide
      iy = m.y \ rch '' integer divide
      plotx = rcw * ix
      ploty = rch * iy
      
      cpy_coeffs coeffs(), coef_rc(ix,iy).coeffs()
      for y as long = 0 to ubound( coef_rc, 2 )
      for x as long = 0 to ubound( coef_rc, 1 )
      coef_rc( x, y ).keep = 0
      cpy_coeffs coef_rc(x,y).coeffs(), coeffs()
      next
      next
      
      ix = ubound( coef_rc, 1 ) \ 2 '' integer divide
      iy = ubound( coef_rc, 2 ) \ 2
      coef_rc( ix, iy ).keep = -1
      updating = true
    endif
  endif
  buttons_previous = m.buttons
end sub

dim as single       image(), constructed1()
dim as single       radon(), radon2()

sub draw_all
  if not updating then exit sub
'  locate 1,1
  for y as long = 0 to ubound( coef_rc, 2 )
  for x as long = 0 to ubound( coef_rc, 1 )
    cpy_coeffs coeffs(), coef_rc(x,y).coeffs()
    if coef_rc(x,y).keep then
    else
      randomize_
      cpy_coeffs coef_rc(x,y).coeffs(), coeffs()
    endif
    _redim constructed1()
      dim sng repurposed_as_angle0 = coef_rc(x,y).coeffs(1)
    Radon_And_FilteredInverse constructed1(), image(), radon(), radon2(), @my_sharpen, repurposed_as_angle0
    _show constructed1(), x*(rcw) + rect_border, y*(rch) + rect_border
  next
  next
  updating = false
end sub

end namespace

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

sub GenerateTestImage( image() as single )
  
  var n_pages = 2
  
    screenres _
  rcw * (ubound( ns_selection.coef_rc, 1 )+1), _
  rch * (ubound( ns_selection.coef_rc, 2 )+1), 32, n_pages
 
  using array2d_util
  wm = IMG_WIDTH - 1
  hm = IMG_HEIGHT - 1
  
  dim as any ptr im = imagecreate( img_width, img_height )
  
  randomize 1 ' seed
  
  for x = 1 to 6
    var a = int(rnd*256)
    circle im, (rnd * img_width, rnd * img_height), 25 + rnd*50, rgb(a,a,a),,,,f
      a = int(rnd*256)
    line im, (rnd*img_width,rnd*img_height)-(rnd*img_width,rnd*img_height), rgb(a,a,a), bf
  next
  
  _redim image()
  
      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
  next
  next
  
  remap image()
  
  if imageinfo(im) = 0 then imagedestroy im
 
end sub

'
' ---------------------------------------------------------------
'
'       Main program
'
' ---------------------------------------------------------------
'

GenerateTestImage ns_selection.image()

helpscreen
sleep 1000

randomize

using array2d_util
using ns_RadonTransform

dim as string kstr

do
  sleep 15
  
    #if 1
  ns_selection.do_mouse
  ns_selection.draw_all
  #endif
  
  kstr = inkey
  dim as long klen = len(kstr)
'  windowtitle str(klen)
  
  if klen = 2 then
    if kstr[1] = 59 then toggle_help: continue do
    ns_selection.updating = true
  elseif klen = 1 then
    select case kstr
    case chr(27) '' ESC
      if b_showhelp then
        toggle_help
      else
        exit do
      endif
    case else
      ns_selection.updating = true
    end select
  endif

loop

locate 2,2
? "fin!"
sleep 1500

deltarho[1859]
Posts: 4345
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Radon

Post by deltarho[1859] »

dafhi" wrote:i often find myself sitting for minutes trying to word something so it's either somewhat clear and hopefully thoughtful.
Writer's block.

I get that. There are many ways to combat it. What works for me is to save what I have done and then do something else for a while. Very typically when I return, the block has gone, and I then continue to 'rattle on'. :)
Luxan
Posts: 241
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: Radon

Post by Luxan »

Often I attempt to do too many tasks at once, the human mind doesn't perform well in that instance.

Lets review what we've been doing.

From a ChatGPT session I introduced code for
the Radon and Inverse Radon method .
Though not verified, the geometry of the source
and detector is most likely a single source and
single detector. Otherwise a single source and
multiple detectors, in a fan beam arrangement.

Then dafhi provided code for an efficient accumulator
based transform, possibly similar to this :

https://en.wikipedia.org/wiki/Hough_transform

I provided some code for digital filtering, via the fft
and inverse fft , upon the rows of the radon transformed
data. I believe I also provided code for the 2d fft and 2d
ifft; in both instances I used the fbmath library.
Mention was made of CUDA to develop a Neural Network for
inference, or possibly a digital filter.


Most recently other authors have provided interesting
techniques to enhance the detection of objects within
the scanner space.

As a slight detour I want to discuss scanner geometry.
Firstly, consider a discrete area, equivalent to a
2d array. The pset command helps us visualize this.

To completely fill a circle the size and number of divisions
around the circumference, of radius r, must be calculated.

Here's some code to illustrate this :

Code: Select all


'
'
'      Line segments & circles , using pset
'
'      Length of a circle = 2*pi*r
'      Therefore, number of [ unit ] divisions is 2*pi*r
'
'   crc_d1.bas
'
'      Radon Transform of a filled circle :   2*sqr(r^2 - t^2) , |t| <= r*cos(theta) 
'                                             0                , otherwise

'
' f1(r,t):=2*sqrt(r^2 - t^2)$
' f2(r,t):=if(abs(t)<r) then 1/f1(r,t) else 1$   normalize radon values prior to inverse ?
'


const pi=4*atn(1)
dim as single r,x,y,z,dt,x1,y1,s,ds,px,py,sum
dim as integer i,j,n,m

r=180
screen 12
line(0,0)-(400,400),11,b
'circle(200,200),r,14,0,6.28
m=8
dt=pi/(m*r)
n=m*r
ds=1/(2*r)

for i=0 to n step 1
    y=200+r*sin(dt*i)
    x=200+r*cos(dt*i)
    pset(x,y),11
    y1=200+r*sin(dt*i+pi)
    x1=200+r*cos(dt*i+pi)
    pset(x1,y1),14
   ' line(x,y)-(x1,y1),10
    sum=0
    for s=0 to 1 step ds
      px=x+(x1-x)*s
      py=y+(y1-y)*s
      pset(px,py),10
      sum=sum+1
    next s
    sleep 10
    
    
    
next i



sleep
end
'
' ======================================================================
'


dafhi
Posts: 1682
Joined: Jun 04, 2005 9:51

Re: Radon

Post by dafhi »

deltarho[1859] wrote: Jun 29, 2024 20:40
dafhi" wrote:i often find myself sitting for minutes trying to word something so it's either somewhat clear and hopefully thoughtful.
Writer's block.

I get that. There are many ways to combat it. What works for me is to save what I have done and then do something else for a while. Very typically when I return, the block has gone, and I then continue to 'rattle on'. :)
an escape like that for me can be automatic. for me it can help also to have multiple projects, aka i'll switch when i've had enough. Luxan's comment about staying focused on a single-project .. heh

Luxan i may know where you are going with a suggested new search pattern. One of the videos i watched mentioned fans as being more efficient than beams

[edit] if you're reading this i updated my coefficients optimizer 5 or so posts up
Last edited by dafhi on Jun 30, 2024 11:22, edited 1 time in total.
Luxan
Posts: 241
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: Radon

Post by Luxan »

There are a few directions we might take when delving into the geometry of detectors.
How the data is collected, and possibly rearranged for filtering, depends heavily upon the geometry of
the scanner.

The original code I wrote, independent of ChatGPT, used multiple parallel beams and adjacent detectors.
I'd need to compare that with the code that has been posted here to determine just how valid that is.

And you're right, switching to an unrelated task can give the mind an opportunity to sort through the present
difficulties. Sometimes there's a backlog of unresolved tasks, in that instance just do something that doesn't
require undue brain power; like pegging out the washing.
Luxan
Posts: 241
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: Radon

Post by Luxan »

Correction, my original code, before the ChatGPT example, used a parallel beam arrangement
with multiple sources and the corresponding , Opposite, detectors .

The exact line integral was replaced with a summation, or numerical integration, as even for relatively
simple objects, 2d functions, the integration was too complicated for the Computer Algebra Systems
I have available.
dafhi
Posts: 1682
Joined: Jun 04, 2005 9:51

Re: Radon

Post by dafhi »

i hope you find the right tooling to satisfy your needs.

my Gen 4 accumulator is done. i now store a decently accurate floating point average in a byte. it's awesome.

though in my random ray shooter, i discovered i'd "need" to track write counts per pixel,
which Gen 4 "doesn't have enough space for,"

but a hack works surprisingly well

Code: Select all

  dim as long n = min( 6, ray_index )

here's Gen 4

Code: Select all

sub addv( byref b as ubyte, v as single, iteration as single )
  /'
    Accumulator (Gen 4) - 2024 July 3  by dafhi
    
    usage:  ns_accum.addv  ret_ubyte, val_0_to_1, iteration
  '/
  
  dim as single delta = int(v * 255.999) - b

  dim as single lerp = 1/iteration
 
  b = (b + delta * lerp)

end sub
Post Reply