#
# Hough Transformation Example
#

package require Tk
package require vectcl
namespace import vectcl::vexpr

set C_RANGE     125                             ;# the radius in the hough volume
set C_WIDTH     [vexpr {2 * C_RANGE}]           ;# width of hough volume
set C_HEIGHT    [vexpr {2 * C_RANGE}]           ;# height of hough volume

set R_MIN       100                             ;# minimum of valid radius
set R_MAX       550                             ;# maximum of valid radius
set R_WIDTH     [vexpr {R_MAX - R_MIN}]

set DEG_STEP     15     ;# search for edges every n degrees

set hotspotX    739     ;# X coordinate in the white spot
set hotspotY    605     ;# Y coordinate in the white spot

set scrnCenterY 675     ;# Y coordinate of the center of the screen
set scrnLeft    0       ;# X coordinate of the left bound of the screen
set scrnRight   0       ;# X coordinate of the left bound of the screen
set scrnWidth   0       ;# width of screen bounds
set scrnBrMin   250     ;# minimum screen brightness = 25 % (per mill)

set threshold   100     ;# minimum brightness difference = 10 % (per mill)

###############################################################################

proc MarkCenter {xc yc} {
  vexpr {
    x1 = xc -25
    x2 = xc +25
    y1 = yc -25
    y2 = yc +25
  }
  .c create line $x1 $yc $x2 $yc -width 1 -fill #ff00ff
  .c create line $xc $y1 $xc $y2 -width 1 -fill #ff00ff
  update
}


proc MarkScreenBounds {} {
  global scrnCenterY scrnLeft scrnRight scrnWidth
 
  puts "screen bounds: $scrnLeft..$scrnRight = $scrnWidth pixel"
  vexpr {
    y  = scrnCenterY
    y1 = y -50
    y2 = y +50
    x1 = scrnLeft -30
    x2 = scrnLeft
    x3 = scrnRight
    x4 = scrnRight +30
  }
  .c create line $x1 $y  $x2 $y  -width 2 -fill cyan -arrow last
  .c create line $x3 $y  $x4 $y  -width 2 -fill cyan -arrow first
  .c create line $x2 $y1 $x2 $y2 -width 1 -fill cyan
  .c create line $x3 $y1 $x3 $y2 -width 1 -fill cyan
}


proc MarkPoint {alpha r xp yp br} {
  puts [format "  alpha=%5.3f r=%3d P=(%4d,%4d) br=%d" $alpha $r $xp $yp $br]
  vexpr {
    x1 = xp - 7
    x2 = xp + 7
    y1 = yp - 7
    y2 = yp + 7
  }
  .c create oval $x1 $y1 $x2 $y2 -outline #ff0000 -fill {} -width 2
  update
}


proc MarkImprint {h_max m_x m_y m_r} {
  puts [format "h_max=%.0f @ P=(%u,%u) R=%u" $h_max $m_x $m_y $m_r]
  vexpr {
    x1 = m_x - m_r
    x2 = m_x + m_r
    y1 = m_y - m_r
    y2 = m_y + m_r
  }
  .c create oval $x1 $y1 $x2 $y2 -outline #00ff00 -fill {} -width 3
  vexpr {
    x1 = m_x - 5
    x2 = m_x + 5
    y1 = m_y - 5
    y2 = m_y + 5
  }
  .c create oval $x1 $y1 $x2 $y2 -fill #00ff00
  update
}


proc ShowHB {hbval} {
  .c create text 50 50 -anchor nw -font {Arial -24 bold} -fill white \
    -text [format "HB=%.1f" $hbval]
}

###############################################################################

namespace import tcl::mathfunc::int
namespace import tcl::mathfunc::round
namespace import tcl::mathfunc::hypot

###############################################################################

proc CalcHB {rPixel} {
  global scrnWidth
  
  vexpr {
    screenDiameter = 135.1                        # diameter of the original screen in mm
    enlargement    = 20                           # optical factor of the HB gauge
    force          = 29419.95                     # pressure force in N
    ballDiameter   = 10                           # ball diameter in mm
    PI             = 2 * asin(1.0)
  
    # calculate the diameter of the imprint in mm
    imprintDiameter = screenDiameter * 2 * rPixel / scrnWidth / enlargement

    imprintDepth = 0.0   # the depth of the imprint
    hbValue      = 0.0   # the brinell hardness value

    if imprintDiameter < ballDiameter {
      imprintDepth = abs(ballDiameter - hypot(ballDiameter, imprintDiameter)) / 2
      if imprintDepth > 0.0 {
        hbValue = force / (9.80665 * PI * ballDiameter * imprintDepth)
      }
    }
  }

  puts [format "D=%.3f mm, d=%.3f mm, h=%.6f mm, HB=%.1f" \
          $ballDiameter $imprintDiameter $imprintDepth $hbValue]
          
  ShowHB $hbValue
  return $hbValue
}

###############################################################################

proc houghtrans {} {
  global C_WIDTH C_HEIGHT R_WIDTH C_RANGE R_MIN R_MAX DEG_STEP \
         scrnCenterY scrnLeft scrnRight scrnWidth scrnBrMin \
         hotspotX hotspotY threshold 

  set t_start [clock milliseconds]
  puts "Houghtrans - Tcl edition"

  puts "reading picture data ..."
  image create photo hbcam -file hbcam.png
  set w [image width  hbcam]
  set h [image height hbcam]
  pack [canvas .c -width $w -height $h] -expand yes -fill both
  .c create image 0 0 -image hbcam -anchor nw
  update

  # copy image into pixel (a vectcl 2D vector)
  .c create line 0 0 $w 0 -width 3 -fill yellow -tags act_row
  vexpr {pixel = constfill(0,h,w)}              ;# create 2D pixel array
  for {set y 0} {$y < $h} {incr y} {            ;# for each row
    if {$y % 5 == 0} {
      .c coords act_row 0 $y $w $y              ;# show where we are (sometimes)
      update
    }
    for {set x 0} {$x < $w} {incr x} {          ;# for each column
      set red [lindex [hbcam get $x $y] 0]
      vexpr {pixel[y,x] = red / 0.255}          ;# set gray value (brightness) in per mill
    }
  }
  .c delete act_row
  MarkCenter $hotspotX $hotspotY

  puts "searching edges ..."
  vexpr {
    # allocate 3D hough volume
    houghdata = constfill(0,C_HEIGHT,C_WIDTH,R_WIDTH)
    
    # search for left screen bound
    searchLimit = w/4                           # search up to 1/4 of the picture
    scrnLeft = 0                                # init left screen bound
    x = 0                                       # start from left side
    while x < searchLimit {                     # search from left to right
      if pixel[scrnCenterY,x] > scrnBrMin {     # search for more than x% brightness
        scrnLeft = x                            # found left bound
        x = searchLimit                         # finish loop
      } else {
        x += 1
      }
    }
    # search for right screen bound
    scrnRight = 0
    searchLimit = 3*w/4                         # search up to 3/4 of the picture
    x = w-1                                     # start from right side
    while x > searchLimit {                     # search from right to left
      if pixel[scrnCenterY,x] > scrnBrMin {     # search for more than x% brightness
        scrnRight = x                           # found right bound
        x = searchLimit                         # finish loop
      } else {
        x -= 1
      }
    }
    scrnWidth = scrnRight - scrnLeft            # this is the reference width
    MarkScreenBounds()                          # show arrows for screen bounds

    # search for edge of impression ...
    fullCircle = 4 * asin(1.0)                  # 2 pi
    delta = fullCircle / (360.0 / DEG_STEP)     # angle step size    
    alpha = 0.0                                 # start at zero degrees
    while alpha < fullCircle {                  # for each angle ...
      br = 0                                    # init brightness
      xp = 0                                    # init x coordinate of resulting point
      yp = 0                                    # init y coordinate of resulting point

      fx = cos(alpha)                           # vector x-scaling
      fy = sin(alpha)                           # vector y-scaling

      br_list = constfill(-1.0, 30)             # make a vector for 30 brightness values

      l = h / 2                                 # not more than half height of picture
      m = l -1                                  # upper radius limit
      r = 0                                     # resulting radius

      for ri = 0:m { 
        x = hotspotX + round(fx * double(ri))   # calculate x coordinate
        y = hotspotY + round(fy * double(ri))   # calculate y coordinate
        
        if x > 0 && x < w && y > 0 && y < h {   # is this point in my area?
          db = 0                                # reset brightness difference
          br_list[29] = pixel[y,x]              # save brightness of point in per mill
          
          if br_list[0] != -1.0 {               # if we have enough samples, calculate the jump
            b1 = mean(br_list[0:24])            # inner brightness
            b2 = mean(br_list[25:29])           # outer brightness
            db = int(b2 - b1)                   # get brightness difference
          }
          br_list[0:28] = br_list[1:29]         # shift the brightness buffer

          if db > threshold {                   # minimum brightness reached
            xp = x                              # set resulting x coordinate
            yp = y                              # set resulting y coordinate
            br = db                             # set resulting brightness
            r = ri                              # set resulting radius
            ri = m                              # finish, leave the loop
          }
        }
      }

      if r > 0 && r < l {                       # Found edge!
        MarkPoint(alpha,r,xp,yp,br)             # show where we are
        
        for y = 0:C_HEIGHT-1 {                  # for each row in the hough volume ...
          for x = 0:C_WIDTH-1 {                 # for each column in the hough volume ...
            dx = xp - (hotspotX + x - C_RANGE)  # x-distance from hotspot point
            dy = (hotspotY + y - C_RANGE) - yp  # y-distance from hotspot point
            ri = int(hypot(dx, dy))             # calculate the radius of this point
            if ri >= R_MIN && ri < R_MAX {
              i = ri - R_MIN                    # radius index goes from R_MIN to R_MAX
              houghdata[y,x,i] += br            # add brightness value for this point
            }
          }
        }
      }
      
      alpha += delta                            # next angle
    }

    # searching the hough volume for it's maximum value ...
    h_max = -1
    m_x = -1
    m_y = -1
    m_r = -1
    for y = 0:C_HEIGHT-1 {                      # for each row in the hough volume ...
      for x = 0:C_WIDTH-1 {                     # for each column in the hough volume ...
        h_val = max(houghdata[y,x,:])           # search new maximum with vexpr max vector function
        if h_val > h_max {
          for r = 0:R_WIDTH-1 {                 # for each radius row in the hough volume ...
            h_val = houghdata[y,x,r]
            if h_val > h_max {                  # look for the greatest added brightness
              h_max = h_val                     # my maximum value (until now)
              m_x = hotspotX + x - C_RANGE      # copy x coordinate
              m_y = hotspotY + y - C_RANGE      # copy y coordinate
              m_r = r + R_MIN                   # copy radius
            }
          }
        }
      }
    }
    
    MarkImprint(h_max,m_x,m_y,m_r)              # show the resulting circle and midpoint
    
    CalcHB(m_r)                                 # calc brinell hardness with the determined radius
  }

  set t_end [clock milliseconds]
  puts [format "time=%.3f sec" [expr ($t_end - $t_start) / 1000.0]]  
}

houghtrans