1. 引入

上一章我们研究了Focal类中最基础的游标(Cursor)类算子,游标类算子的核心思想代表了大多数Focal类算子的实现方法.

当然,通用的往往更慢,Geotrellis中还有其他三类为特定使用场景优化的Focal类算子:

• 核(kernel)算子:为卷积运算优化的游标类算子
• 像元尺度(cellwise)类算子:为简单邻域优化的游标类算子
• 表面点(surfacepoint)类算子:为山体运算优化的Focal类算子

我们先从核算子入手.

2. 核算子-以Convolve算子为例

核算子是为卷积运算专门优化的游标算子,因此在定义上也与游标算子十分类似:

代码位于[geotrellis/raster/mapalgebra/focal/FocalCalculation.scala]

abstract class KernelCalculation[T](tile: Tile, kernel: Kernel, val analysisArea: Option[GridBounds[Int]], target: TargetCell)
    extends FocalCalculation[T](tile, kernel, analysisArea, target)
{
  // 默认使用线形策略(基于实验结果,线形策略的在卷积运算中表现更好)
  def traversalStrategy = ScanLineTraversalStrategy

  def execute(): T = {
    // 创建了一个核游标对象
    val cursor = new KernelCursor(tile, kernel, bounds)
    val calcFunc: () => Unit =
      target match {
        // ... 省略与游标算子类似的逻辑
      }
    CursorStrategy.execute(cursor, calcFunc, bounds, traversalStrategy)
    result
  }

  def calc(r: Tile, kernel: KernelCursor): Unit
}

对比普通游标类算子,核算子有两个主要区别:

• 游标算子传入了一个邻域对象(neightborhood),而核算子传入了一个卷积核(kernel)对象
• 游标算子根据邻域对象创建了一个游标对象(Cursor),核算子根据卷积核对象创建了一个核游标对象(KernelCursor)

我们先从卷积核开始分析.

2.1 卷积核

卷积核是一种特殊的邻域:

代码位于[geotrellis/raster/mapalgebra/focal/Kernel.scala]

case class Kernel(tile: Tile) extends Neighborhood {
  val extent = (tile.rows / 2).toInt
  val hasMask = false
  def cellType = tile.cellType
}

与普通邻域不同,卷积核有一个表示"权重"的tile对象(权重矩阵),卷积核的尺寸定义也来自于该权重矩阵.

在Geotrellis中定义权重矩阵有3种方式:

• 通过邻域定义
• 定义高斯核(径向基)
• 定义圆形核

我们以通过邻域定义权重矩阵为例:

代码位于[geotrellis/raster/mapalgebra/focal/Kernel.scala]

object Kernel {
  // 隐式转换,所有需要kernel的地方如果传入的是tile,自动转换为kernel
  implicit def tile2Kernel(r: Tile): Kernel = Kernel(r)

  def apply(nbhd: Neighborhood): Kernel = {
    // 因为没有遮罩的正方形邻域将产生全为1的权重矩阵,没有实际意义
    require(nbhd.hasMask, "Neighborhood must have a mask method")
    val w = 2 * nbhd.extent + 1
    val size = w * w
    // tile是使用一维矩阵表示二维矩阵
    val tile = IntArrayTile(new Array[Int](size), w, w)
    // 遍历该邻域的每一个位置
    var r = 0
    while (r < w) {
      var c = 0
      while (c < w) {
        // 有值为1,无值为0
        tile.set(c, r, if (nbhd.mask(c, r)) 0 else 1)
        c = c + 1
      }
      r = r + 1
    }
    new Kernel(tile)
  }
  // ...省略其他两种定义权重矩阵
}

通过邻域定义的方式将会产生这样一种权重矩阵:权重只有1和0之分,原始遮罩中参与计算的区域将设为1,不参与计算的部分将设为0.

我们接下来再看一下核游标是如何操作卷积核的.

2.2 核游标

核游标是特殊的游标:

代码位于[geotrellis/raster/mapalgebra/focal/KernelCursor.scala]

class KernelCursor(r: Tile, kernel: Kernel, analysisArea: GridBounds[Int])
    extends Cursor(r, analysisArea, kernel.extent)
    with MacroIterableTile // 注意这里混入了宏特质
    with Serializable {
  private val ktileArr = kernel.tile.toArray
  private val kcols = kernel.tile.cols

  def foreachWithWeight(f: (Int, Int, Int) => Unit): Unit =
    macro TileMacros.intForeach_impl

  // 实现宏特质约定的foreachIntVisitor方法
  def foreachIntVisitor(f: IntTileVisitor): Unit = {
    var y = rowmin
    var x = 0
    while(y <= rowmax) {
      x = colmin
      while(x <= colmax) {
        // 计算偏移量以返回对应的权重值
        val kcol = focusCol + extent - x
        val krow = focusRow + extent - y
        val w = ktileArr(krow * kcols + kcol)
        f(x, y, w)
        x += 1
      }
      y += 1
    }
  }

  // ...省略全部处理Double类型数据的类似逻辑
}

核游标继承了游标对移动(move)的逻辑,但由于具有特殊的权重矩阵,因此独立实现了专用的遍历方法foreachWithWeight,而不再使用游标中的foreachAdded等方法.

值得注意的是,在核游标中,似乎"多此一举"的使用了宏实现了foreachWithWeight方法,类似的逻辑也可以在游标定义中看见:

代码位于[geotrellis/raster/mapalgebra/focal/Cursor.scala]

// ... 省略
val addedCells = new CellSet {
    def foreach(f: (Int, Int)=>Unit) = { Cursor.this.foreachAdded(f) }
}
// ... 省略

// 实际调用时:
cursor.addedCells.foreach { /* ...省略... */ }

为何不能直接实现foreach方法,非要包上一层呢?因为诸如foreach方法中包含了对private级属性的操作,因为scala严格的访问隔离,因此不能直接在外部调用.

想在外部获取/操作这些private级的字段需要使用闭包的方式将对内部引用的操作包装,如游标将foreach操作包装在一个轻量的内部匿名CellSet类中,而核游标则通过宏将foreach操作包装在匿名IntTileVisitor类中.

2.3 Convolve算子的实现

卷积算子本质上也是一个游标类算子,与其他普通游标类算子的差异就是如何遍历其中的聚焦值了:

代码位于[geotrellis/raster/mapalgebra/focal/Convolve.scala]

object Convolve {
  def calculation(tile: Tile, kernel: Kernel, bounds: Option[GridBounds[Int]] = None, target: TargetCell = TargetCell.All): KernelCalculation[Tile] = {
    if (tile.cellType.isFloatingPoint) {
      if(kernel.cellType.isFloatingPoint) {
        new KernelCalculation[Tile](tile, kernel, bounds, target) with ArrayTileResult {
          def calc(t: Tile, cursor: KernelCursor) = {
            var s = Double.NaN
            // 遍历每一个位置和对应位置的权重值
            cursor.foreachWithWeightDouble { (x, y, w) =>
              val v = t.getDouble(x, y)
              // 最终结果=值*权重
              if(isData(v)) {
                if(isData(s)) { s = s + (v * w) }
                else { s = v * w }
              }
            }
            resultTile.setDouble(cursor.col, cursor.row, s)
          }
        }
      } else {
              // ... 省略处理卷积核为Int类型的类似逻辑
      }
    } else {
        // ... 省略处理源数据为Int类型的类似逻辑
    }
}

至此,核算子/卷积算子的实现就比较清晰了.

3. 像元尺度类算子-以Sum算子为例

像元尺度类算子是对游标类算子在简单邻域场景下的简化,当某个算子仅需要获取增加/移除的像元位置即可实现算法时,就能用像元尺度算子替代.

相比使用标准的游标类算子,它们实现的功能是一样的,但速度会更快,最多可能有40%左右的速度优势.因此,部分统计类的游标类算子也具有像元尺度类算子的实现.

3.1 抽象像元尺度类算子的定义

代码位于[geotrellis/raster/mapalgebra/focal/FocalCalculation.scala]

abstract class CellwiseCalculation[T] (
    r: Tile, n: Neighborhood, analysisArea: Option[GridBounds[Int]], target: TargetCell)
  extends FocalCalculation[T](r, n, analysisArea, target)
{
  // 屏蔽掉默认的三种滑动策略 
  def traversalStrategy: Option[TraversalStrategy] = None

  def execute(): T = n match {
    case s: Square =>
      val calcSetValue: (Int, Int, Int, Int) => Unit =
        target match {
          case TargetCell.All =>
            { (_, _, x, y) => setValue(x, y) }
          // ... 省略 ...
        }
      
      // 包装成对象,作为参数传入CellwiseStrategy方法中
      val strategyCalc =
        new CellwiseStrategyCalculation {
          def add(r: Tile, x: Int, y: Int): Unit = CellwiseCalculation.this.add(r, x, y)
          def remove(r: Tile, x: Int, y: Int): Unit = CellwiseCalculation.this.remove(r, x, y)
          def reset(): Unit = CellwiseCalculation.this.reset()
          def setValue(focusCol: Int, focusRow: Int, x: Int, y: Int): Unit =
            calcSetValue(focusCol, focusRow, x, y)
        }
      // 执行滑动操作
      CellwiseStrategy.execute(r, s, strategyCalc, bounds)
      result
    case _ => sys.error("Cannot use cellwise calculation with this traversal strategy.")
  }

  def add(r: Tile, x: Int, y: Int): Unit
  def remove(r: Tile, x: Int, y: Int): Unit
  def reset(): Unit
  def setValue(x: Int, y: Int): Unit
}

trait CellwiseStrategyCalculation {
  def add(r: Tile, x: Int, y: Int): Unit
  def remove(r: Tile, x: Int, y: Int): Unit
  def reset(): Unit
  def setValue(focusCol: Int, focusRow: Int, x: Int, y: Int): Unit
}

对比游标类算子,像元尺度类算子最主要的区别是不再有游标的概念,而是把游标的操作逻辑和滑动策略融合在一起.

因为限定邻域为最简单的正方形邻域,因此像元尺度类算子主要的逻辑都在其滑动策略中.

3.2 像元尺度滑动策略的定义

在游标类算子中,游标策略负责处理游标的移动,哪些像元位置被添加/移除则由游标内部负责维护.而在像元尺度类算子中,游标策略身兼两职,像元位置添加/移除的逻辑被包含在移动策略的内部:

代码位于[geotrellis/raster/mapalgebra/focal/FocalStrategy.scala]

object CellwiseStrategy {
  def execute(r: Tile, n: Square, calc: CellwiseStrategyCalculation, analysisArea: GridBounds[Int]): Unit =
    handleScanLine(r, n.extent, calc, analysisArea)

  // 定制了专用的线形轨迹滑动逻辑
  private def handleScanLine(r: Tile, n: Int, calc: CellwiseStrategyCalculation, analysisArea: GridBounds[Int]) = {
    // 固定值初始化过程
    
    // 行列移动的上下限
    val rowMin = analysisArea.rowMin
    val colMin = analysisArea.colMin
    val rowMax = analysisArea.rowMax
    val colMax = analysisArea.colMax
    // 实际数据的行列上下限,如果不限制分析区域,则该值与rowMax/colMax一致
    val rowBorderMax = r.rows - 1
    val colBorderMax = r.cols - 1
    // 行列号偏移,如果不限制分析区域,则该值与rowMin/colMin一致
    val analysisOffsetCols = analysisArea.colMin
    val analysisOffsetRows = analysisArea.rowMin
    
    // 随移动变化的值
    
    // 当前聚焦的行
    var focusRow = rowMin
    while (focusRow <= rowMax) {
      // 当前滑动窗口的行范围
      val curRowMin = max(0, focusRow - n)
      val curRowMax = min(rowBorderMax, focusRow + n )
      // 重置计算核心
      calc.reset()
      // 当前滑动窗口的列范围
      val curColMin = max(0, colMin - n)
      val curColMax = min(colBorderMax, colMin + n)
      // 遍历当前滑动窗口下的每一个位置,并调用add方法
      var curRow = curRowMin
      while (curRow <= curRowMax) {
        var curCol = curColMin
        while (curCol <= curColMax) {
          calc.add(r, curCol, curRow)
          curCol += 1
        }
        curRow += 1
      }
      // 给第一列赋值
      calc.setValue(colMin, focusRow, 0, focusRow - rowMin)
      // 处理其他列
      var focusCol = colMin + 1
      // 处理滑动中移除的列
      while (focusCol <= colMax) {
        val oldWestCol = focusCol - n - 1
        if (oldWestCol >= 0) {
          var yy = curRowMin
          while (yy <= curRowMax) {
            calc.remove(r, oldWestCol, yy)
            yy += 1
          }
        }
        // 处理滑动中添加的列
        val newEastCol = focusCol + n
        if (newEastCol <= colBorderMax) {
            var yy = curRowMin
            while (yy <= curRowMax) {
              calc.add(r, newEastCol, yy)
              yy += 1
            }
          }
        // 给当前行剩余的每一列赋值
        calc.setValue(focusCol, focusRow, focusCol - colMin, focusRow - rowMin)
        focusCol += 1
      }
      // 再处理下一行
      focusRow += 1
    }
  }
}

像元尺度策略是线形策略,具体描述如下:

1. 开始处理第一行.先处理第一列,因为第一列只有add逻辑
2. 再处理剩余且他列,因为非首列就既有add逻辑,又有remove逻辑
3. 切换到下一行的第一列.以此类推直到处理完所有行

在这里还有一点需要考虑,如果同样使用正方形邻域(最简单邻域),针对同样的算子(比如Sum算子),为何依旧是像元尺度类算子比游标类的实现更快呢?

其实原理比较简单:因为游标类算子为了兼容各种复杂算子,内部使用了更多的逻辑判断和临时变量运算,这些逻辑的消耗虽然很小,但与之相比,使用更少判断和变量运算的像元尺度算子,确实会一直保持对游标类算子的优势,尽管这个优势是线性而非指数级.

所以,针对同一个算子,如果既有游标类实现,又有像元尺度类实现,在使用简单邻域的情况下应优先使用像元尺度类实现.

3.3 Sum算子的实现

Sum算子同时实现了游标算子和像元尺度算子:

代码位于[geotrellis/raster/mapalgebra/focal/Sum.scala]

object Sum {
  def calculation(tile: Tile, n: Neighborhood, bounds: Option[GridBounds[Int]] = None, target: TargetCell = TargetCell.All): FocalCalculation[Tile] = {
    if (tile.cellType.isFloatingPoint) n match {
      // 在正方形邻域时会主动切换到像元尺度算子
      case Square(ext) => new CellwiseDoubleSumCalc(tile, n, bounds, target)
      // 复杂邻域则切换到游标算子
      case _ =>           new CursorDoubleSumCalc(tile, n, bounds, target)
    } else n match {
      // ... 省略处理Int类型的类似逻辑
    }
  }

  def apply(tile: Tile, n: Neighborhood, bounds: Option[GridBounds[Int]] = None, target: TargetCell = TargetCell.All): Tile =
    calculation(tile, n, bounds, target).execute
}

class CellwiseDoubleSumCalc(r: Tile, n: Neighborhood, bounds: Option[GridBounds[Int]], target: TargetCell)
    extends CellwiseCalculation[Tile](r, n, bounds, target) with ArrayTileResult {

  // 计数器
  var count: Int = 0
  var total: Double = Double.NaN

  def add(r: Tile, x: Int, y: Int) = {
    val v = r.getDouble(x, y)
    if (isData(v)) {
      if (isData(total)) { total += v; count += 1 }
      else { total = v; count = 1 }
    }
  }

  def remove(r: Tile, x: Int, y: Int) = {
    val v = r.getDouble(x, y)
    if (isData(v)) {
      count -= 1
      // 保证在全部移除时设置为NaN而不是浮点计算误差值
      if (count == 0) total = Double.NaN else total -= v
    }
  }

  def reset() = total = Double.NaN
  def setValue(x: Int, y: Int) = resultTile.setDouble(x, y, total)
}

可见,一个标准的像元尺度类算子基本上通过定义add和remove方法来实现特定的算法.当然,类似Max算法,就没法仅仅通过这两个方法实现目标算法,因此无法使用像元尺度类算子实现.

至此,像元尺度类算子的功能与实现也比较清晰了,与核算子一样,都是特定场景下对性能的优化.

4. 表面点类算子-以Aspect算子为例

对表面点类相关算法的解释在Arcgis官网中有详细的图文描述.

在表面点类算子中,已经很难看到游标类算子的影子了,因为其中具有相当多与山体计算相关的针对性方法,但我们依然可以抓住游标类算子的一般规律--关注滑动窗口/算子核心/移动策略--来对比分析表面点类算子.

我们再来回顾一下游标类算子,核算子与像元尺度类算子:

算子类型	滑动窗口	算子核心	滑动策略
游标类算子	任意形状任意尺寸的邻域	游标	之字形策略
核算子	任意尺寸的正方形邻域	核游标	线形策略
像元尺度类算子	任意尺寸的正方形邻域	无(融入策略中)	自定义的线形策略

我们再从这三个角度分析表面点类算子.

4.1 抽象表面点类算子的定义

表面点类算子更接近于像元尺度算子:

• 移动窗口为3x3尺寸的正方形邻域(自身算法性质决定)
• 没有算子核心(融入了移动策略中)
• 使用自定义的线形策略作为移动策略

不过虽然都是自定义的线形移动策略,表面点类算子的实现还是与像元尺度类算子的有很大的不同.表面点类算子将整个Tile分为9个大区,分别进行对应的处理:

代码位于[geotrellis/raster/mapalgebra/focal/hillshade/SurfacePointCalculation.scala]

abstract class SurfacePointCalculation[T](r: Tile, n: Neighborhood, analysisArea: Option[GridBounds[Int]], val cellSize: CellSize, target: TargetCell = TargetCell.All)
  extends FocalCalculation[T](r, n, analysisArea, target)
{
  // 使用特质以避免装箱造成的性能损失
  private trait SetValue { def apply(z: Double, x: Int, y: Int): Unit }

  var lastY = -1

  var cellWidth = 0.0
  var cellHeight = 0.0
  
  // 表面点类算子的滑动窗口,定义了左中右三列的3x3矩阵
  var west = new Array[Double](3)
  var base = new Array[Double](3)
  var east = new Array[Double](3)
  
  // 表面点对象,用于记录某个位置的各种表面点参数
  val s = new SurfacePoint

  // 实际算子需要实现这个函数
  def setValue(x: Int, y: Int, s: SurfacePoint): Unit

  private val _setValue: SetValue =
    target match {
      case TargetCell.All =>
        new SetValue {
          def apply(z: Double, x: Int, y: Int): Unit = {
            // 计算表面点相关参数值
            calcSurface()
            // 调用实际算子,实现赋值给结果集对应位置
            setValue(x, y, s)
          }
        }
      // ... 省略...
    }
  
  // 因为设计为线形移动策略,因此仅需向右移动
  // 向右移动的本质是原先的右侧变为中侧/中侧变为左侧,左侧变为右侧
  // 因为变为右侧的左侧列马上会被新值覆盖,因此不需要清空值
  def moveRight() = {
    val tmp = west
    west = base
    base = east
    east = tmp
  }
  
  // 计算表面点值
  protected def calcSurface(): Unit = {
    // base(1)就是3x3滑动窗口中最中心位置的值
    // 如果该值不存在,则当前位置的表面点无意义,无论周围是否有值
    if(isNoData(base(1))) {
      s.`dz/dx` = Double.NaN
      s.`dz/dy` = Double.NaN
    } else {
      // 计算中心点周围8位置的值,如果值为Nodata,则用中心点值替代
      // 取值策略也是表面点类算子与像元尺度类算子的不同之处
      val neValue = if(isData(east(0))) east(0) else base(1)
      val eValue = if(isData(east(1))) east(1) else base(1)
      val seValue = if(isData(east(2))) east(2) else base(1)
      val nValue = if(isData(base(0))) base(0) else base(1)
      val sValue = if(isData(base(2))) base(2) else base(1)
      val nwValue = if(isData(west(0))) west(0) else base(1)
      val wValue = if(isData(west(1))) west(1) else base(1)
      val swValue = if(isData(west(2))) west(2) else base(1)
      // 给表面点的两个偏导数赋值,它们是相关山体算法的关键
      s.`dz/dx` = (neValue + 2*eValue + seValue - nwValue - 2*wValue - swValue) / (8 * cellWidth)
      s.`dz/dy` = (swValue + 2*sValue + seValue - nwValue - 2*nValue - neValue) / (8 * cellHeight)
    }
  }

  def execute(): T = {
    val colMin = bounds.colMin
    val colMax = bounds.colMax

    val rowMin = bounds.rowMin
    val rowMax = bounds.rowMax

    val colBorderMax = r.cols - 1
    val rowBorderMax = r.rows - 1
    
    // 设置像元实际尺寸
    cellWidth = cellSize.width
    cellHeight = cellSize.height

    def getValSafe(col: Int, row: Int, focalVal: Double) = {
      if(col < 0 || colBorderMax < col || row < 0 || rowBorderMax < row) {
        focalVal
      } else {
        r.getDouble(col, row)
      }
    }

    var focalValue = r.getDouble(colMin, rowMin)

    /// 左上
    // 如果越界,则用当前聚焦点的值替代
    west(0) = getValSafe(colMin-1, rowMin-1, focalValue)
    west(1) = getValSafe(colMin-1, rowMin  , focalValue)
    west(2) = getValSafe(colMin-1, rowMin+1, focalValue)
    base(0) = getValSafe(colMin  , rowMin-1, focalValue)
    base(1) = focalValue
    base(2) = r.getDouble(colMin, rowMin + 1)
    east(0) = getValSafe(colMin+1, rowMin-1, focalValue)
    east(1) = r.getDouble(colMin + 1, rowMin)
    east(2) = r.getDouble(colMin + 1, rowMin + 1)
    _setValue(focalValue, 0, 0)

    var col = colMin + 1

    /// 中上
    while (col < colMax) {
      moveRight()
      focalValue = r.getDouble(col, rowMin)
      west(0) = getValSafe(col-1, rowMin-1, focalValue)
      base(0) = getValSafe(col  , rowMin-1, focalValue)
      east(0) = getValSafe(col+1, rowMin-1, focalValue)
      east(1) = r.getDouble(col+1, rowMin)
      east(2) = r.getDouble(col+1, rowMin + 1)
      _setValue(focalValue, col-colMin, 0)
      col += 1
    }

    /// 右上
    moveRight()
    focalValue = r.getDouble(col, rowMin)
    west(0) = getValSafe(col-1, rowMin-1, focalValue)
    base(0) = getValSafe(col  , rowMin-1, focalValue)
    east(0) = getValSafe(col+1, rowMin-1, focalValue)
    east(1) = getValSafe(col+1, rowMin  , focalValue)
    east(2) = getValSafe(col+1, rowMin+1, focalValue)
    _setValue(focalValue, col-colMin, 0)

    var row = rowMin + 1
    
    // 省略处理左中,中部,右中,左下,中下,右下的类似逻辑

    result
  }
}

在代码中需要操作表面点对象(SurfacePoint),为其两个字段赋值.

不过表面点对象并非是类似游标的东西,它是一个1x1的点,记录了某个位置山体算法所需要的必要参数:

代码位于[geotrellis/raster/mapalgebra/focal/hillshade/SurfacePointCalculation.scala]

class SurfacePoint() {
  var isNaN = false
  
  // 为保持语义,使用了字面符标识作为变量名
  var `dz/dx` = Double.NaN
  var `dz/dy` = Double.NaN
  
  // 角度坡向
  def aspectAzimuth() = {
    if (Math.abs(`dz/dx`) <= 1E-10 && Math.abs(`dz/dy`) <= 1E-10) {
      -1.0
    } else {
      var aAzimuth = toDegrees(atan2(`dz/dy`, `dz/dx`) - Pi * 0.5)
      if (aAzimuth < 0) { aAzimuth += 360.0 }
      if (Math.abs(aAzimuth - 360.0) <= 1E-10) { aAzimuth = 0.0 }
      aAzimuth
    }
  }
    
  // 弧度坡向
  def aspect() = {
    var a = atan2(`dz/dy`, -`dz/dx`)
    if (`dz/dx` == 0 && `dz/dy` == 0) {
      a = Double.NaN
    } else {
      if (a < 0) { a += 2*Pi }
    }
    if (a == 2*Pi) { a = 0.0 }
    a
  }
  
  // 坡度
  def slope(zFactor: Double): Double = atan(zFactor * sqrt(`dz/dx` * `dz/dx` + `dz/dy` * `dz/dy`))
  def slope(): Double = slope(1.0)

  // ... 省略坡度/坡向的相关三角函数
}

有了表面点对象,我们就可以在实际的算子中调用它了.

4.2 Aspect算子的实现

代码位于[geotrellis/raster/mapalgebra/focal/aspect.scala]

object Aspect {

  def apply(tile: Tile, n: Neighborhood, bounds: Option[GridBounds[Int]], cs: CellSize, target: TargetCell = TargetCell.All): Tile = {
    new SurfacePointCalculation[Tile](tile, n, bounds, cs, target)
      with DoubleArrayTileResult
    {
      def setValue(x: Int, y: Int, s: SurfacePoint): Unit = {
        resultTile.setDouble(x, y, s.aspectAzimuth)
      }
    }
  }.execute()
}

从坡向算子可以看出,因为具有固定的尺寸,表面点类算子的实现较为简单,操作的表面点对象是单一点,从而无需使用其他Focal类算子中必须通过诸如foreach/foreachAdded方法求值的设计.

5. 总结

我们从通用性最强的游标类算子开始分析,得到了Focal类算子通用的三大核心点:

• 滑动窗口:描述Focal操作每次聚焦的是哪片区域
• 滑动策略:描述如何在Tile中滑动窗口
• 算子核心:描述如何从滑动窗口中获得我们想要的结果值

基于这3点,我们分析了特化了滑动窗口和算子核心的核算子,简化了滑动策略和算子核心的像元尺度算子,和特化滑动窗口且简化滑动策略和算子核心的表面点算子,它们都是针对特定领域应用而进行的性能优化.

至此,我们就分析完了4大类Focal算子的结构与实现.