getEuclideanDistance.h

Go to the documentation of this file.

#ifndef _DLRCOMPUTERVISION_EUCLIDEANDISTANCE_H_
#define _DLRCOMPUTERVISION_EUCLIDEANDISTANCE_H_

#include <list>
#include <dlrComputerVision/imageFormat.h>
#include <dlrComputerVision/imageIO.h>
#include <dlrComputerVision/image.h>

namespace dlr {

  namespace computerVision {
    
    template<ImageFormat FORMAT>
    Array2D<double>
    getEuclideanDistance(const Image<FORMAT>& inputImage,
                         size_t maxNumberOfPasses=10);


    template<ImageFormat FORMAT>
    Array2D<double>
    getEuclideanDistance(const Image<FORMAT>& inputImage,
                         size_t maxNumberOfPasses,
                         size_t& numberOfPassesUsed);
  
  } // namespace computerVision

} // namespace dlr


/* ============ Definitions of inline & template functions ============ */


#include <cmath>
#include <dlrNumeric/index2D.h>

namespace dlr {

  namespace computerVision {
    
    namespace privateCode {

      inline bool
      eucDistPropagate(size_t toIndex,
                       size_t fromIndex,
                       Array2D<double>& distanceMap,
                       Array2D<Index2D>& referentMap,
                       int row,
                       int column)
      {
        // Inexpensive check to see if we _might_ need to update.
        if(distanceMap(fromIndex) < distanceMap(toIndex)) {
          Index2D& referent = referentMap(fromIndex);
          double deltaU = column - referent.getColumn();
          double deltaV = row - referent.getRow();
          double newDistance = deltaU * deltaU + deltaV * deltaV;
          if(newDistance < distanceMap(toIndex)) {
            distanceMap(toIndex) = newDistance;
            referentMap(toIndex) = referent;
            return true;
          }
        }
        return false;
      }

    } // namespace privateCode

  
    template<ImageFormat FORMAT>
    Array2D<double>
    getEuclideanDistance(const Image<FORMAT>& inputImage,
                         size_t maxNumberOfPasses)
    {
      size_t numberOfPassesUsed;
      return getEuclideanDistance(inputImage, maxNumberOfPasses,
                                  numberOfPassesUsed);
    }


    template<ImageFormat FORMAT>
    Array2D<double>
    getEuclideanDistance(const Image<FORMAT>& inputImage,
                         size_t maxNumberOfPasses,
                         size_t& numberOfPassesUsed)
    {
      Array2D<double> distanceMap(inputImage.rows(), inputImage.columns());
      Array2D<Index2D> referentMap(inputImage.rows(), inputImage.columns());

      // Initialize the distance matrix.
      double maxSqDistance = (inputImage.rows() * inputImage.rows()
                              + inputImage.columns() * inputImage.columns());
      size_t index0 = 0;
      for(size_t row = 0; row < inputImage.rows(); ++row) {
        for(size_t column = 0; column < inputImage.columns(); ++column) {
          if(inputImage(index0)) {
            distanceMap(index0) = 0.0;
            referentMap(index0) = Index2D(static_cast<int>(row), static_cast<int>(column));
          } else {
            distanceMap(index0) = maxSqDistance;
          }
          ++index0;
        }
      }

      // This is hard to do efficiently closed-form, so we'll iterate.
      // If numberOfPasses is set to 1, we'll still get a reasonable
      // map.
      size_t columns = distanceMap.columns();
      size_t rows = distanceMap.rows();
      size_t passNumber;
      for(passNumber = 0; passNumber < maxNumberOfPasses; ++passNumber) {
        // This variable will tell us if we can quit early because the
        // distances are all correct.
        bool isChanged = false;
      
        // === Propagate distances East. === 

        // Propagate distances East along top row.
        index0 = 1;
        for(size_t column = 1; column < inputImage.columns(); ++column) {
          isChanged |= privateCode::eucDistPropagate(
            index0, index0 - 1, distanceMap, referentMap,
            0, static_cast<int>(column));
          isChanged |= privateCode::eucDistPropagate(
            index0, index0 + columns - 1, distanceMap, referentMap,
            0, static_cast<int>(column));
          ++index0;
        }

        // Propagate distances East through the bulk of the image.
        for(size_t row = 1; row < inputImage.rows() - 1; ++row) {
          index0 = row * inputImage.columns() + 1;
          for(size_t column = 1; column < inputImage.columns(); ++column) {
            isChanged |= privateCode::eucDistPropagate(
              index0, index0 - columns - 1, distanceMap, referentMap,
              static_cast<int>(row), static_cast<int>(column));
            isChanged |= privateCode::eucDistPropagate(
              index0, index0 - 1, distanceMap, referentMap,
              static_cast<int>(row), static_cast<int>(column));
            isChanged |= privateCode::eucDistPropagate(
              index0, index0 + columns - 1, distanceMap, referentMap,
              static_cast<int>(row), static_cast<int>(column));
            ++index0;
          }
        }

        // Propagate distances East along bottom row.
        index0 = (inputImage.rows() - 1) * inputImage.columns() + 1;
        for(size_t column = 1; column < inputImage.columns(); ++column) {
          isChanged |= privateCode::eucDistPropagate(
            index0, index0 - columns - 1, distanceMap, referentMap,
            static_cast<int>(rows) - 1, static_cast<int>(column));
          isChanged |= privateCode::eucDistPropagate(
            index0, index0 - 1, distanceMap, referentMap,
            static_cast<int>(rows) - 1, static_cast<int>(column));
          ++index0;
        }

        // === Propagate distances South. ===
        index0 = columns + 1;
        for(size_t row = 1; row < inputImage.rows(); ++row) {
          index0 = row * inputImage.columns();

          // Update first pixel of the row.
          isChanged |= privateCode::eucDistPropagate(
            index0, index0 - columns, distanceMap, referentMap,
            static_cast<int>(row), 0);
          isChanged |= privateCode::eucDistPropagate(
            index0, index0 - columns + 1, distanceMap, referentMap,
            static_cast<int>(row), 0);
          ++index0;

          // Update the bulk of the pixels in the current row.
          for(size_t column = 1; column < inputImage.columns() - 1; ++column) {
            isChanged |= privateCode::eucDistPropagate(
              index0, index0 - columns - 1, distanceMap, referentMap,
              static_cast<int>(row), static_cast<int>(column));
            isChanged |= privateCode::eucDistPropagate(
              index0, index0 - columns, distanceMap, referentMap,
              static_cast<int>(row), static_cast<int>(column));
            isChanged |= privateCode::eucDistPropagate(
              index0, index0 - columns + 1, distanceMap, referentMap,
              static_cast<int>(row), static_cast<int>(column));
            ++index0;
          }

          // Update last pixel of the row.
          isChanged |= privateCode::eucDistPropagate(
            index0, index0 - columns - 1, distanceMap, referentMap,
            static_cast<int>(row), static_cast<int>(columns) - 1);
          isChanged |= privateCode::eucDistPropagate(
            index0, index0 - columns, distanceMap, referentMap,
            static_cast<int>(row), static_cast<int>(columns) - 1);
          ++index0;
        }

        // === Propagate distances West. ===

        // Propagate distances West along top row.
        index0 = columns - 2;

        // Trick here... note that column wraps around to a very large
        // number once instead of going negative.
        for(size_t column = inputImage.columns() - 2; column < columns;
            --column) {
          isChanged |= privateCode::eucDistPropagate(
            index0, index0 + 1, distanceMap, referentMap,
            0, static_cast<int>(column));
          isChanged |= privateCode::eucDistPropagate(
            index0, index0 + columns + 1, distanceMap, referentMap,
            0, static_cast<int>(column));
          --index0;
        }
      
        // Propagate distances West through the bulk of the image.
        for(size_t row = 1; row < inputImage.rows() - 1; ++row) {
          index0 = (row + 1) * inputImage.columns() - 2;
          for(size_t column = inputImage.columns() - 2; column < columns;
              --column) {
            isChanged |= privateCode::eucDistPropagate(
              index0, index0 - columns + 1, distanceMap, referentMap,
              static_cast<int>(row), static_cast<int>(column));
            isChanged |= privateCode::eucDistPropagate(
              index0, index0 + 1, distanceMap, referentMap,
              static_cast<int>(row), static_cast<int>(column));
            isChanged |= privateCode::eucDistPropagate(
              index0, index0 + columns + 1, distanceMap, referentMap,
              static_cast<int>(row), static_cast<int>(column));
            --index0;
          }
        }
      
        // Propagate distances West along bottom row.
        index0 = rows * columns - 2;
        for(size_t column = inputImage.columns() - 2; column < columns; --column) {
          isChanged |= privateCode::eucDistPropagate(
            index0, index0 - columns + 1, distanceMap, referentMap,
            static_cast<int>(rows - 1), static_cast<int>(column));
          isChanged |= privateCode::eucDistPropagate(
            index0, index0 + 1, distanceMap, referentMap,
            static_cast<int>(rows - 1), static_cast<int>(column));
          --index0;
        }


        // === Propagate distances North. ===

        // Remember that row wraps around to a very large number instead
        // of going negative.
        for(size_t row = inputImage.rows() - 2; row < rows; --row) {
          index0 = row * columns;

          // Update first pixel of the row.
          isChanged |= privateCode::eucDistPropagate(
            index0, index0 + columns, distanceMap, referentMap,
            static_cast<int>(row), 0);
          isChanged |= privateCode::eucDistPropagate(
            index0, index0 + columns + 1, distanceMap, referentMap,
            static_cast<int>(row), 0);
          ++index0;

          // Update the bulk of the pixels in the current row.
          for(size_t column = 1; column < inputImage.columns() - 1; ++column) {
            isChanged |= privateCode::eucDistPropagate(
              index0, index0 + columns - 1, distanceMap, referentMap,
              static_cast<int>(row), static_cast<int>(column));
            isChanged |= privateCode::eucDistPropagate(
              index0, index0 + columns, distanceMap, referentMap,
              static_cast<int>(row), static_cast<int>(column));
            isChanged |= privateCode::eucDistPropagate(
              index0, index0 + columns + 1, distanceMap, referentMap,
              static_cast<int>(row), static_cast<int>(column));
            ++index0;
          }

          // Update last pixel of the row.
          isChanged |= privateCode::eucDistPropagate(
            index0, index0 + columns - 1, distanceMap, referentMap,
            static_cast<int>(row), static_cast<int>(columns) - 1);
          isChanged |= privateCode::eucDistPropagate(
            index0, index0 + columns, distanceMap, referentMap,
            static_cast<int>(row), static_cast<int>(columns) - 1);
          ++index0;
        }

        // If no pixels needed changing, we're done.
        if(!isChanged) {
          break;
        }
      }

      // Convert squared distance to distance.
      for(index0 = 0; index0 < distanceMap.size(); ++index0) {
        distanceMap(index0) = std::sqrt(distanceMap(index0));
      }

      numberOfPassesUsed = passNumber;
      return distanceMap;
    }
  
  } // namespace computerVision
  
} // namespace dlr

#endif /* #ifndef _DLRCOMPUTERVISION_EUCLIDEANDISTANCE_H_ */