Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

Cleanum and document.

  • Loading branch information...
commit b2db6928d6cdc7bf6678924fcee8a045caf5700b 1 parent 2f5a8e0
Edward Rosten authored
View
2  Makefile.in
@@ -380,7 +380,7 @@ tests/%.out: tests/%.test
.PHONY: examples
-EXAMPLES=colourmaps
+EXAMPLES=colourmaps distance_transform
EXAMPLE_PROGS=$(patsubst %,examples/%, $(EXAMPLES))
View
361 cvd/distance_transform.h
@@ -12,221 +12,209 @@
namespace CVD {
-
template <class Precision = double>
- class DistanceTransformEuclidean {
-
- // Implements
- // Distance Transforms of Sampled Functions
- // Pedro F. Felzenszwalb Daniel P. Huttenlocher
-
-
-
- public:
- DistanceTransformEuclidean()
- :sz(ImageRef(-1,-1)),
- big_number(1e9) //Hmm, why doesn't HUGE_VAL work?
- //Anyway, hilariously small number here so it works with int, too.
- {}
-
- private:
- inline void transform_row(const int n) {
- v[0] = 0; //std::numeric_limits<Precision>::infinity();
- z[0] = -big_number; //std::numeric_limits<Precision>::infinity();
- z[1] = +big_number; // std::numeric_limits<Precision>::infinity();
- int k = 0;
- for (int q = 1; q < n; q++) {
- Precision s = ((f[q]+(q*q))-(f[v[k]]+(v[k]*v[k])))/(2*q-2*v[k]);
- while (s <= z[k]) {
- k--;
- s = ((f[q]+(q*q))-(f[v[k]]+(v[k]*v[k])))/(2*q-2*v[k]);
- }
- k++;
- v[k] = q;
- z[k] = s;
- z[k+1] = +big_number; //std::numeric_limits<Precision>::infinity();
+ class DistanceTransformEuclidean
+ {
+ // Implements
+ // Distance Transforms of Sampled Functions
+ // Pedro F. Felzenszwalb Daniel P. Huttenlocher
+
+
+
+ public:
+ DistanceTransformEuclidean()
+ :sz(ImageRef(-1,-1)),
+ big_number(1e9) //Hmm, why doesn't HUGE_VAL work?
+ //Anyway, hilariously small number here so it works with int, too.
+ {}
+
+ private:
+ inline void transform_row(const int n) {
+ v[0] = 0; //std::numeric_limits<Precision>::infinity();
+ z[0] = -big_number; //std::numeric_limits<Precision>::infinity();
+ z[1] = +big_number; // std::numeric_limits<Precision>::infinity();
+ int k = 0;
+ for (int q = 1; q < n; q++) {
+ Precision s = ((f[q]+(q*q))-(f[v[k]]+(v[k]*v[k])))/(2*q-2*v[k]);
+ while (s <= z[k]) {
+ k--;
+ s = ((f[q]+(q*q))-(f[v[k]]+(v[k]*v[k])))/(2*q-2*v[k]);
}
- k = 0;
- for (int q = 0; q < n; q++) {
- while (z[k+1] < q) {
- k++;
- }
- d[q] = ((q - v[k]) * (q - v[k])) + f[v[k]];
- pos[q] = q > v[k] ? (q - v[k]) : (v[k] - q);
- //cout << "n=" << n << " q=" << q << " d[q]=" << d[q] << " v[k]=" << v[k] << " f[v[k]]=" << f[v[k]] << endl;
+ k++;
+ v[k] = q;
+ z[k] = s;
+ z[k+1] = +big_number; //std::numeric_limits<Precision>::infinity();
+ }
+ k = 0;
+ for (int q = 0; q < n; q++) {
+ while (z[k+1] < q) {
+ k++;
}
+ d[q] = ((q - v[k]) * (q - v[k])) + f[v[k]];
+ pos[q] = q > v[k] ? (q - v[k]) : (v[k] - q);
+ //cout << "n=" << n << " q=" << q << " d[q]=" << d[q] << " v[k]=" << v[k] << " f[v[k]]=" << f[v[k]] << endl;
}
+ }
- void transform_image(SubImage<Precision> &DT) {
- const ImageRef img_sz(DT.size());
+ void transform_image(SubImage<Precision> &DT) {
+ const ImageRef img_sz(DT.size());
+ for (int x = 0; x < img_sz.x; x++) {
+ for (int y = 0; y < img_sz.y; y++) {
+ f[y] = DT[y][x];
+ }
+ transform_row(img_sz.y);
+ for (int y = 0; y < img_sz.y; y++) {
+ DT[y][x] = d[y];
+ }
+ }
+ //#if 0
+ for (int y = 0; y < img_sz.y; y++) {
for (int x = 0; x < img_sz.x; x++) {
- for (int y = 0; y < img_sz.y; y++) {
- f[y] = DT[y][x];
- }
- transform_row(img_sz.y);
- for (int y = 0; y < img_sz.y; y++) {
- DT[y][x] = d[y];
- }
+ f[x] = DT[y][x];
}
- //#if 0
+ transform_row(img_sz.x);
+ for (int x = 0; x < img_sz.x; x++) {
+ DT[y][x] = d[x];
+ }
+ }
+ //#endif
+ }
+
+ template <class Functor>
+ /// Perform distance transform with reverse lookup.
+ /// @param ADT which edge pixel is closest to this one
+ void transform_image_with_ADT(SubImage <Precision> &DT, SubImage <ImageRef> &ADT, const Functor& func) {
+ const ImageRef img_sz(DT.size());
+ const double maxdist = img_sz.x * img_sz.y;
+ for (int x = 0; x < img_sz.x; x++) {
for (int y = 0; y < img_sz.y; y++) {
- for (int x = 0; x < img_sz.x; x++) {
- f[x] = DT[y][x];
- }
- transform_row(img_sz.x);
- for (int x = 0; x < img_sz.x; x++) {
- DT[y][x] = d[x];
- }
+ f[y] = DT[y][x];
+ }
+ transform_row(img_sz.y);
+ for (int y = 0; y < img_sz.y; y++) {
+ DT[y][x] = d[y];
}
- //#endif
}
-
- template <class Functor>
- /// Perform distance transform with reverse lookup.
- /// @param ADT which edge pixel is closest to this one
- void transform_image_with_ADT(SubImage <Precision> &DT, SubImage <ImageRef> &ADT, const Functor& func) {
- const ImageRef img_sz(DT.size());
- const double maxdist = img_sz.x * img_sz.y;
+ for (int y = 0; y < img_sz.y; y++) {
for (int x = 0; x < img_sz.x; x++) {
- for (int y = 0; y < img_sz.y; y++) {
- f[y] = DT[y][x];
+ f[x] = DT[y][x];
+ }
+ transform_row(img_sz.x);
+ for (int x = 0; x < img_sz.x; x++) {
+ DT[y][x] = d[x];
+ const double hyp = d[x];
+ if (hyp >= maxdist) {
+ ADT[y][x] = ImageRef(-1, -1);
+ continue;
}
- transform_row(img_sz.y);
- for (int y = 0; y < img_sz.y; y++) {
- DT[y][x] = d[y];
+ const double dx = pos[x];
+ const double dy = sqrt(hyp - dx * dx);
+ const int ddy = dy;
+ const ImageRef candA(x - dx, y - ddy);
+ const ImageRef candB(x - dx, y + ddy);
+ const ImageRef candC(x + dx, y - ddy);
+ const ImageRef candD(x + dx, y + ddy);
+ /** cerr << "hyp=" << hyp << " dx="<< dx << " dy=" << dy << " ddy=" << ddy
+ << " A=" << candA << " B=" << candB << " C=" << candC << " D=" << candD << endl;*/
+ if (DT.in_image(candA) && func(candA)) {
+ ADT[y][x] = candA;
}
- }
- for (int y = 0; y < img_sz.y; y++) {
- for (int x = 0; x < img_sz.x; x++) {
- f[x] = DT[y][x];
+ else if (DT.in_image(candB) && func(candB)) {
+ ADT[y][x] = candB;
+ }
+ else if (DT.in_image(candC) && func(candC)) {
+ ADT[y][x] = candC;
+ }
+ else if (DT.in_image(candD) && func(candD)) {
+ ADT[y][x] = candD;
}
- transform_row(img_sz.x);
- for (int x = 0; x < img_sz.x; x++) {
- DT[y][x] = d[x];
- const double hyp = d[x];
- if (hyp >= maxdist) {
- ADT[y][x] = ImageRef(-1, -1);
- continue;
- }
- const double dx = pos[x];
- const double dy = sqrt(hyp - dx * dx);
- const int ddy = dy;
- const ImageRef candA(x - dx, y - ddy);
- const ImageRef candB(x - dx, y + ddy);
- const ImageRef candC(x + dx, y - ddy);
- const ImageRef candD(x + dx, y + ddy);
- /** cerr << "hyp=" << hyp << " dx="<< dx << " dy=" << dy << " ddy=" << ddy
- << " A=" << candA << " B=" << candB << " C=" << candC << " D=" << candD << endl;*/
- if (DT.in_image(candA) && func(candA)) {
- ADT[y][x] = candA;
- }
- else if (DT.in_image(candB) && func(candB)) {
- ADT[y][x] = candB;
- }
- else if (DT.in_image(candC) && func(candC)) {
- ADT[y][x] = candC;
- }
- else if (DT.in_image(candD) && func(candD)) {
- ADT[y][x] = candD;
- }
- else {
- //ADT[y][x] = ImageRef(-1,-1);
- /**cerr << hyp << " " << ImageRef(x, y);*/
- throw std::string("feature with onval set not found");
- }
+ else {
+ //ADT[y][x] = ImageRef(-1,-1);
+ /**cerr << hyp << " " << ImageRef(x, y);*/
+ throw Exceptions::Vision::BadInput("DistanceTransformEuclidean: no points set");
}
}
}
-
- void resize(const ImageRef &new_size) {
- if (new_size != sz) {
- sz = new_size;
- int m(std::max(new_size.x, new_size.y));
- d.resize(m);
- v.resize(m);
- f.resize(m);
- pos.resize(m);
- z.resize(m+1);
- }
+ }
+
+ void resize(const ImageRef &new_size) {
+ if (new_size != sz) {
+ sz = new_size;
+ int m(std::max(new_size.x, new_size.y));
+ d.resize(m);
+ v.resize(m);
+ f.resize(m);
+ pos.resize(m);
+ z.resize(m+1);
}
+ }
- public:
- template <class T>
- void transform_ADT(const SubImage<T>& feature, SubImage<Precision> &DT, SubImage<ImageRef> &ADT) {
- if(feature.size() != DT.size())
- throw Exceptions::Vision::IncompatibleImageSizes(__FUNCTION__);
-
- if(feature.size() != ADT.size())
- throw Exceptions::Vision::IncompatibleImageSizes(__FUNCTION__);
+ public:
+ template <class T>
+ void transform_ADT(const SubImage<T>& feature, SubImage<Precision> &DT, SubImage<ImageRef> &ADT) {
+ if(feature.size() != DT.size())
+ throw Exceptions::Vision::IncompatibleImageSizes(__FUNCTION__);
+
+ if(feature.size() != ADT.size())
+ throw Exceptions::Vision::IncompatibleImageSizes(__FUNCTION__);
- resize(feature.size());
+ resize(feature.size());
- apply_functor(DT, NotZero<T>(feature));
+ apply_functor(DT, NotZero<T>(feature));
- transform_image_with_ADT(DT, ADT, NotZero<T>(feature));
- /*for (int y = 0; y < DT.size().y; y++) {
- for (int x = 0; x < DT.size().x; x++) {
- DT[y][x] = sqrt(DT[y][x]);
- }
- }*/
- }
+ transform_image_with_ADT(DT, ADT, NotZero<T>(feature));
+ }
- void transform(SubImage<Precision> &out) {
- resize(out.size());
- //DistanceTransformPreProcess<T, Precision> preprocess;
- //preprocess(feature, onval, out);
- transform_image(out);
- /*for (int y = 0; y < out.size().y; y++) {
- for (int x = 0; x < out.size().x; x++) {
- out[y][x] = sqrt(out[y][x]);
- //cout << out[y][x] << ",";
- }
- }*/
- }
+ void transform(SubImage<Precision> &out) {
+ resize(out.size());
+ transform_image(out);
+ }
- template<class C>
- struct NotZero
- {
- NotZero(const SubImage<C>& im_)
- :im(im_)
- {}
-
- const SubImage<C>& im;
- bool operator()(const ImageRef& i) const
- {
- return im[i] != 0;
- }
- };
-
- template<class Out, class Functor>
- void apply_functor(SubImage<Out>& out, const Functor& f)
+ template<class C>
+ struct NotZero
+ {
+ NotZero(const SubImage<C>& im_)
+ :im(im_)
+ {}
+
+ const SubImage<C>& im;
+ bool operator()(const ImageRef& i) const
{
- for (int y = 0; y < out.size().y; y++)
- for (int x = 0; x < out.size().x; x++)
- if (f(ImageRef(x, y)))
- out[y][x] = 0;
- else
- out[y][x] = big_number;
+ return im[i] != 0;
}
-
- private:
- ImageRef sz;
- double big_number;
- std::vector <Precision> d;
- std::vector <int> v;
- std::vector <Precision> z;
- std::vector <Precision> f;
- std::vector <Precision> pos;
- };
+ };
+
+ template<class Out, class Functor>
+ void apply_functor(SubImage<Out>& out, const Functor& f)
+ {
+ for (int y = 0; y < out.size().y; y++)
+ for (int x = 0; x < out.size().x; x++)
+ if (f(ImageRef(x, y)))
+ out[y][x] = 0;
+ else
+ out[y][x] = big_number;
+ }
+
+ private:
+ ImageRef sz;
+ double big_number;
+ std::vector <Precision> d;
+ std::vector <int> v;
+ std::vector <Precision> z;
+ std::vector <Precision> f;
+ std::vector <Precision> pos;
+ };
///Compute squared Euclidean distance transform using the Felzenszwalb & Huttenlocher algorithm.
+ ///Example in examples/distance_transform.cc
///@ingroup gVision
///@param in input image: thresholded so anything &gt; 0 is on the object
///@param out output image is euclidean distance of input image.
+ ///@throws Exceptions::Vision::BadInput Throws if the input contains no points.
template <class T, class Q>
void euclidean_distance_transform_sq(const SubImage<T> &in, SubImage<Q> &out) {
if(in.size() != out.size())
@@ -238,9 +226,12 @@ namespace CVD {
///Compute squared Euclidean distance transform using the Felzenszwalb & Huttenlocher algorithm.
+ ///Example in examples/distance_transform.cc
///@ingroup gVision
///@param in input image: thresholded so anything &gt; 0 is on the object
///@param out output image is euclidean distance of input image.
+ ///@param lookup_DT For each output pixel, this is the location of the closest input pixel.
+ ///@throws Exceptions::Vision::BadInput Throws if the input contains no points.
template<class T, class Q>
void euclidean_distance_transform_sq(const SubImage<T> &in, SubImage<Q> &out, SubImage<ImageRef>& lookup_DT) {
if(in.size() != out.size())
@@ -250,9 +241,11 @@ namespace CVD {
}
///Compute Euclidean distance transform using the Felzenszwalb & Huttenlocher algorithm.
+ ///Example in examples/distance_transform.cc
///@ingroup gVision
///@param in input image: thresholded so anything &gt; 0 is on the object
///@param out output image is euclidean distance of input image.
+ ///@throws Exceptions::Vision::BadInput Throws if the input contains no points.
template<class T, class Q>
void euclidean_distance_transform(const SubImage<T> &in, SubImage<Q> &out) {
euclidean_distance_transform_sq(in, out);
@@ -263,9 +256,12 @@ namespace CVD {
///Compute Euclidean distance transform using the Felzenszwalb & Huttenlocher algorithm.
+ ///Example in examples/distance_transform.cc
///@ingroup gVision
///@param in input image: thresholded so anything &gt; 0 is on the object
///@param out output image is euclidean distance of input image.
+ ///@param lookup_DT For each output pixel, this is the location of the closest input pixel.
+ ///@throws Exceptions::Vision::BadInput Throws if the input contains no points.
template<class T, class Q>
void euclidean_distance_transform(const SubImage<T> &in, SubImage<Q> &out, SubImage<ImageRef>& lookup_DT) {
euclidean_distance_transform_sq(in, out, lookup_DT);
@@ -308,8 +304,13 @@ namespace CVD {
///@ingroup gVision
///@param in input image: thresholded so anything &gt; 0 is on the object
///@returns output image is euclidean distance of input image.
+ ///@throws Exceptions::Vision::BadInput Throws if the input contains no points.
template <class T>
Image euclidean_distance_transform(const SubImage<T> &in);
#endif
+
+
+ ///@example distance_transform.cc
+ ///Example of how to use CVD::euclidean_distance_transform_sq
};
#endif
View
8 cvd/vision_exceptions.h
@@ -31,6 +31,14 @@ namespace CVD {
what = "Input ImageRefs not in image in " + function;
};
};
+
+
+ struct BadInput: public All{
+ BadInput(const std::string& function)
+ {
+ what = "Bad input in " + function;
+ };
+ };
};
}
}
View
49 examples/distance_transform.cc
@@ -0,0 +1,49 @@
+#include <cvd/distance_transform.h>
+#include <cvd/random.h>
+#include <cvd/image_io.h>
+#include <algorithm>
+
+using namespace CVD;
+using namespace std;
+
+int main()
+{
+ //Create a blank image.
+ Image<byte> im(ImageRef(128, 128), 0);
+
+ //Scatter down 7 points at random.
+ for(int i=1; i < 8; i++)
+ im[rand() % im.size().y][rand() % im.size().x] = i;
+
+ Image<int> dt(im.size());
+ Image<ImageRef> inverse_dt(im.size());
+
+ //Perfom the distance transform
+ euclidean_distance_transform_sq(im, dt, inverse_dt);
+
+ //Create an output which is the distance transfom of the input,
+ //but coloured according to which pixel is closest.
+ int largest_distance = *max_element(dt.begin(), dt.end());
+
+ Image<Rgb<byte> > out(im.size());
+
+ for(int y=0; y < im.size().y; y++)
+ for(int x=0; x < im.size().x; x++)
+ {
+ int c = floor(sqrt(dt[y][x]*1.0/largest_distance) * 255 + .5);
+
+ Rgb<byte> r(0,0,0);
+ if(im[inverse_dt[y][x]]&1)
+ r.red = c;
+ if(im[inverse_dt[y][x]]&2)
+ r.green = c;
+ if(im[inverse_dt[y][x]]&4)
+ r.blue = c;
+
+ out[y][x] = r;
+ }
+
+
+ img_save(out, "distance_transform_result.png");
+
+}
Please sign in to comment.
Something went wrong with that request. Please try again.