https://doi.org/10.5201/ipol.2021.296
skeleton.h
/*
Copyright (c) 2020 Yuchen He <yhe306@gatech.edu> and Luis Alvarez <lalvarez@ulpgc.es>
This program is free software: you can redistribute it and/or modify it under
the terms of the GNU Affero 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 Affero General Public License for more
details.
You should have received a copy of the GNU Affero General Public License along
with this program. If not, see <http://www.gnu.org/licenses/>.
*/
/**
* \file skeleton.h
* \brief functions needed for HJS
* \author Yuchen He and Luis Alvarez
*/
#ifndef SKELETON_H
#define SKELETON_H
#include <cmath>
#include "auxiliary/image2D.h"
/// HJS ALGORITHM
/// DISTANCE TRANSFORMATIONS
template <class T>
image2D<float> distance_function2D_contour(const image2D<T> &Is, vector <point2d> &contour); /* brute_force */
template <class T>
image2D<float> distance_function2D_fast_sweeping(const image2D<T> &Is, vector <point2d> &contour); /* fast_sweeping */
template <class T>
image2D<float> distance_function2D_F_H(const image2D<T> &Is, vector <point2d> &contour);
/// F-H related codes
// compute the slope of f between p and q
static float obtain_slope(float *f, int q, int p);
// compute the squared distances to a binary 1-D image
static void squared_distances_1d(float *d, float *f, int n);
// compute the squared distances to 2-D image
static void squared_distances_2d(float *f, int w, int h);
void squared_euclidean_distance_to_nonzeros(float *x, int w,int h);
void build_positive_distance(float *d, float *m, int w, int h);
/// FLUX COMPUTATION USING IJCV PAPER "Hamilton-Jacobi Skeletons"
image2D<float> flux(image2D<float> &Ix, image2D<float> &Iy);
/// HOMOTOPY PRESERVING THINNING USING IJCV PAPER "Hamilton-Jacobi Skeletons"
template <class T>
vector<pair<point2d,float>> homotopy_thinning(image2D<T> &Is, const image2D<float> &D, image2D<float> &flux_field, vector<point2d> contour, float flux_threshold);
template <class T>
bool is_simple(image2D<T> &Is, int &index, vector<int> &N);
template <class T>
bool is_end_point(image2D<T> &Is, int &index, vector<int> &N);
/// RECONSTRUCT SHAPE FROM UNION OF DISKS
template<class T>
void shape_from_skeleton(image2D<T> &canvas, const vector<pair<point2d,float>> &medial_axis, const float &epsilon);
//////////////////////////////////////////////////////////////////////////////////////////////
/**
* \brief brute_force method to compute distance map
*/
/// COMPUTATION OF 2D SIGNED DISTANCE FUNCTION FROM AN IMAGE AS REGION INDICATOR
template <class T>
image2D<float> distance_function2D_contour(
const image2D<T> &Is, /// REGION INDICATOR Is[m]==0 MEANS VOXEL m IS INSIDE DE REGION, Is[m]==1 OUTSIDE THE REGION
vector <point2d> &contour) /// 2D POINTS DEFINING THE REGION CONTOUR (IF ITS SIZE IS 0, IT IS COMPUTED FROM THE REGION INDICATOR Is
{
/// AUXILIARY VARIABLES
int width=Is.width();
int height=Is.height();
image2D<float> u(width,height,-1.);
/// DEFINITION OF VOXEL NEIGBORHOOD AND DISTANCE
vector<int> N=u.Neighborhood();
/// COMPUTATION OF THE SILHOUETTE BOUNDARY IF REQUIRED
if(contour.size()==0){
for(int x=1;x<width-1;x++){
for(int y=1;y<height-1;y++){
int m=x+y*width;
if(Is[m]==0.0){
for(size_t k=0;k<N.size();k++){
if(Is[m+N[k]]==1.0){
u[m]=0.;
contour.push_back(point2d((float) x,(float) y));
break;
}
}
}
}
}
}
/// COMPUTE THE DISTANCE TO THE CONTOUR BY DEFINITION
#ifdef _OPENMP
#pragma omp parallel for
#endif
for(int x=0;x<width;x++){
for(int y=0;y<height;y++){
int m=x+y*width;
if(Is[m]==1.0|| u[m]>=0 ) continue;
float min_error=width*width+height*height;
for(size_t n=0;n<contour.size();n++){
float error= (float) (x-contour[n].x)*(x-contour[n].x)+(y-contour[n].y)*(y-contour[n].y);
if(error<min_error) min_error=error;
}
u[m]=sqrtf(min_error);
}
}
return u;
}
//////////////////////////////////////////////////////////////////////////////////////////////
/**
* \brief fast_sweeping method to compute distance map
*/
/// COMPUTATION OF 2D SIGNED DISTANCE FUNCTION FROM AN IMAGE AS REGION INDICATOR
template <class T>
image2D<float> distance_function2D_fast_sweeping(
const image2D<T> &Is, /// REGION INDICATOR Is[m]==0 MEANS VOXEL m IS INSIDE DE REGION, Is[m]==1 OUTSIDE THE REGION
vector <point2d> &contour) /// 2D POINTS DEFINING THE REGION CONTOUR (IF ITS SIZE IS 0, IT IS COMPUTED FROM THE REGION INDICATOR Is
{
/// AUXILIARY VARIABLES
int width=Is.width();
int height=Is.height();
/// INITIALIZE THE VOXELS WITH LARGE POSITIVE VALUES
image2D<float> u(width,height, (float) width*width+height*height);
/// DEFINITION OF VOXEL NEIGBORHOOD AND DISTANCE
vector<int> N=u.Neighborhood();
/// COMPUTATION OF THE SILHOUETTE BOUNDARY IF REQUIRED
/// AND THE DISTANCE ON THE BOUNDARY IS UPDATED TO 0
if(contour.size()==0){
for(int x=1;x<width-1;x++){
for(int y=1;y<height-1;y++){
int m=x+y*width;
if(Is[m]==0.0){
for(size_t k=0;k<N.size();k++){
if(Is[m+N[k]]==1.0){
u[m]=0.;
contour.push_back(point2d((float) x,(float) y));
break;
}
}
}
}
}
}
float min_x;
float min_y;
float u_new;
for(int iter=0; iter<4; iter++){
for(int x=0;x<width;x++){
for(int y=0;y<height;y++){
int m;
int x_index;
int y_index;
if(iter==0){
// forward x; forward y
x_index = x;
y_index = y;
}
else if(iter==1){
// backward x; forward y
x_index = width-x-1;
y_index = y;
}
else if(iter==2){
// backward x; backward y
x_index = width-x-1;
y_index = height-y-1;
}
else{//iter==3
// forward x; backward y
x_index = x;
y_index = height-y-1;
}
m = x_index+y_index*width;
if(x_index==0){
min_x = u[m+1];
}
else if(x_index==width-1){
min_x = u[m-1];
}
else{
min_x = (u[m-1]<u[m+1])?u[m-1]:u[m+1];
}
if(y_index==0){
min_y = u[m+width];
}
else if(y_index==height-1){
min_y = u[m-width];
}
else{
min_y = (u[m-width]<u[m+width])?u[m-width]:u[m+width];
}
if(fabs(min_x-min_y)<1){
u_new = (min_x+min_y+sqrtf(2-(min_x-min_y)*(min_x-min_y)))/2;
}
else{
u_new = (min_x<min_y?min_x:min_y)+1;
}
u[m] = u_new<u[m]?u_new:u[m];
}
}
}
for(int m =0; m<u.size(); m++){
u[m] = (Is[m]==1.0)?-1.0:u[m];
}
return u;
}
//////////////////////////////////////////////////////////////////////////////////////////////
/**
* \brief Homopotopy preserving thinning part
*/
/// HOMOTOPY THINNING ACCORDING TO THE AVERAGE FLUX COMPUTED ABOVE
template <class T>
vector<pair<point2d,float>> homotopy_thinning(image2D<T> &Is, const image2D<float> &D, image2D<float> &flux_field, vector<point2d> contour, float flux_threshold){
int image_width = flux_field.width();
vector<int> N=flux_field.Neighborhood();
using point_flux = pair<vector<int>,float>;
// Define the comparing rule for the heap
struct my_comparator
{
bool operator()(point_flux const& pf1, point_flux const& pf2) const
{
return pf1.second < pf2.second;
}
};
using my_priority_queue = std::priority_queue<point_flux, vector<point_flux>, my_comparator>;
my_priority_queue point_flux_pairs;
// Insert simple boundary points into the heap
for(size_t k=0; k<contour.size();k++){
int p_x = (int) contour[k].x;
int p_y = (int) contour[k].y;
if(p_x==0 || p_x==Is.width()-1 || p_y==0 || p_y==Is.height()-1) continue;
int index = p_x + p_y*image_width;
if(is_simple(Is, index, N)){
Is[index]= 2.0;
point_flux_pairs.push(point_flux ({p_x, p_y},flux_field[index]));
}
}
// Iteratively sifting the points
while(!point_flux_pairs.empty()){
int p_x = point_flux_pairs.top().first[0];
int p_y = point_flux_pairs.top().first[1];
float flux = point_flux_pairs.top().second;
point_flux_pairs.pop();
int index = p_x + p_y*image_width;
Is[index]= 0.0;
if(is_simple(Is, index, N)){
if((!is_end_point(Is, index, N))||(flux>flux_threshold)){
// Refill the heap with neighboring points
Is[index] = 1.0;
for(size_t m=0;m<N.size();m++){
int cur_index = index+N[m];
if(Is[cur_index]==0.0){
int p_y_new = cur_index/image_width;
int p_x_new = cur_index-p_y_new*image_width;
if(p_x_new==0 || p_x_new==Is.width()-1 || p_y_new==0 || p_y_new==Is.height()-1) continue;
if(is_simple(Is, cur_index, N)){
Is[cur_index] = 2.0;
point_flux_pairs.push(point_flux ({p_x_new, p_y_new},flux_field[cur_index]));
}
}
}
}
else
{
Is[index] = 0.0;
}
}
}
// Record skeleton point coordinates
vector<pair<point2d,float>> skeleton_point;
for(int m=0; m<Is.size(); m++){
if(Is[m] == 0.0){
int p_y = m/image_width;
int p_x = m-p_y*image_width;
skeleton_point.push_back(pair<point2d,float> (point2d((float) p_x,(float) p_y),D[m]));
}
}
return skeleton_point;
}
//////////////////////////////////////////////////////////////////////////////////////////////
/**
* \brief Subprocedure to check simple points used in homotopy_thinning
*/
/// Check simple points
template <class T>
bool is_simple(image2D<T> &Is, int &index, vector<int> &N) {
int width = Is.width();
int p_y = index/width;
int p_x = index-p_y*width;
if(p_x==0 || p_x==Is.width()-1 || p_y==0 || p_y==Is.height()-1) return true;
set<int> neighbor_vertices;
int num_vertices = 0;
int num_edges = 0;
// Reindexing the neighbors
for(int m=0;m<8;m++){
if(Is[index+N[m]]!=1.0){
switch(m){
case 0:
neighbor_vertices.insert(3);
break;
case 1:
neighbor_vertices.insert(7);
break;
case 2:
neighbor_vertices.insert(5);
break;
case 3:
neighbor_vertices.insert(1);
break;
case 4:
neighbor_vertices.insert(4);
break;
case 5:
neighbor_vertices.insert(2);
break;
case 6:
neighbor_vertices.insert(6);
break;
case 7:
neighbor_vertices.insert(0);
break;
}
}
}
// Count the edges while avoiding the 3-circles
for (auto it = neighbor_vertices.begin(); it != neighbor_vertices.end(); it++) {
int prev_vertex = (*it-1>=0)?(*it-1):(*it+7);
int post_vertex = (*it+1<=7)?(*it+1):(*it-7);
if((neighbor_vertices.find(post_vertex)!= neighbor_vertices.end()) ||
(*it==1&&(neighbor_vertices.find(3)!= neighbor_vertices.end()))||
(*it==3&&(neighbor_vertices.find(5)!= neighbor_vertices.end()))||
(*it==5&&(neighbor_vertices.find(7)!= neighbor_vertices.end()))||
(*it==7&&(neighbor_vertices.find(1)!= neighbor_vertices.end()))){
num_edges += 1;
num_vertices += 1;
}else{
num_vertices += 1;
}
if(*it==0||*it==2||*it==4||*it==6){
if((neighbor_vertices.find(prev_vertex)!= neighbor_vertices.end()) &&
(neighbor_vertices.find(post_vertex)!= neighbor_vertices.end())){
num_edges -= 1;
num_vertices -= 1;
}
}
}
return (num_vertices-num_edges == 1);
}
//////////////////////////////////////////////////////////////////////////////////////////////
/**
* \brief Subprocedure to check endpoints used in homotopy_thinning
*/
/// Check endpoints
template <class T>
bool is_end_point(image2D<T> &Is, int &index, vector<int> &N) {
vector<int> neighbor_vertices;
for(size_t m=0;m<N.size();m++){
if(Is[index+N[m]]!=1.0){
switch(m){
case 0:
neighbor_vertices.push_back(3);
break;
case 1:
neighbor_vertices.push_back(7);
break;
case 2:
neighbor_vertices.push_back(5);
break;
case 3:
neighbor_vertices.push_back(1);
break;
case 4:
neighbor_vertices.push_back(4);
break;
case 5:
neighbor_vertices.push_back(2);
break;
case 6:
neighbor_vertices.push_back(6);
break;
case 7:
neighbor_vertices.push_back(0);
break;
}
if(neighbor_vertices.size()>2) return 0;
}
}
if(neighbor_vertices.size()==1){
return 1;
}else{// neighbor_vertices.size()==2
if((abs(neighbor_vertices[0]-neighbor_vertices[1])==1)||
(neighbor_vertices[0]==7 && neighbor_vertices[1]==0)||
(neighbor_vertices[0]==0 && neighbor_vertices[1]==7)){
return 1;
}
}
return 0;
}
//////////////////////////////////////////////////////////////////////////////////////////////
/**
* \brief Reconstruct the shape from skeleton by union of disks
*/
/// shape_from_skeleton
template <class T>
void shape_from_skeleton(image2D<T> &canvas, const vector<pair<point2d,float>> &medial_axis, const float &epsilon){
int image_width = canvas.width();
for(size_t i = 0; i < medial_axis.size(); i++){
int mid_x = medial_axis[i].first.x;
int mid_y = medial_axis[i].first.y;
float medial_dist = medial_axis[i].second;
float inner_dist = medial_dist/sqrt(2);
/// WE INCREASE THE DISK RADIUS TO COMPENSATE BIAS INTRODUCED BY PIXELIZATION
medial_dist+=epsilon;
/// RECONSTRUCT SHAPE BY TAKING UNION OF DISKS
int min_x=round(mid_x-medial_dist),max_x=round(mid_x+medial_dist);
int min_y=round(mid_y-medial_dist),max_y=round(mid_y+medial_dist);
if(min_x<0) min_x=0;
if(min_y<0) min_y=0;
if(max_x>=image_width) max_x=image_width-1;
if(max_y>=canvas.height()) max_y=canvas.height()-1;
float medial_dist2=medial_dist*medial_dist;
for(int y=min_y;y<=max_y;y++){
for(int x=min_x;x<=max_x;x++){
if(fabs(x-mid_x)>inner_dist || fabs(y-mid_y)>inner_dist){ /* compute the distance outside a cube inscribed in the maximal disk */
if( pow(x - mid_x,2) + pow(y - mid_y,2) <= medial_dist2){
canvas[y*image_width+x] = 0.0;
}
}else{
canvas[y*image_width+x] = 0.0;
}
}
}
}
}
template <class T>
image2D<float> distance_function2D_F_H(const image2D<T> &Is, vector <point2d> &contour) /* F-H */
{
unsigned w = Is.width(); // image width
unsigned h = Is.height(); // image height
float *X = (float *) malloc(w*h*sizeof*X);
image2D<float> D(w,h,-1.0);
/// DEFINITION OF VOXEL NEIGBORHOOD AND DISTANCE
vector<int> N=D.Neighborhood();
/// COMPUTATION OF THE SILHOUETTE BOUNDARY IF REQUIRED
image2D<T> mod_Is = Is;
if(contour.size()==0){
for(int x=1;x<w-1;x++){
for(int y=1;y<h-1;y++){
int m=x+y*w;
if(Is[m]==0.0){
for(size_t k=0;k<N.size();k++){
if(Is[m+N[k]]==1.0){
mod_Is[m] = 1.0;
contour.push_back(point2d((float) x,(float) y));
break;
}
}
}
}
}
}
for (size_t i = 0; i < w*h; i++)
X[i] = mod_Is[i];
squared_euclidean_distance_to_nonzeros(X, w, h);
// for (size_t i = 0; i < h*w; i++)
// {
// D[i] = Is[i]==0?sqrt(X[i]):-1;
// }
for (size_t i = 0; i < h*w; i++)
{
D[i] = sqrt(X[i]);
}
// D = (D-D.min())/(D.max()-D.min());
return D;
}
#endif
