Image box-counting

extensions [
  gis
  matrix ]
globals
[divisions fractal-dim
  len walk-count old-box-count
  new-box-count x-axis-list y-axis-list
  box-size explore? iteration
  iterations-to-go slope lin-reg-eq
  any-new-boxes?
  r-square
  automatic-bcd?
]
breed         ;used only for taking box counting dimension of fractal examples.
[boxes box]
boxes-own
[past? present? future?]
patches-own [fractal?]

;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;import image;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;

to setup
  clear-all
  if image = "coastline"
  [
    import-pcolors  "Coastline.png"
    ask patches with [pcolor != white]
    [set pcolor black]
  ]
  if image = "tree"
   [
    import-pcolors  "Tree.png"

    ask patches with [pcolor != white]
    [set pcolor black]
  ]

  if image = "custom-image"
  [
    import-pcolors  "custom-image.png"   ;make sure your image is a square with colors that represent what you're interested in dimensioning (black and white is easiest)
  ]

  set explore? false
  set automatic-bcd? true
  set x-axis-list [ ]
  set y-axis-list [ ]
  ask turtles [die]
  set box-size initial-box-length
  set iteration 0

  reset-ticks
end 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;Box-Counting-Dim;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to bcd-go
  set len initial-box-length
  bcd-begin
end 

to bcd-begin
  ask boxes ; clears screen in case there are any old boxes
    [die]

  if automatic-bcd?
    [
      if box-size >= initial-box-length
      [ set automatic-bcd? false]
    ]

  if ticks != 0
    [set box-size box-size + increment
      set iteration iteration + 1
      set iterations-to-go 91 - iteration
    ]
  set old-box-count 0
  set new-box-count 1    ;eliminates an error for first round
  make-initial-box
  make-neighbors-count-patches
end 

to setup-bcd
  set box-size initial-box-length
  make-initial-box
end 

;makes a starter box at the beginning of each run with present? = true.
;This box will then be used to make boxes with future? = true
;which will be used for the next run.

to make-initial-box
  create-boxes 1
  ask boxes [
    set shape "square"
    set size box-size
    setxy  -60 0
    set heading 90
    set color red
    set past? false
    set present? true
    set future? false
  ]
end 

;makes a Moore neighborhood around the starter box and counts patches below each new box (exploit).
;If there are no new boxes with patches under them for a given run a box is sent outside the neighborhhod
;to cover non-contiguous patches (explore). If this box finds no new patches the run is complete.

to make-neighbors-count-patches
  ask boxes with [present? = true ]
    [make-neighbors]

  ask boxes with [future?  = true]
    [exploit]
  count-new-boxes
  if any-new-boxes?     = false
    [explore]

  if any-new-boxes? = false
    [calculate-box-counting-dimension]
  if any-new-boxes? = false and automatic-bcd? [
    bcd-begin
    stop ]
  update-past-present-future-states
  tick
  if any-new-boxes? = true
    [make-neighbors-count-patches]
end 

to make-neighbors
  hatch 1 [fd box-size rt 90
    set present? false set future? true
    hatch 1 [fd box-size rt 90
      set present? false set future? true
      hatch 1 [fd box-size
        set present? false set future? true
        hatch 1 [fd box-size rt 90
          set present? false set future? true
          hatch 1 [fd box-size
            set present? false set future? true
            hatch 1 [fd box-size rt 90
              set present? false set future? true
              hatch 1 [fd box-size
                set present? false set future? true
                hatch 1 [fd box-size
                  set present? false set future? true
                ]]]]]]]]
end 

to exploit
  if count boxes in-radius (box-size / 2) > 1  ; eliminates duplicates
    [die]

  if count patches-under-me = 0
    [ die ]
end 

to-report patches-under-me
  report  patches in-radius  ( (1.4 * size ) / 2 )  with [pcolor = color-value]
  ;   [let my-n-edge ycor + (0.5 * size)  ;; if my shape is a square
  ;   let my-s-edge ycor - (0.5 * size)
  ;   let my-w-edge xcor - (0.5 * size)
  ;   let my-e-edge xcor + (0.5 * size)
  ;   report patches with [pcolor = green and
  ;   (pycor - 0.5 < my-n-edge) and (pycor + 0.5 > my-s-edge) and
  ;   (pxcor + 0.5 > my-w-edge) and (pxcor - 0.5 < my-e-edge)
  ;   ]
  ;  ]
end 

to explore
  if count boxes with [present? = true] > 0 [
    ask patches with [pcolor = color-value] [
      ifelse count boxes in-radius  (  box-size ) = 0
      [set explore? true]
      [stop]
    ]
  ]

  if explore? [
    ask one-of boxes with [present? = true] [
      hatch 1 [
        set present? false set future? true
        move-to min-one-of patches with [pcolor = color-value and count boxes in-radius  ( box-size ) = 0 ]
        [distance myself]]
    ]
  ]
  count-new-boxes
  set explore? false
end 

to count-new-boxes
  set old-box-count new-box-count
  set new-box-count count boxes
  ifelse old-box-count = new-box-count
  [set any-new-boxes? false]
  [set any-new-boxes?  true]
end 

to update-past-present-future-states
  ask boxes [
    if present? = true
    [set past? true set present? false]
    if future?   = true
    [set future? false set present? true]
  ]
end 

to calculate-box-counting-dimension

  if count boxes >= 1 [     ; eliminates potential error message if setup is pressed during a box-counting procedure
    set-current-plot "Box Counting Plot"
    set-current-plot-pen "default"
    let no-boxes log (count boxes ) 10
    let scale (log ( 1 / box-size ) 10 )
    plotxy scale no-boxes
    set x-axis-list lput scale x-axis-list
    set y-axis-list lput no-boxes y-axis-list
  ]

  stop
end 

to fractal-dimension
  if ticks = 0 [
    if divisions > 0 [ ; eliminates potential error message
      let line-segments count turtles
      set fractal-dim precision(( log  line-segments  10 / log  divisions 10 ))3
      show line-segments
      show divisions
      show fractal-dim
    ]
  ]
  stop
end 


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;Linear Reg;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;Linear regression is used to find the 'best-fit' line
;through all the tuples (box-size,number of boxes) plotted in the "Number of Boxes vs. Scale" plot.
;The slope of the line is the box counting dimension.

to linear-regression

  if count boxes >= 1 [  ; eliminates potential error message if setup is pressed during a box-counting procedure [

    let regression matrix:regress matrix:from-column-list (list y-axis-list  x-axis-list)   ;using the regression tool from the matrix extension
                                                                                            ;setting y-intercept and slope (measure of goodness of fit)
    let y-intercept item 0 (item 0 regression)
    set slope item 1 (item 0 regression)
    set r-square item 0 (item 1 regression)


    ; set the equation to the appropriate string
    set lin-reg-eq (word (precision slope 3) " * x + " (precision y-intercept 3))


    set-current-plot "Box Counting Plot"
    set-current-plot-pen "pen-4"
    plot-pen-reset
    auto-plot-off
    plotxy plot-x-min (plot-x-min * slope + y-intercept)
    plot-pen-down
    plotxy plot-x-max (plot-x-max * slope + y-intercept)
    plot-pen-up

  ]
end