'''
Module Name:  Meanshift pseudo code

Created on May 24, 2016

@author: Bo Yang, Moritz Lennert, Markus Metz

GRASS GSoC2016 project, All Rights Reserved
'''

#==================================================================================
#               input parameters
#==================================================================================

# sampling window radius = spatial bandwidth in number of pixels
radius = 3.5

# range bandwidth (spectral bandwidth) for normalized difference between old and new values
# must be 0 < range bandwidth < 1
range_bandwidth = 0.1

# threshold to control when iterations should stop
# stop if the maximum difference between old and new values is < threshold
threshold_factor= 0.01

# maximum number of iterations
max_iter = 20


#===================================================================================
#                       input  Raster
#===================================================================================

# use existing functions in i.segment


#===================================================================================
#               gaussian weight difference (S=||B1-A1||)
#===================================================================================

# window size in number of pixels
window_size = (int)radius * 2 + 1

# initialize

# copy input raster to output raster
rastout = rastin

# iteration counter
iter = 0

# number of changed pixels with change > threshold
n_changes = 1

while (iter < max_iter && n_changes > 0 {

  # make output raster of last iteration input raster of this iteration
  rasttmp = rastin
  rastin = rastout
  rastout = rasttmp

  # go through all rows
  for (row = 0; row < nrows; row++) {

    # go through all columns
    for (col = 0; col < ncols; col++) {
    
      # figure out moving window, clip to region if necessary
      mwrow1 = row - (int)radius
      mwrow2 = mwrow1 + window_size
      if (mwrow1 < 0)
        mwrow1 = 0
      if (mwrow2 > nrows)
        mwrow2 = nrows

      
      mwcol1 = col - (int)radius
      mwcol2 = mwcol1 + window_size
      if (mwcol1 < 0)
        mwcol1 = 0
      if (mwcol2 > ncols)
        mwcol2 = ncols

      # read band values from rastin for row, col

      sum_of_weights = 0
      # go through the window
      for (mwrow = mwrow1; mwrow < mwrow2; mwrow++) {
        for (mwcol = mwcol1; mwcol < mwcol2; mwcol++) {
        
          # check spatial distance
          if the spatial distance to focus pixel is <= radius {

            # set the weight, here 1, can be changed to gaussian weight
            weight = 1
          
            # read band values from rastin for mrow, mcol as oldvalue

            # check spectral distance
            if the spectral distance to focus pixel is <= range_bandwidth {

              # adjust the weight, here 1, can be changed to gaussian weight
              # the weight is a function of the spectral distance, one weight for all bands
              weight *= 1
              sum_of_weights += weight
              for all bands
                newvalue = oldvalue * weight
            }
          }
        }
      }
      # end of moving window

      newvalue /= sum_of_weights

      write newvalue to rastout

      if the spectral distance between old and new value is > threshold
        # not yet converged
        n_changes++

    }
    # end columns
  }
  # end rows
}
# end iterations