1. 引入

上一章我們研究了Geotrellis中Local算子的實現.Local類算子種類較多,實現也較簡單,其核心為:Tile中某個位置的值由另一個Tile中相同位置的值通過計算得到.

Geotrellis中也存在更加復雜的Focal類算子,其核心為:Tile中某個像元位置的值由該像元周圍若干像元(鄰域)的值經過特定計算得到.可見,Focal類算子有兩個關鍵點:

• 定義鄰域
• 定義如何根據鄰域內容計算結果

我們首先看鄰域是如何定義的.

2. 鄰域定義

在Geotrellis中,鄰域描述了某個像元周圍哪些位置的像元參與到計算中.這些周圍區(qū)域的形狀可能簡單也可能復雜,為了不針對各種形狀制定復雜的結構從而影響計算效率,Geotrellis通過這樣一種方式描述鄰域:

代碼位于[geotrellis/raster/mapalgebra/focal/Neighborhood.scala]

trait Neighborhood extends Serializable {
  // 聚焦范圍,范圍為1則對應3x3的聚焦鄰域,即向中心點周圍8方向各擴展1個單位
  val extent:Int 

  // 是否具有遮罩(只有最基礎的正方形鄰域沒有遮罩)
  val hasMask:Boolean

  /** 
   * 定義具體的鄰域:
   * 通過定義不同位置的返回值確定鄰域形狀,返回true代表被遮罩覆蓋,不參與計算
   * (0,0) 代表左上角, (extent*2+1,extent*2+1) 代表右下角
   */
  def mask(col:Int,row:Int):Boolean = { false }
}

Geotrellis定義了5種鄰域(#符號表示在鄰域中最終參與計算,.符號表示不參與計算,extent設為4):

正方形(Square)鄰域      十字形(Nesw)鄰域        圓形(Circle)鄰域

# # # # # # # # #       . . . . # . . . .       . . . . # . . . .
# # # # # # # # #       . . . . # . . . .       . . # # # # # . . 
# # # # # # # # #       . . . . # . . . .       . # # # # # # # . 
# # # # # # # # #       . . . . # . . . .       . # # # # # # # .
# # # # # # # # #       # # # # # # # # #       # # # # # # # # #
# # # # # # # # #       . . . . # . . . .       . # # # # # # # .
# # # # # # # # #       . . . . # . . . .       . . # # # # # . .
# # # # # # # # #       . . . . # . . . .       . . # # # # # . .
# # # # # # # # #       . . . . # . . . .       . . . . # . . . . 

環(huán)形(Annulus)鄰域       90°-180°楔形(Wedge)鄰域   180°-90°的楔形鄰域

. . . . # . . . .       . . . . # . . . .       . . . . # . . . . 
. . # # # # # . .       . . # # # . . . .       . . . . # # # . . 
. # # # # # # # .       . # # # # . . . .       . . . . # # # # . 
. # # . . . # # .       . # # # # . . . .       . . . . # # # # . 
# # # . . . # # #       # # # # # . . . .       # # # # # # # # # 
. # # . . . # # .       . . . . . . . . .       . # # # # # # # . 
. # # # # # # # .       . . . . . . . . .       . # # # # # # # . 
. . # # # # # . .       . . . . . . . . .       . . # # # # # . . 
. . . . # . . . .       . . . . . . . . .       . . . . # . . . . 

了解完鄰域的定義,我們就需要看看Focal類算子如何根據鄰域中的內容計算結果.

3. Focal類算子的整體架構

Focal算子可以劃分為4大類:

• CursorCalculation(游標類算子):
  • Max(最大值)/Min(最小值)
  • StandardDeviation(標準差)
  • Mean(均值)/Median(中值)/Mode(眾數)/sum(總和)/moran(莫蘭指數)
• KernelCalculation(核算子)
  • Convolve(卷積)
• CellwiseCalculation(像元尺度算子)
  • Conway(生命游戲)
  • Mean(均值)/Median(中值)/Mode(眾數)/sum(總和)/moran(莫蘭指數)
• SurfacePointCalculation(表面點算子)
  • aspect(坡向)
  • slope(坡度)
  • hillshade(山體陰影)

其中部分統計算子同時具有游標類算子和像元尺度算子兩種實現.

這4大類算子都繼承自抽象的Focal算子類(FocalCalculation):

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

trait Resulting[T] {
  val copyOriginalValue: (Int, Int, Int, Int) => Unit
  def result: T
}

abstract class FocalCalculation[T](
    val r: Tile, n: Neighborhood, analysisArea: Option[GridBounds[Int]], val target: TargetCell)
  extends Resulting[T]
{
  // 計算范圍.默認為Tile自身的范圍 
  val bounds: GridBounds[Int] = analysisArea.getOrElse(GridBounds(r))

  def execute(): T
}

其中包含了全部Focal算子需要具備的三個核心:

• 輸入:定義了支持的輸入參數
• 運算:約定實現execute方法
• 輸出:約定實現對結果集操作的方法
  • 不同的算子按需混入不同類型的Resulting特質,實現操作具體數據類型的結果集
  • 我們以最具有一般性的ArrayTileResult為例(ArrayTileResult適用于任何不需要改變原始數據類型的算子):

// 通過自類型注釋(self-type annotation),限定自己必須是FocalCalculation類型
// 保證其中具有r,bounds等對象
trait ArrayTileResult extends Resulting[Tile] { self: FocalCalculation[Tile] =>
  def resultCellType: DataType with NoDataHandling = r.cellType
  val cols: Int = bounds.width
  val rows: Int = bounds.height
  // 定義返回的對象,類型為傳入tile對象的數據類型
  val resultTile: MutableArrayTile = ArrayTile.empty(resultCellType, cols, rows)
  val copyOriginalValue: (Int, Int, Int, Int) => Unit =
    if(!r.cellType.isFloatingPoint) {
      { (focusCol: Int, focusRow: Int, col: Int, row: Int) =>
        resultTile.set(col, row, r.get(focusCol, focusRow))
      }
    } else {
      { (focusCol: Int, focusRow: Int, col: Int, row: Int) =>
        resultTile.setDouble(col, row, r.getDouble(focusCol, focusRow))
      }
    }
  // 將結果集指向約定值,在算子中我們操作的也是這個值
  // 使用def預定義再通過后期指向的方式是為了通過編譯
  def result = resultTile
}

4大類中最關鍵的算子是游標類算子(CursorCalculation),核算子是其特化,而像元尺度算子和表面點算子可以看做其簡化.因此我們從游標類算子開始分析.

4. 游標類算子-以Max算子為例

4.1 抽象游標類算子的定義

首先來分析抽象的游標類算子是如何定義的:

代碼位于[geotrellis/raster/mapalgebra/focal/FocalMethods.scala]

abstract class CursorCalculation[T](tile: Tile, n: Neighborhood, val analysisArea: Option[GridBounds[Int]], target: TargetCell)
  extends FocalCalculation[T](tile, n, analysisArea, target)
{
  def traversalStrategy = TraversalStrategy.DEFAULT

  def execute(): T = {
    // 1.
    val cursor = Cursor(tile, n, bounds)
    // 2.
    val calcFunc: () => Unit =
      target match {
        case TargetCell.All =>
          { () => calc(tile, cursor) }
        case TargetCell.Data =>
          { () =>
            calc(tile, cursor)
            if(!isData(r.get(cursor.focusCol, cursor.focusRow))){
              copyOriginalValue(cursor.focusCol, cursor.focusRow, cursor.col, cursor.row)
            }
          }
        case TargetCell.NoData =>
          { () =>
            calc(tile, cursor)
            if(!isNoData(r.get(cursor.focusCol, cursor.focusRow))) {
              copyOriginalValue(cursor.focusCol, cursor.focusRow, cursor.col, cursor.row)
            }
          }
      }
    // 3.
    CursorStrategy.execute(cursor, calcFunc, bounds, traversalStrategy)
    result
  }

  def calc(tile: Tile, cursor: Cursor): Unit
}

游標算子的核心是execute方法,在該方法中主要執(zhí)行了如下3步:

1. 定義一個游標對象(Cursor)
2. 定義calcFunc方法
3. 將一系列參數傳入CursorStrategy.execute方法中執(zhí)行

我們按順序先從游標對象的定義開始分析.

4.2 游標對象的定義

Focal類算子是一種滑動窗口型算子,而滑動的那個"窗口"就是游標對象,它定義了這個窗口的大小與形狀,是我們進行滑動和讀取數據的基本粒度.

4.2.1 游標對象的初始化與滑動

游標對象進行滑動前首先需要將其初始化:

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

class Cursor(r: Tile, analysisArea: GridBounds[Int], val extent: Int) {
 
 // 記錄總的行列數
 private val rows = r.rows
 private val cols = r.cols

 // 偏移量,如果不指定限制范圍,則偏移量為0
 val analysisOffsetCols = analysisArea.colMin
 val analysisOffsetRows = analysisArea.rowMin
 
 // 游標的邊長
 private val d = 2 * extent + 1

 // 記錄當前游標所在區(qū)域的行列范圍
 private var _colmin = 0
 private var _colmax = 0
 private var _rowmin = 0
 private var _rowmax = 0

 // 相當于對上面的min/max定義了一個公開的"get"方法
 protected def colmin = _colmin
 protected def colmax = _colmax
 protected def rowmin = _rowmin
 protected def rowmax = _rowmax

 // 記錄在當前移動過程中,加入的/移除的行列號
 private var addedCol = 0
 private var removedCol = 0
 private var addedRow = 0
 private var removedRow = 0
   
 // 當前移動方向
 var movement = NoMovement

 // 當前聚焦位置及其"get"方法
 private var _col = 0
 private var _row = 0
 def focusCol = _col
 def focusRow = _row

 def isReset = movement == NoMovement

 // 當前聚焦的行列號的真實值
 def col = _col - analysisOffsetCols
 def row = _row - analysisOffsetRows

 // 設置起點
 def centerOn(col: Int, row: Int) = {
   movement = NoMovement
   // 設置當前要處理的行列號為起點行列號
   _col = col
   _row = row
   // 計算在該起點下,最大最小行列號都是多少
   _colmin = max(0, _col - extent)
   _colmax = min(cols - 1, _col + extent)
   _rowmin = max(0, _row - extent)
   _rowmax = min(rows - 1, _row + extent)
 }
 
 // ... 省略其他
}

初始化過程實現了三個功能:

• 限定行列范圍
• 設定滑動起點
• 定義輔助變量

我們以指定參數[Tile=5x4,extent=1,不指定分析范圍,初始位置為(0,0)]為例:

初始化后各個變量為:

初始化完成后就可以進行滑動(move)操作了:

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

sealed trait Movement { val isVertical: Boolean }
// 定義滑動方向標志位
object Movement {
 val Up = new Movement { val isVertical = true }
 val Down = new Movement { val isVertical = true }
 val Left = new Movement { val isVertical = false }
 val Right = new Movement { val isVertical = false }
 val NoMovement = new Movement { val isVertical = false }
}

class Cursor(r: Tile, analysisArea: GridBounds[Int], val extent: Int) {
   // ... 省略初始化階段
   
   // 滑動游標
   def move(m: Movement) = {
       movement = m
       m match {
         case Up =>
           addedRow = _rowmin - 1
           removedRow = _row + extent
           _row -= 1
         // ... 省略類似的向下移動邏輯
         case Left =>
           addedCol = _colmin - 1
           removedCol = _col + extent
           _col -= 1
         // ... 省略類似的向右移動邏輯
         case _ =>
       }
       
       // 重新計算滑動后的新值
       _colmin = max(0, _col - extent)
       _colmax = min(cols - 1, _col + extent)
       _rowmin = max(0, _row - extent)
       _rowmax = min(rows - 1, _row + extent)
 }
   
   // ...省略其他
 
}

我們以向右滑動為例:

此時對比初始值:

這些變量都是為了方便進行遍歷操作而準備的.

4.2.2 簡單鄰域下的遍歷操作

遍歷操作使我們可以獲取該游標中的每一個像元的值.在上文中我們介紹的幾種鄰域中,最簡單的是正方形鄰域,因為沒有需要跳過的位置,因此處理起來也最簡單.我們以正方形鄰域為例,看看遍歷操作的一般邏輯.

針對游標對象,有兩類遍歷方式:

• 完全遍歷:遍歷游標中的全部位置.
• 變化區(qū)域遍歷:只遍歷上一次滑動操作后新加入/退出的行列,這在某些算子中可以避免重復計算

核心代碼如下:

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

class Cursor(r: Tile, analysisArea: GridBounds[Int], val extent: Int) {
    
    // ... 省略初始化階段
    
    // 遍歷全部位置
    protected def foreach(f: (Int, Int)=>Unit): Unit = {
          var y = _rowmin
          var x = 0
          // 遍歷全部的行列
          while(y <= _rowmax) {
            x = _colmin
            while(x <= _colmax) {
              f(x, y)
              x += 1
            }
            y += 1
          }
        }
    
    // 只遍歷移動后新加入的位置
    protected def foreachAdded(f: (Int, Int)=>Unit): Unit = {
        // 如果不移動,則與foreach無異
        if(movement == NoMovement) {
          foreach(f)
        } else if (movement.isVertical) {
          // 垂直方向移動后,只遍歷新加入的那行
          // 排除無意義的addedRow值
          if(0 <= addedRow && addedRow < rows) {
              var x = _colmin
              while(x <= _colmax) {
                f(x, addedRow)
                x += 1
              }
          }
        } else { // 水平方向移動后,只遍歷新加入的那列
          if(0 <= addedCol && addedCol < cols) {
              var y = _rowmin
              while(y <= _rowmax) {
                f(addedCol, y)
                y += 1
              }
          }
        }
      }
  
    // ... 省略與foreachAdded相似的foreachRemoved
    
    // ... 省略其他

4.2.3 復雜鄰域下的遍歷操作與游標遮罩的定義

游標遮罩(CursorMask)是一個輔助遍歷器,輔助游標遍歷非正方形的復雜鄰域.

相比正方形鄰域,復雜鄰域需要面對的問題就是其中有若干區(qū)域在遍歷中應被遮罩遮蔽,不參與計算.比較笨的方法是每個位置都進行判斷,但必然會產生大量重復計算.為避免浪費計算資源,游標遮罩采取空間換時間的策略,設計了"快查表"機制,盡可能在不占用很多額外空間的前提下只計算需要計算的位置.

在當前語境下,"遮罩"等同于"不參與計算","非遮罩"等同于"參與計算"

4.2.3.1 "快查表"的數據結構

游標遮罩為"快查表"設計了一種便于計算的用一維數組表示二維數據的數據結構:遮罩集合(MaskSet):

代碼位于[geotrellis/raster/mapalgebra/focal/CursorMask.scala]

private class MaskSet {
    // d=2*extent+1,與Cursor中意義一致
    private var data = Array.ofDim[Int](d*(d+1))
    // 插入一個值的時候需要更新兩個值
    def add(i:Int,v:Int) = {
      val len = data(i*(d+1)) + 1
      data(i*(d+1) + len) = v
      data(i*(d+1)) = len
    }
    // 遍歷時需要指定所在行
    def foreachAt(i:Int)(f:Int=>Unit) = {
      val len = data(i*(d+1))
      if(len > 0) {
        var x = 0
        while(x < len) {
          f(data(i*(d+1) + x + 1))
          x += 1
        }
      }
    }
}

遮罩集合雖然使用1維數組存儲數據,但它與游標的尺寸是對應的.一個尺寸為3x3的游標,對應的遮罩集合中data數組的長度就是3x(3+1)=12,即我們可以將其看做一個3x4的二維矩陣疊放到一個一維數組中,額外申請的一列用來存儲每一行中元素的數量.這種數據結構的好處是:在一個數據稀疏的矩陣中,可以避免大量不必要的判斷.

下面我們來分析在游標遮罩對象中,各種"快查表"都是如何被初始化的.

4.2.3.2 "快查表"的初始化

代碼位于[geotrellis/raster/mapalgebra/focal/CursorMask.scala]

class CursorMask(d:Int,f:(Int,Int)=>Boolean) {
    // ...省略MaskSet的定義
    
    // 下文中的最左/最右/位置等都是指在游標范圍內的相對位置
    // 需要通過反算得到在游標移動后在tile中的真實位置
    
    // 記錄所有參與計算的非遮罩位置
    private val unmasked = new MaskSet
    
    // 記錄最左列/最右列哪一行參與計算,因為只有1列,所以無需定義為遮罩集合的形式
    private var westColumnUnmasked = Array.ofDim[Int](d+1)
    private var eastColumnUnmasked = Array.ofDim[Int](d+1)
    
    // 記錄向某個方向移動后改變遮罩性質的全部位置
    private val maskedAfterMoveLeft = new MaskSet
    private val unmaskedAfterMoveLeft = new MaskSet
    private val maskedAfterMoveUp = new MaskSet
    private val unmaskedAfterMoveUp = new MaskSet
    
    var x = 0
    var y = 0
    var len = 0
    var isMasked = false
    var leftIsMasked = false
    // 遍歷游標中的全部位置
    while(y < d) {
        x = 0
        len = 0
        while(x < d) {
          isMasked = f(x,y)
          if(!isMasked) {
            // 記錄所有參與計算的非遮罩位置
            unmasked.add(y,x)
            // 記錄最左/最右兩列的信息,因為他們比較特殊,所以特別記錄
            if(x == 0) {
              val len = westColumnUnmasked(0) + 1
              westColumnUnmasked(len) = y
              westColumnUnmasked(0) = len
            }
            if(x == d - 1) {
              val len = eastColumnUnmasked(0) + 1
              eastColumnUnmasked(len) = y
              eastColumnUnmasked(0) = len
            }
          }
        
          // 記錄隨著向左移動,是否發(fā)生遮罩狀態(tài)變化
          if(x > 0) {
            val leftIsMasked = f(x-1,y)
            // 向左移動將變?yōu)榉钦谡譅顟B(tài)
            if(leftIsMasked && !isMasked) {
              unmaskedAfterMoveLeft.add(y,x)
            } else if(!leftIsMasked && isMasked) {
              maskedAfterMoveLeft.add(y,x)
            }
          }
          // ... 省略類似的向上移動邏輯
          x += 1
        }
        y += 1
    }
    
    // ... 省略其他
}

我們以十字型鄰域(3x3)為例:

 . # .
 # # #
 . # .

初始化后,各個快查表的內容為:

unmasked
[1, 1, null, null, 3, 0, 1, 2, 1, 1, null, null]
westColumnUnmasked
[1, 1, 0, 0]
eastColumnUnmasked
[1, 1, 0, 0]
maskedAfterMoveUp
[null, null, null, null, null, null, null, null, 2, 0, 2, null]
unmaskedAfterMoveUp
[null, null, null, null, 2, 0, 2, null, null, null, null, null]
maskedAfterMoveLeft
[1, 2, null, null, null, null, null, null, 1, 2, null, null]
unmaskedAfterMoveLeft
[1, 1, null, null, null, null, null, null, 1, 1, null, null]

這些快查表中內容現在看比較抽象,我們將在遍歷操作中看看它們如何起作用.

4.2.3.3 使用快查表遍歷數據

與簡單鄰域的遍歷一樣,復雜鄰域也有完全遍歷和變化區(qū)域遍歷兩種.我們先從完全遍歷開始:

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

object Cursor {
  def apply(r: Tile, n: Neighborhood, analysisArea: GridBounds[Int]): Cursor = {
    val result = new Cursor(r, analysisArea, n.extent)
    // 復雜鄰域時設置游標遮罩對象
    if(n.hasMask) { result.setMask(n.mask) }
    result
  }
}

class Cursor(r: Tile, analysisArea: GridBounds[Int], val extent: Int) {

    // ... 省略
    
    def setMask(f: (Int, Int) => Boolean) = {
        hasMask = true
        mask = new CursorMask(d, f)
    }
  
    protected def foreach(f: (Int, Int)=>Unit): Unit = {
        if(!hasMask) {
            // ... 省略簡單鄰域的情況
        } else {
            var y = 0
            while(y < d) {
                // 遍歷某一行中的全部列中參與計算的位置
                mask.foreachX(y) { x =>
                  // 反算回真實位置
                  val xTile = x + (_col - extent)
                  val yTile = y + (_row - extent)
                  // 避免邊緣位置越界
                  if(_colmin <= xTile && xTile <= _colmax && _rowmin <= yTile && yTile <= _rowmax) {
                    f(xTile, yTile)
                  }
                }
                y += 1
            }
        }
    }
    
    // ... 省略其他
}

核心為對游標遮罩對象foreachX方法的調用:

代碼位于[geotrellis/raster/mapalgebra/focal/CursorMask.scala]

class CursorMask(d:Int,f:(Int,Int)=>Boolean) {
    // ... 省略
    
    // 固定行,遍歷列
    def foreachX(row:Int)(f:Int=>Unit) = {
        // 遍歷全部非遮罩位置
        unmasked.foreachAt(row)(f)
    }
    
    // 省略
}

變化區(qū)域遍歷則復雜許多,因為不僅要處理邊緣部分的增減,還要處理因為遮罩區(qū)域變化引起的增減.我們以foreachAdded為例:

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

class Cursor(r: Tile, analysisArea: GridBounds[Int], val extent: Int) {

    // ... 省略

    protected def foreachAdded(f: (Int, Int)=>Unit): Unit = {
        // 先處理邊緣部分
        // 垂直移動的情況
        if (movement.isVertical) {
          if(0 <= addedRow && addedRow < rows) {
            if(!hasMask) {
                // ... 省略簡單鄰域的情況
            } else {
              // addedRow - (_row - extent)=在Cursor中row的相對變化值
              mask.foreachX(addedRow - (_row - extent)) { x =>
                // 再從相對值反算回來
                val xTile = x + (_col - extent)
                if(0 <= xTile && xTile <= cols) {
                  f(xTile, addedRow)
                }
              }
              
              // 等同于下面的方法
              /*
              if(movement == Up) {
                mask.foreachX(0) { x =>
                  val xTile = x + (_col - extent)
                  if(0 <= xTile && xTile < cols) {
                    f(xTile, addedRow)
                  }
                }
              }
              else {
                mask.foreachX(d-1) { x =>
                  val xTile = x + (_col - extent)
                  if(0 <= xTile && xTile < cols) {
                    f(xTile, addedRow)
                  }
                }
              }
              */
            }
          }
        } else { // 水平移動的情況
          if(0 <= addedCol && addedCol < cols) {
            if(!hasMask) {
                // ... 省略簡單鄰域的情況
            } else {
              if(movement == Left) {
               // 向左移動,左側的非遮罩位置會變?yōu)樾录尤胛恢?                mask.foreachWestColumn { y =>
                  // 反算回實際行
                  val yTile = y + (_row - extent)
                  if(0 <= yTile && yTile < rows) {
                    f(addedCol, yTile)
                  }
                }
              } else { // 向右移動
                // ...省略相似邏輯
              }
            }
          }
        }
        
        // 處理游標內部的其他變化區(qū)域
        if(hasMask) {
          mask.foreachUnmasked(movement) { (x, y) =>
            val xTile = x + (_col - extent)
            val yTile = y + (_row - extent)
            if(0 <= xTile && xTile < cols && 0 <= yTile && yTile < rows) {
              f(xTile, yTile)
            }
          }
        }
    }
    
    // ... 省略其他

}

核心為對游標遮罩foreachWestColumn和foreachUnmasked方法的調用:

代碼位于[geotrellis/raster/mapalgebra/focal/CursorMask.scala]

class CursorMask(d:Int,f:(Int,Int)=>Boolean) {
    // ... 省略
    
    def foreachWestColumn(f:Int=>Unit) = {
        val len = westColumnUnmasked(0)
        if(len > 0) {
          var y = 1
          // 遍歷游標最左列中每一個非遮罩位置
          while(y <= len) {
            f(westColumnUnmasked(y))
            y += 1
          }
        }
    }
    
    private def foreach(xOffset:Int,yOffset:Int,startY:Int,set:MaskSet)(f:(Int,Int)=>Unit) = {
        var y = startY
        while(y < d) {
          set.foreachAt(y) { x => f(x - xOffset, y - yOffset) }
          y += 1
        }
    }
    
    def foreachUnmasked(mv:Movement)(f:(Int,Int)=>Unit): Unit = {
        mv match {
          case Left => foreach(0,0,0,unmaskedAfterMoveLeft)(f)
          case Right => foreach(1,0,0,maskedAfterMoveLeft)(f)
          case Up => foreach(0,0,1,unmaskedAfterMoveUp)(f)
          case Down => foreach(0,1,1,maskedAfterMoveUp)(f)
          case _ =>
        }
    }
    
    // ... 省略
}

變化區(qū)域遍歷的流程是這樣的:

1. 通過遍歷"邊緣快查表",得到最左邊/最右邊新增的像元位置
2. 通過遍歷"變化快查表",得到其他位置新增的像元位置

為什么不能直接通過遍歷"變化快查表"得到包括左右邊緣的變化信息呢?因為當初定義"變化快查表"的時候,沒法處理邊緣區(qū)域這些移動后會越界的位置,因此無法記錄邊緣區(qū)域的變化.

在上面定義"快查表"的時候,我們可能有些疑惑:為什么只需要定義諸如unmaskedAfterMoveLeft之類向左/向上的方法,卻不需要定義對應的unmaskedAfterMoveRight等向右/向下的方法呢?我們從foreachUnmasked方法中就可以看到答案:unmaskedAfterMoveRight與偏移一列后的maskedAfterMoveLeft表達了同樣的意思.

至此,我們就大致了解了游標對象提供了哪些功能,以及實現原理,基于這些原理與機制,我們需要看看如何在這些基礎之上,構建各種實際的算子.

4.3 calcFunc方法與Max算子

我們再回到抽象Focal類的execute方法中.在定義完游標對象后,又定義了calcFunc方法.

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

sealed trait TargetCell extends Serializable

object TargetCell {
  object NoData extends TargetCell
  object Data extends TargetCell
  object All extends TargetCell
}

// ... 省略 

val calcFunc: () => Unit =
  target match {
    case TargetCell.All =>
      { () => calc(tile, cursor) }
    case TargetCell.Data =>
      { () =>
        calc(tile, cursor)
        // 把原始的Nodata值復制回來,避免算子使Nodata值發(fā)生改變
        if(!isData(r.get(cursor.focusCol, cursor.focusRow))){
          copyOriginalValue(cursor.focusCol, cursor.focusRow, cursor.col, cursor.row)
        }
      }
    case TargetCell.NoData =>
      // ... 省略與TargetCell.Data類似的邏輯
}

// ... 省略

其中的核心是對calc方法的調用,它們對應著具體算子的實現.我們以Max算子為例:

代碼位于[geotrellis/raster/mapalgebra/focal/Max.scala]

def calc(r: Tile, cursor: Cursor) = {
  var m = Int.MinValue
  cursor.allCells.foreach { (x, y) =>
    val v = r.get(x, y)
    if(v > m) { m = v }
  }
  // 注意這里的col/row (col = focusCol - analysisOffsetCols)
  // 即結果的寫入永遠從0,0開始,沒有偏移量
  resultTile.set(cursor.col, cursor.row, m)
}

Max算子的實現很簡單:遍歷當前游標中的每一個像元,得到最大值,然后將最大值賦值給游標當前的聚焦位置.

calcFunc還受目標策略(TargetCell)影響,如果用戶指定當前算子只針對非Nodata值或Nodata值,在計算完后還需要用舊值將原先的非Nodata值或Nodata值的位置恢復,避免算子使Nodata值發(fā)生改變,造成與預期不一樣的結果.這部分的邏輯主要調用了copyOriginalValue方法,我們再來分析方法的實現:

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

val copyOriginalValue: (Int, Int, Int, Int) => Unit =
if(!r.cellType.isFloatingPoint) {
    { (focusCol: Int, focusRow: Int, col: Int, row: Int) =>
        resultTile.set(col, row, r.get(focusCol, focusRow))
    }
} else {
    // ... 省略類似的邏輯
}

傳入的值分兩種:focusCol/focusRow,col/row (col = focusCol - analysisOffsetCols),我們將其還原到calcFunc的場景(假定為非浮點類數值):

case TargetCell.Data =>
  { () =>
    calc(tile, cursor)
    if(!isData(r.get(cursor.focusCol, cursor.focusRow))){
      resultTile.set(cursor.col, cursor.row, r.get(cursor.focusCol, cursor.focusRow))
    }
  }

恢復Nodata值的邏輯是這樣:

1. 檢測當前游標聚焦點位置是否為Nodata值
2. 如果是Nodata值,則將其賦予給當前游標聚焦點減去分析區(qū)偏移的位置.這也與結果的寫入邏輯相一致.

4.4 游標策略的定義

如果說游標對象和calc方法定義了游標在滑動到Tile中某一個位置時的靜態(tài)行為,游標策略(CursorStrategy)則描述了游標如何在整個Tile上滑動的動態(tài)過程.因此在抽象游標算子的最后一步,會將一系列參數傳入CursorStrategy.execute方法中執(zhí)行.我們來分析CursorStrategy.execute方法的定義

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

// 定義了3種游標滑動策略
sealed trait TraversalStrategy
// 之字形
case object ZigZagTraversalStrategy extends TraversalStrategy
// 線型
case object ScanLineTraversalStrategy extends TraversalStrategy
// 螺旋型
case object SpiralZagTraversalStrategy extends TraversalStrategy
// 默認策略為之字形
object TraversalStrategy {
  def DEFAULT: TraversalStrategy = ZigZagTraversalStrategy
}

object CursorStrategy {
  
  // ... 省略
  
  def execute(
    cursor: Cursor,
    calc: () => Unit,
    analysisArea: GridBounds[Int],
    traversalStrategy: TraversalStrategy
  ): Unit = {
    traversalStrategy match {
      case ZigZagTraversalStrategy => handleZigZag(analysisArea, cursor, calc)
      case ScanLineTraversalStrategy => handleScanLine(analysisArea, cursor, calc)
      case SpiralZagTraversalStrategy => handleSpiralZag(analysisArea, cursor, calc)
    }
  }
  
  // ...省略其他
}

可見,execute方法的核心就是使用不同的策略操游標進行滑動的過程.我們以默認的之字形策略為例:

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

  private def handleZigZag(analysisArea: GridBounds[Int], cursor: Cursor, calc: () => Unit) = {
    val colMax = analysisArea.colMax
    val rowMax = analysisArea.rowMax
    val colMin = analysisArea.colMin
    val rowMin = analysisArea.rowMin

    var col = colMin
    var row = rowMin

    var direction = 1
    // 初始化游標的位置
    cursor.centerOn(col, row)

    while(row <= rowMax) {
      // 每滑動一次都進行一次計算,并記錄結果
      calc()
      // 先向右滑動
      col += direction
      // 移動越界了,就開始換行并轉向
      if(col < colMin || colMax < col) {
        // 滑動到正下方位置
        direction *= -1
        row += 1
        col += direction
        cursor.move(Movement.Down)
      } else {
        if(direction == 1) { cursor.move(Movement.Right) }
        else { cursor.move(Movement.Left) }
      }
    }
  }

不同策略的滑動軌跡如下:

之字形策略              線型策略                螺旋型策略
→   →   →   ↓           →   →   →   →           →   →   →   ↓
↓   ←   ←   ←           →   →   →   →           →   →   →   ↓
→   →   →   →           →   →   →   →           ↑   ←   ←   ←

5. 總結

我們再整體梳理一下Max算子:

代碼位于[geotrellis/raster/mapalgebra/focal/FocalMethods.scala]

object Max {
  def calculation(tile: Tile, n: Neighborhood, bounds: Option[GridBounds[Int]] = None, target: TargetCell = TargetCell.All): FocalCalculation[Tile] = {
    if (tile.cellType.isFloatingPoint) {
      new CursorCalculation[Tile](tile, n, bounds, target)
        with ArrayTileResult
      {
        def calc(r: Tile, cursor: Cursor) = {
            // ... 省略
        }
      }

    } else {
      // ... 省略處理Int類型數據時類似的邏輯
    }
  }

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

1. Max算子接收原始Tile和鄰域信息作為參數,輸出結果Tile.本質上是對一個Focal類算子執(zhí)行execute方法.
2. 根據最大值算法的特點,Max算子實現為Focal類算子的一個子類:游標類算子,并且根據其返回值的特性,混入了ArrayTileResult特質.
3. 為成為一個實際的游標類算子,Max實現了預定的calc方法,Max算法的核心就在這里.
4. 用戶實際調用Max算子時,游標類算子會創(chuàng)建一個滑動游標,在用戶傳入的Tile上按照指定的策略滑動,每次滑動后都會計算當前滑動窗口內通過定義的calc方法得到的結果值,并存入結果Tile中.
5. 當滑動完成,用戶會得到最終的結果Tile.

Max算子通過類似Local類算子的方法掛載:

代碼位于[geotrellis/raster/mapalgebra/focal/FocalMethods.scala]

trait FocalMethods extends MethodExtensions[Tile] {
  // ... 省略
  def focalMax(n: Neighborhood, bounds: Option[GridBounds[Int]] = None, target: TargetCell = TargetCell.All): Tile = {
    Max(self, n, bounds, target)
  }
  // ... 省略
}

至此, 我們也大致了解了游標類算子的結構與運行流程.