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;
}
if (newY <= floor(image.rows/2) && (newX >= source.x-newY && newX <= source.x+newY) || (newY > floor(image.rows/2) && (newX >=source.x-(goal.y-newY) && newX <=source.x+(goal.y-newY))))
is a nightmare to parse. \$\endgroup\$ – Emily L. Aug 7 '15 at 13:31