locate_cliff_toe.cpp

Go to the documentation of this file.

 
/* ==============================================================================================================================
 
   This file is part of CoastalME, the Coastal Modelling Environment.
 
   CoastalME is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation; either version 3 of the License, or (at your option) any later
version.
 
   This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
 
   You should have received a copy of the GNU General Public License along with
this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave,
Cambridge, MA 02139, USA.
 
==============================================================================================================================*/
#include <assert.h>
#include <cfloat>
 
#include <cpl_port.h>
#include <iostream>
using std::cerr;
using std::cout;
using std::endl;
 
#include <iomanip>
using std::resetiosflags;
using std::setiosflags;
using std::setprecision;
using std::setw;
 
#include <stdint.h>
 
#include <sstream>
using std::stringstream;
 
#include "cme.h"
#include "raster_grid.h"
#include "simulation.h"
 
/*===============================================================================================================================
 
 Calculates slope at every cell throughout the grid using finite difference
method
 
===============================================================================================================================*/
void CSimulation::nCalcSlopeAtAllCells(void) {
  for (int nX = 0; nX < m_nXGridSize; nX++) {
    for (int nY = 0; nY < m_nYGridSize; nY++) {
 
      // Now lets look at surounding cells
      if (nX > 0 && nX < m_nXGridSize - 1 && nY > 0 && nY < m_nYGridSize - 1) {
        double dElevLeft =
            m_pRasterGrid->m_Cell[nX - 1][nY].dGetSedimentTopElev();
        double dElevRight =
            m_pRasterGrid->m_Cell[nX + 1][nY].dGetSedimentTopElev();
        double dElevUp =
            m_pRasterGrid->m_Cell[nX][nY - 1].dGetSedimentTopElev();
        double dElevDown =
            m_pRasterGrid->m_Cell[nX][nY + 1].dGetSedimentTopElev();
        // Calculate slope using finite difference method
        double dSlopeX = (dElevRight - dElevLeft) / (2.0 * m_dCellSide);
        double dSlopeY = (dElevDown - dElevUp) / (2.0 * m_dCellSide);
        double dSlope = sqrt(dSlopeX * dSlopeX + dSlopeY * dSlopeY);
        m_pRasterGrid->m_Cell[nX][nY].SetSlope(dSlope);
      }
    }
  }
}
 
/*===============================================================================================================================
 
Highligts cells with slope greater than threshold
 
===============================================================================================================================*/
void CSimulation::nLocateCliffCell(void) {
  for (int nX = 0; nX < m_nXGridSize; nX++) {
    for (int nY = 0; nY < m_nYGridSize; nY++) {
      double dSlope = m_pRasterGrid->m_Cell[nX][nY].dGetSlope();
      if (dSlope >= m_dCliffSlopeLimit) {
        m_pRasterGrid->m_Cell[nX][nY].SetAsCliff(true);
      }
    }
  }
}
/*===============================================================================================================================
 
Removes small cliff islands below a given number of cellsusing flood fill
algorithm
 
===============================================================================================================================*/
void CSimulation::nRemoveSmallCliffIslands(int const dMinCliffCellThreshold) {
  // unsigned int m_nXGridSize = static_cast<unsigned int>(m_nXGridSize);
  // unsigned int m_nYGridSize = static_cast<unsigned int>(m_nYGridSize);
 
  // Create a 2D array to track which cells have been visited during flood fill
  vector<vector<bool>> bVisited(static_cast<unsigned int>(m_nXGridSize),
                                vector<bool>(static_cast<unsigned int>(m_nYGridSize), false));
 
  // Vector to store cells that belong to small cliff islands to be removed
  vector<pair<int, int>> VSmallIslandCells;
 
  // Loop through all cells to find unvisited cliff cells
  for (unsigned int nX = 0; nX < static_cast<unsigned int>(m_nXGridSize); nX++) {
    for (unsigned int nY = 0; nY < static_cast<unsigned int>(m_nYGridSize); nY++) {
      // Check if this is an unvisited cliff cell
      if (!bVisited[nX][nY] && m_pRasterGrid->m_Cell[nX][nY].bIsCliff()) {
        // Found the start of a new cliff region - use flood fill to find all
        // connected cliff cells
        vector<pair<int, int>> VCurrentCliffRegion;
 
        // Stack for iterative flood fill algorithm
        vector<pair<int, int>> VStack;
        VStack.push_back(std::make_pair(nX, nY));
 
        // Flood fill to find all connected cliff cells
        while (!VStack.empty()) {
          pair<int, int> currentCell = VStack.back();
          VStack.pop_back();
 
          size_t nCurX = static_cast<size_t>(currentCell.first);
          size_t nCurY = static_cast<size_t>(currentCell.second);
 
          // Skip if already visited or out of bounds
          if (nCurX >= static_cast<unsigned int>(m_nXGridSize) ||
              nCurY >= static_cast<unsigned int>(m_nYGridSize) || bVisited[nCurX][nCurY] ||
              !m_pRasterGrid->m_Cell[nCurX][nCurY].bIsCliff()) {
            continue;
          }
 
          // Mark as visited and add to current cliff region
          bVisited[nCurX][nCurY] = true;
          VCurrentCliffRegion.push_back(std::make_pair(nCurX, nCurY));
 
          // Add neighboring cells to stack (8-connectivity: N, NE, E, SE, S,
          // SW, W, NW)
          VStack.push_back(std::make_pair(nCurX - 1, nCurY));     // North
          VStack.push_back(std::make_pair(nCurX - 1, nCurY + 1)); // Northeast
          VStack.push_back(std::make_pair(nCurX, nCurY + 1));     // East
          VStack.push_back(std::make_pair(nCurX + 1, nCurY + 1)); // Southeast
          VStack.push_back(std::make_pair(nCurX + 1, nCurY));     // South
          VStack.push_back(std::make_pair(nCurX + 1, nCurY - 1)); // Southwest
          VStack.push_back(std::make_pair(nCurX, nCurY - 1));     // West
          VStack.push_back(std::make_pair(nCurX - 1, nCurY - 1)); // Northwest
        }
 
        // Calculate area of this cliff region (number of cells * cell area)
        int dCliffRegionArea = static_cast<int>(VCurrentCliffRegion.size());
 
        // If area is below threshold, mark all cells in this region for removal
        if (dCliffRegionArea < dMinCliffCellThreshold) {
          VSmallIslandCells.insert(VSmallIslandCells.end(),
                                   VCurrentCliffRegion.begin(),
                                   VCurrentCliffRegion.end());
        }
      }
    }
  }
 
  // Remove cliff designation from all small island cells
  for (const auto &cell : VSmallIslandCells) {
    m_pRasterGrid->m_Cell[cell.first][cell.second].SetAsCliff(false);
  }
}
 
/*===============================================================================================================================
 
Traces the seaward extent of cliff cells using wall following algorithm
 
===============================================================================================================================*/
void CSimulation::nTraceSeawardCliffEdge(void) {
  // Clear previous cliff edges
  m_VCliffEdge.clear();
 
  // Find all possible cliff edge start points by scanning for cliff/non-cliff
  // transitions
  vector<CGeom2DIPoint> V2DIPossibleStartCell;
  vector<bool> VbPossibleStartCellHandedness;
  vector<int> VnSearchDirection;
 
  // Scan all cells to find cliff toe edges by checking elevation patterns
  for (int nX = 2; nX < m_nXGridSize - 2; nX++) {
    for (int nY = 2; nY < m_nYGridSize - 2; nY++) {
      if (m_pRasterGrid->m_Cell[nX][nY].bIsCliff()) {
 
        // East direction (check if this is a seaward-facing cliff toe)
        if (!m_pRasterGrid->m_Cell[nX][nY + 1].bIsCliff()) {
          V2DIPossibleStartCell.push_back(CGeom2DIPoint(nX, nY));
          VbPossibleStartCellHandedness.push_back(true);
          VnSearchDirection.push_back(EAST);
        }
 
        // South direction
        if (!m_pRasterGrid->m_Cell[nX + 1][nY].bIsCliff()) {
          V2DIPossibleStartCell.push_back(CGeom2DIPoint(nX, nY));
          VbPossibleStartCellHandedness.push_back(true);
          VnSearchDirection.push_back(SOUTH);
        }
 
        // West direction
        if (!m_pRasterGrid->m_Cell[nX][nY - 1].bIsCliff()) {
          V2DIPossibleStartCell.push_back(CGeom2DIPoint(nX, nY));
          VbPossibleStartCellHandedness.push_back(true);
          VnSearchDirection.push_back(WEST);
        }
 
        // North direction
        if (!m_pRasterGrid->m_Cell[nX - 1][nY].bIsCliff()) {
          V2DIPossibleStartCell.push_back(CGeom2DIPoint(nX, nY));
          VbPossibleStartCellHandedness.push_back(true);
          VnSearchDirection.push_back(NORTH);
        }
      }
    }
  }
 
  // Create a 2D array to track which cells have been used in cliff edge tracing
  vector<vector<bool>> bUsedInCliffTrace(m_nXGridSize,
                                         vector<bool>(m_nYGridSize, false));
 
  int nCliffEdgesTraced = 0;
 
  // For each possible start point, attempt to trace a cliff edge
  for (size_t nStartPoint = 0; nStartPoint < V2DIPossibleStartCell.size();
       nStartPoint++) {
    int nXStart = V2DIPossibleStartCell[nStartPoint].nGetX();
    int nYStart = V2DIPossibleStartCell[nStartPoint].nGetY();
 
    // Skip if this cell has already been used in another cliff trace
    if (bUsedInCliffTrace[nXStart][nYStart]) {
      continue;
    }
 
    // Begin cliff edge tracing using wall following algorithm
    vector<CGeom2DIPoint> VCliffEdge;
    int nSearchDirection = VnSearchDirection[nStartPoint];
 
    int nX = nXStart;
    int nY = nYStart;
    int nLength = 0;
    const int nMaxLength = m_nXGridSize * m_nYGridSize; // Safety limit
 
    do {
      // Add current point to cliff edge
      VCliffEdge.push_back(CGeom2DIPoint(nX, nY));
      bUsedInCliffTrace[nX][nY] = true;
      nLength++;
 
      // Find next cell using wall following algorithm (right-hand rule)
      int nXSeaward, nYSeaward, nXStraightOn, nYStraightOn;
      int nXAntiSeaward, nYAntiSeaward, nXGoBack, nYGoBack;
 
      // Calculate candidate positions based on search direction
      switch (nSearchDirection) {
      case NORTH:
        nXSeaward = nX + 1;
        nYSeaward = nY; // Right (East)
        nXStraightOn = nX;
        nYStraightOn = nY - 1; // Straight (North)
        nXAntiSeaward = nX - 1;
        nYAntiSeaward = nY; // Left (West)
        nXGoBack = nX;
        nYGoBack = nY + 1; // Back (South)
        break;
      case EAST:
        nXSeaward = nX;
        nYSeaward = nY + 1; // Right (South)
        nXStraightOn = nX + 1;
        nYStraightOn = nY; // Straight (East)
        nXAntiSeaward = nX;
        nYAntiSeaward = nY - 1; // Left (North)
        nXGoBack = nX - 1;
        nYGoBack = nY; // Back (West)
        break;
      case SOUTH:
        nXSeaward = nX - 1;
        nYSeaward = nY; // Right (West)
        nXStraightOn = nX;
        nYStraightOn = nY + 1; // Straight (South)
        nXAntiSeaward = nX + 1;
        nYAntiSeaward = nY; // Left (East)
        nXGoBack = nX;
        nYGoBack = nY - 1; // Back (North)
        break;
      case WEST:
        nXSeaward = nX;
        nYSeaward = nY - 1; // Right (North)
        nXStraightOn = nX - 1;
        nYStraightOn = nY; // Straight (West)
        nXAntiSeaward = nX;
        nYAntiSeaward = nY + 1; // Left (South)
        nXGoBack = nX + 1;
        nYGoBack = nY; // Back (East)
        break;
      default:
        nXSeaward = nXStraightOn = nXAntiSeaward = nXGoBack = nX;
        nYSeaward = nYStraightOn = nYAntiSeaward = nYGoBack = nY;
        break;
      }
 
      // Try to move using wall following priority: seaward, straight,
      // anti-seaward, back
      bool bFoundNextCell = false;
 
      // 1. Try seaward (right turn)
      if (bIsWithinValidGrid(nXSeaward, nYSeaward) &&
          m_pRasterGrid->m_Cell[nXSeaward][nYSeaward].bIsCliff()) {
        nX = nXSeaward;
        nY = nYSeaward;
        // Update search direction (turn right)
        switch (nSearchDirection) {
        case NORTH:
          nSearchDirection = EAST;
          break;
        case EAST:
          nSearchDirection = SOUTH;
          break;
        case SOUTH:
          nSearchDirection = WEST;
          break;
        case WEST:
          nSearchDirection = NORTH;
          break;
        }
        bFoundNextCell = true;
      }
      // 2. Try straight ahead
      else if (bIsWithinValidGrid(nXStraightOn, nYStraightOn) &&
               m_pRasterGrid->m_Cell[nXStraightOn][nYStraightOn].bIsCliff()) {
        nX = nXStraightOn;
        nY = nYStraightOn;
        // Direction stays the same
        bFoundNextCell = true;
      }
      // 3. Try anti-seaward (left turn)
      else if (bIsWithinValidGrid(nXAntiSeaward, nYAntiSeaward) &&
               m_pRasterGrid->m_Cell[nXAntiSeaward][nYAntiSeaward].bIsCliff()) {
        nX = nXAntiSeaward;
        nY = nYAntiSeaward;
        // Update search direction (turn left)
        switch (nSearchDirection) {
        case NORTH:
          nSearchDirection = WEST;
          break;
        case EAST:
          nSearchDirection = NORTH;
          break;
        case SOUTH:
          nSearchDirection = EAST;
          break;
        case WEST:
          nSearchDirection = SOUTH;
          break;
        }
        bFoundNextCell = true;
      }
      // 4. Try going back (U-turn)
      else if (bIsWithinValidGrid(nXGoBack, nYGoBack) &&
               m_pRasterGrid->m_Cell[nXGoBack][nYGoBack].bIsCliff()) {
        nX = nXGoBack;
        nY = nYGoBack;
        // Update search direction (U-turn)
        nSearchDirection = nGetOppositeDirection(nSearchDirection);
        bFoundNextCell = true;
      }
 
      if (!bFoundNextCell) {
        break; // No valid next cell found, end this cliff edge trace
      }
 
    } while ((nX != nXStart || nY != nYStart) && nLength < nMaxLength);
 
    // Only keep cliff edges that have a reasonable length
    if (VCliffEdge.size() > 2) {
      nCliffEdgesTraced++;
 
      // Convert grid coordinates to external CRS for smoothing
      CGeomLine CliffEdgeExtCRS;
      for (const auto &point : VCliffEdge) {
        CliffEdgeExtCRS.Append(dGridCentroidXToExtCRSX(point.nGetX()),
                               dGridCentroidYToExtCRSY(point.nGetY()));
        bUsedInCliffTrace[point.nGetX()][point.nGetY()] = true;
      }
 
      // Apply cliff edge specific smoothing
      if (m_nCliffEdgeSmooth == SMOOTH_RUNNING_MEAN)
        CliffEdgeExtCRS = LSmoothCoastRunningMean(&CliffEdgeExtCRS);
      else if (m_nCliffEdgeSmooth == SMOOTH_SAVITZKY_GOLAY)
        CliffEdgeExtCRS = LSmoothCoastSavitzkyGolay(&CliffEdgeExtCRS, 0, 0);
 
      // Store the smoothed cliff edge in external CRS
      m_VCliffEdge.push_back(CliffEdgeExtCRS);
    }
  }
}
 
/*===============================================================================================================================
 
Validates cliff edges to ensure they represent cliff toes, truncating edges that
trace cliff tops
 
===============================================================================================================================*/
void CSimulation::nValidateCliffToeEdges(void) {
  vector<CGeomLine> ValidatedCliffEdges;
 
  // Process each traced cliff edge
  for (size_t nEdge = 0; nEdge < m_VCliffEdge.size(); nEdge++) {
    CGeomLine &CliffEdge = m_VCliffEdge[nEdge];
 
    // Try validating in forward direction first
    CGeomLine ForwardValidated = nValidateCliffToeDirection(CliffEdge, false);
 
    // If we got a very short result (broke early), try reverse direction
    CGeomLine ReverseValidated;
    if (ForwardValidated.nGetSize() < CliffEdge.nGetSize() * 0.2) {
      ReverseValidated = nValidateCliffToeDirection(CliffEdge, true);
    }
 
    // Use whichever result is longer
    CGeomLine BestValidated;
    if (ReverseValidated.nGetSize() > ForwardValidated.nGetSize()) {
      BestValidated = ReverseValidated;
    } else {
      BestValidated = ForwardValidated;
    }
 
    // Only keep edges that have a reasonable length after validation
    if (BestValidated.nGetSize() > 2) {
      ValidatedCliffEdges.push_back(BestValidated);
    }
  }
 
  // Replace the original cliff edges with the validated ones
  m_VCliffEdge = ValidatedCliffEdges;
}
 
CGeomLine CSimulation::nValidateCliffToeDirection(CGeomLine &CliffEdge, bool bReverse) {
  CGeomLine ValidatedEdge;
 
  int nConsecutiveFailures = 0;
  const int nMaxConsecutiveFailures = 2;
 
  int nSize = CliffEdge.nGetSize();
 
  // Check each point along the cliff edge
  for (int i = 0; i < nSize - 1; i++) {
    // Determine which point to process based on direction
    int nPoint = bReverse ? (nSize - 1 - i) : i;
    int nNextPoint = bReverse ? (nSize - 2 - i) : (i + 1);
 
    // Skip if we've gone beyond bounds
    if (nNextPoint < 0 || nNextPoint >= nSize)
      continue;
 
    // Get current and next point in grid coordinates
    int nX =
        static_cast<int>((CliffEdge.dGetXAt(nPoint) - m_dGeoTransform[0]) /
                         m_dGeoTransform[1]);
    int nY =
        static_cast<int>((CliffEdge.dGetYAt(nPoint) - m_dGeoTransform[3]) /
                         m_dGeoTransform[5]);
 
    int nNextX = static_cast<int>(
        (CliffEdge.dGetXAt(nNextPoint) - m_dGeoTransform[0]) /
        m_dGeoTransform[1]);
    int nNextY = static_cast<int>(
        (CliffEdge.dGetYAt(nNextPoint) - m_dGeoTransform[3]) /
        m_dGeoTransform[5]);
 
    // Ensure coordinates are within grid bounds
    nX = std::max(0, std::min(m_nXGridSize - 1, nX));
    nY = std::max(0, std::min(m_nYGridSize - 1, nY));
    nNextX = std::max(0, std::min(m_nXGridSize - 1, nNextX));
    nNextY = std::max(0, std::min(m_nYGridSize - 1, nNextY));
 
    // Calculate direction of travel
    int nDirX = nNextX - nX;
    int nDirY = nNextY - nY;
 
    // Get perpendicular directions (90 degrees to travel direction)
    int nLeftX = nX - nDirY; // Rotate direction 90 degrees counter-clockwise
    int nLeftY = nY + nDirX;
    int nRightX = nX + nDirY; // Rotate direction 90 degrees clockwise
    int nRightY = nY - nDirX;
 
    bool bIsValidToe = true;
 
    // Check if perpendicular cells are within bounds
    if (bIsWithinValidGrid(nLeftX, nLeftY) &&
        bIsWithinValidGrid(nRightX, nRightY)) {
      bool bLeftIsCliff = m_pRasterGrid->m_Cell[nLeftX][nLeftY].bIsCliff();
      bool bRightIsCliff = m_pRasterGrid->m_Cell[nRightX][nRightY].bIsCliff();
 
      // One should be cliff and one should be not cliff for a valid cliff
      // edge
      if (bLeftIsCliff != bRightIsCliff) {
        // Get the elevation of these two adjacent cells
        double dLeftElev =
            m_pRasterGrid->m_Cell[nLeftX][nLeftY].dGetSedimentTopElev();
        double dRightElev =
            m_pRasterGrid->m_Cell[nRightX][nRightY].dGetSedimentTopElev();
 
        // Determine which is the cliff side and which is the non-cliff side
        double dCliffElev, dNonCliffElev;
        if (bLeftIsCliff) {
          dCliffElev = dLeftElev;
          dNonCliffElev = dRightElev;
        } else {
          dCliffElev = dRightElev;
          dNonCliffElev = dLeftElev;
        }
 
        // If the non-cliff cell is higher than the cliff cell, this is not
        // the toe
        if (dNonCliffElev > dCliffElev) {
          bIsValidToe = false;
        }
      }
    }
 
    if (bIsValidToe) {
      // Reset consecutive failure counter
      nConsecutiveFailures = 0;
      // Add this point to the validated edge
      ValidatedEdge.Append(CliffEdge.dGetXAt(nPoint),
                           CliffEdge.dGetYAt(nPoint));
    } else {
      // Increment consecutive failure counter
      nConsecutiveFailures++;
 
      // If we haven't hit the threshold yet, still add the point
      if (nConsecutiveFailures < nMaxConsecutiveFailures) {
        ValidatedEdge.Append(CliffEdge.dGetXAt(nPoint),
                             CliffEdge.dGetYAt(nPoint));
      } else {
        // Too many consecutive failures - truncate here
        break;
      }
    }
  }
 
  return ValidatedEdge;
}
 
/*===============================================================================================================================
 
 Locates and traces the cliff toe
 
===============================================================================================================================*/
int CSimulation::nLocateCliffToe(void) {
  // First step: calculate slope at every cell throughout the grid
  nCalcSlopeAtAllCells();
  nLocateCliffCell();
  nRemoveSmallCliffIslands(50);
  nTraceSeawardCliffEdge();
  nValidateCliffToeEdges();
 
  return RTN_OK;
}