Shortest path in image

Question

This code takes an image and detects a global shortest path from the top to bottom row, with the requirement that top and bottom column index be the same. For this, it scans through each element on the top row using Dijkstra's algorithm towards its corresponding point on the bottom row. The shortest path is defined as the sum of all passed pixel values. The constraint here is that each pixel only has three neighbors in the next lower row; lower left, straight down, and lower right. This makes any expansion triangular. Each of those three neighborhood pixels gets the column index of the current pixel as a predecessor (we automatically know the predecessor row is one higher).

I ran a brute-force approach, calculating the shortest path by running along all rows and columns sequentially (pseudocode):

for every element (0,r) in the first row
   for rows i (1 : end)
       for cols j (r-i : r+i) //triangular expansion
          get distances of three elements in row above, so d (dist = infinite if not part of the triangular expansion)
          d[i,j] <- min(d[i-1,j-1],d[i-1,j],d[i-1,j+1]) + image[i,j]
          predecessor[i,j] <- column index of smallest of the three
   once final row is expanded get d[end,r]
   update global distance if d < previous global distance and corresponding global r
backtrack along predecessor matrix of smallest d starting from element (end, global r)

which takes about 3 seconds per image, while using Dijkstra's algorithm takes about 12 seconds:

#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/opencv.hpp" 
#include <stdlib.h>
#include <stdio.h>
#include <math.h> ///round
#include <iostream> ///cout
#include <time.h>
#include <fstream>
#include <queue>

using namespace std;
using namespace cv;

const int inf = 0x7F800000;

struct Node
{
    Point index;              // index of node in graph
    double distance;  // distance from source (only allow positive distances)
};

struct CompareDist
{
    bool operator()(Node const& n1, Node const& n2)
    {
        return (n1.distance > n2.distance);
    }
};

/// Function headers
double dijkstra(Point, Point, const cv::Mat&, cv::Mat&);

int main( int argc, char** argv )
{

//initialize
Mat image= imread( "C:\\pics\\test.tif" ); //LV input image

//shorten calculations by setting an upper bound for shortest paths. 
//the global shortest path can not be larger than the sum of the pixels running in straight line top-down.
Mat columnwiseSum(cv::Mat::zeros(1,image.cols,CV_32F));;
reduce(image, columnwiseSum, 0, CV_REDUCE_SUM, -1); //gives a row of summed columns
double min, max;
Point min_loc, max_loc;
double upperLimit;
minMaxLoc(columnwiseSum, &min, &max, &min_loc, &max_loc);
upperLimit = columnwiseSum.at<float>(min_loc); // this is the first upper bound

//initialize
Mat predecessorMatrixTemp = Mat::zeros(image.size(),CV_32F); //gets updated each run of dijkstra
Mat predecessorMatrix = Mat::zeros(image.size(),CV_32F); //gets updated only if a new global shortest path is found
double dist = inf; //global shortest path distance
double distTemp = inf; //shortest path distance for each dikstra run
int index = 0; //  corresponds to global shortest path column 


//run through every point   
for (int i = 0; i<image.cols; i++){
    Point source = Point(i,0);  // set start and end points in current column i
    Point goal = Point(i,image.rows-1);
    distTemp= dijkstra(source, goal, image,predecessorMatrixTemp, upperLimit); //shortest path for this column


    //if new path is smaller than previous shortest global path, update .
    if (distTemp < dist){
        index = i;
        dist = distTemp;
        upperLimit = distTemp; //set new upper bound 
        predecessorMatrixTemp.copyTo(predecessorMatrix);
    }
}


cout << "Final distance: " << dist << " at column " << index << endl;

//here's code to backtrack along the shortest global path, this works


return 0;
}



// priority queue dijkstra
double dijkstra(Point source, Point goal, const cv::Mat& image, cv::Mat& predecessor, double upperLimit)
{   
//reset predecessor matrix to zero
predecessor = Mat::zeros(image.size(),CV_32F);  

// initialize the distance of each node to infinity
Mat distance = Mat::ones(image.size(),CV_32F);
multiply(distance,inf,distance);
cv::Rect rect(cv::Point(), image.size()); //bounding rectangle ensure s that indexes don't get out of image

// the distance of the source is its value
distance.at<float>(source) = image.at<float>(source); 

// create priority queue structure
std::priority_queue< Node, std::vector< Node >, CompareDist> pQueue;

// enqueue source node with beginning distance
Node first = { source, image.at<float>(source) };
pQueue.push(first);

// take lowest distance priority queue node
while(!pQueue.empty()) 
{
    Node tempNode = pQueue.top();
    pQueue.pop(); //remove this element
    Point nodeIndex = tempNode.index; // get element index 
    if (nodeIndex == goal){ // found the path to goal
        return tempNode.distance;
        }

    int newX, newY; //indices for neighborhood node
    for(int i = -1; i < 2; i++) //for every neighborhood element 
        {
        newY = nodeIndex.y+1; //new row
        newX = nodeIndex.x+i; //new col

        // Update the distance if it is smaller than the current distance

        // this constrains the expansion within a diamond shape, since only certain neighboring nodes are allowed 
        if (newY <= floor(image.rows/2) && (newX >= source.x-newY && newX <= source.x+newY) || //upper half 
        (newY > floor(image.rows/2) && (newX >=source.x-(goal.y-newY) && newX <=source.x+(goal.y-newY))))   // lower half 
            {
                // if new node is inside image boundaries
             if (rect.contains(Point(newX,newY)) ) {   
                double tempDist = tempNode.distance + image.at<float>(Point(newX,newY)); //temporary distance of neighbor = current distance + pixel value 
                //if the new distance is smaller than previous distances to this pixel, and smaller than the upper bound, add new node to queue
                if(distance.at<float>(Point(newX,newY)) > tempDist && tempDist <= upperLimit) 
                {  
                    distance.at<float>(Point(newX,newY)) = tempDist; 
                    Node newNode;
                    newNode.index = Point(newX,newY);
                    newNode.distance = tempDist  ; 
                    predecessor.at<float>(Point(newX,newY)) = nodeIndex.x; //new node gets pointed to predecessor node column
                    pQueue.push(newNode);

                }
            }
        }

    }
}
//if no useful shortest path found
return inf;
}

Answer 1

Style

Using using namespace std is bad. Same goes for using namespace cv obviously.

Your indentation is off, but I'm guessing that this is a copy paste error.

Comments

I find your comments to not be helpful. They should explain why something is done, not what is done. If you need comments to explain what you're doing, then you need to write clearer code with better variable and function names. Possibly breaking out more functions.

Combine definition and initialization

Code like this:

double upperLimit;
upperLimit = columnwiseSum.at<float>(min_loc); // this is the first upper bound

has poor readability. You should structure it like this (again better name on variables):

double upperBound = columnwiseSum.at<float>(min_loc);

See how I combined the definition and initialization of upperBound? This makes your code easier to follow, easier to verify that the variable is always initialized and also shorter to write.

There are multiple occasions of this.

Node class

The name index is wildly misleading, rather it's the position in the image. I would suggest renaming it to position. If you need to comment a member variable, its name is not good enough. Consider renaming distance to distanceFromSource.

Instead of providing the explicit comparator CompareDist you could simply provide an overload for operator > and use std::greater<Node> in your priority queue.

In summary:

struct Node
{
    Point position;
    double distanceFromSource;

    bool operator > (const Node& n) const{
        return distance > n.distance;
    }
};

...

std::priority_queue< Node, std::vector< Node >, std::greater<Node>>

But that's not infinite...

This really isn't infinity...

const int inf = 0x7F800000;

you are also assigning this to double when you want infinity:

double dist = inf; //global shortest path distance
double distTemp = inf; //shortest path distance for each dikstra run

you should just use: std::numeric_limits<double>::infinity() and while you're at it, if you need to comment what a variable is, the name is bad. How about this:

double shortestPathDist = std::numeric_limits<double>::infinity();
double shortestPathDistCurrent = std::numeric_limits<double>::infinity();

Numerical inaccuracy

You are using float to accumulate the distances, after each floating point operation you introduce a truncation error on the order of x*E-6 where x is the stored value. If you'd simply use CV_32SC1 and image.at<int32_t>() you wouldn't have to worry about truncation errors affecting your graph search. You would probably also notice a speed up.

At any rate you're comparing floats to doubles and doing a lot of type-casting which is unnecessary and can take some time if done in a hot loop.

Mat initialisation

This:

// initialize the distance of each node to infinity
Mat distance = Mat::ones(image.size(),CV_32F);
multiply(distance,inf,distance);

going with the above advice to use numeric_limits should be:

Mat distance(image.size(), CV_32F);
distance = std::numeric_limits<float>::infinity();

Dijkstra implementation

Okay the name pQueue tells me something about the type but not what it is used for. A much better name would be unvisitedNodes.

The variable first is only used once to push it into the pQueue better to just do it like this:

unvisitedNodes.emplace(source, image.at<float>(source)); 

Always strive to reduce the number of variables, this makes it easier to reason about the state of the code (but not to absurdum ofcourse). Notice that this could be more readable if we renamed source to start which is more logical.

unvisitedNodes.emplace(start, image.at<float>(start));

Now this tells me that I construct a new unvisited node from the start point. Great! Do note that dijkstra's algorithm starts with a '0' distance. But what you're doing is equivalent to starting at a node with only one edge, to the source node. So this won't affect your results.

The name tempNode again is a terrible name. Better names are currentNode or shortestPathNode.

Please do fix up your indentation:

Point nodeIndex = tempNode.index; // get element index 
if (nodeIndex == goal){ // found the path to goal
    return tempNode.distance;
    }

int newX, newY; //indices for neighborhood node
for(int i = -1; i < 2; i++) //for every neighborhood element 
    {
    newY = nodeIndex.y+1; //new row
    newX = nodeIndex.x+i; //new col

This is a proper nightmare to read. Again if you need to comment the variables, choose better names. Trust me, you'll save more time when you debug the code that you spend writing slightly longer names.

Instead of int newY, newX you should just have Point newPosition. Would simplify lots when addressing the images.

This:

     // this constrains the expansion within a diamond shape, since only certain neighboring nodes are allowed 
    if (newY <= floor(image.rows/2) && (newX >= source.x-newY && newX <= source.x+newY) || //upper half 
    (newY > floor(image.rows/2) && (newX >=source.x-(goal.y-newY) && newX <=source.x+(goal.y-newY))))   // lower half 
        {

piece of readability massacre kept me from reviewing this code for a few days.

Lets, start with floor(image.rows/2). First of image.rows is an integer type so the division by 2 will always truncate towards zero, as the size is always positive this is exactly the same behavior as removing the call to floor.

What you are actually doing is:

1. dividing by two (automatic floor here)
2. convert to float
3. call floor
4. convert newY to float
5. compare newY to result of 3.

If you just remove the floor the result will be exactly the same but the steps will reduce to:

1. divide by two (automatic floor here)
2. compare newY to result of 1.

Now after much head scratching I have managed to decipher that you limit the search to a diamond (rhomb) shape with the apex and nadir at the start and end points. The rhomb has two lines making a hat shape at the top, which the point must be below to be inside and then two lines making v shape at the bottom which the point must be above to be inside.

The code can then be expressed in a much more readable form as:

const bool insideULEdge = (source.x - newX) >= -newY;
const bool insideUREdge = (source.x - newX) <=  newY;
const bool insideLLEdge = (goal.x - newX) >= image.rows - newY - 1;
const bool insideLREdge = (goal.x - newX) <= image.rows - newY - 1;
const bool insideDiamond = insideULEdge && insideUREdge && insideLLEdge && insideLREdge;
if(insideDiamond && rect.contains(Point(newX, newY)){

(note that goal.x == source.x so the above would simplify further.)

This should be extracted into a function called bool isInsideSearchSpace(const Point& start, const Point& end, const Mat& image).

Performance

You have implemented the Dijkstra's method wrong. You are visiting nodes more frequently than you have to. Please read the algorithm here: Dijkstra's algorithm.

You need to keep track of visited nodes separately from the distance to each node.