Permalink
Browse files

Merge branch 'canny-and-chamfer' of github.com:deads/libcvd into cann…

…y-and-chamfer
  • Loading branch information...
2 parents 80e99bb + 1717b04 commit 512f1ca730ecba3de4a562ae4c1c1566bc288b2d Edward Rosten committed Jan 8, 2013
Showing with 391 additions and 0 deletions.
  1. +187 −0 cvd_src/chamfer.hpp
  2. +204 −0 cvd_src/distance_transform.hpp
View
@@ -0,0 +1,187 @@
+#ifndef GRANTA_CHAMFER_HPP
+#define GRANTA_CHAMFER_HPP
+
+#include <cmath>
+#include <vector>
+#include <cvd/image.h>
+#include <cvd/image_io.h>
+#include <cvd/vector_image_ref.h>
+#include <TooN/TooN.h>
+#include <iostream>
+
+#include "edge.hpp"
+
+namespace Granta {
+
+ using namespace std;
+ using namespace CVD;
+
+ /*
+ Represents a contour fragment for the object model.
+ */
+ class ContourTemplate {
+ private:
+ std::vector<TooN::Vector<2> > edgels;
+ std::vector<float> directions;
+
+ public:
+ ContourTemplate() { }
+
+ ContourTemplate(std::vector<TooN::Vector<2> > &edgels, std::vector<float> &directions) : edgels(edgels), directions(directions) {
+ GRANTA_ASSERT(edgels.size() == directions.size());
+ }
+
+ void add_edgel(const TooN::Vector<2> &edgel, const double dir = 0) {
+ edgels.push_back(edgel);
+ directions.push_back(dir);
+ }
+
+ size_t size() const {
+ GRANTA_ASSERT(edgels.size() == directions.size());
+ return edgels.size();
+ }
+
+ TooN::Vector<2> &get_edgel(const size_t i) {
+ return edgels[i];
+ }
+
+ const TooN::Vector<2> &get_edgel(const size_t i) const {
+ return edgels[i];
+ }
+
+ float &get_direction(const size_t i) {
+ return directions[i];
+ }
+
+ const float &get_direction(const size_t i) const {
+ return directions[i];
+ }
+
+ /*
+ Computes all the edgels of all the contours in a masked image
+ using a Canny edge detector. The edgels are rescaled [0,1] but
+ the aspect ratio of the bounding box that encloses all edgels is
+ preserved. Finally, all edgels are shifted so that their centroid
+ lies at (0,0), making the final scale [-1,1]
+ */
+ template <class T>
+ void compute_whole(const Image<T> &image, const Image<byte> &mask) {
+ Image <T> mag(image.size());
+ const double pi = std::atan(1.0)*4;
+ canny2(image, mag);
+ img_save(mag, "Test-compute-canny.txt");
+ img_save(mask, "Test-compute-mask.png");
+ int max_x(0);
+ int max_y(0);
+ int min_x(image.size().x-1);
+ int min_y(image.size().y-1);
+ int sum_x(0);
+ int sum_y(0);
+ int n(0);
+ for (int y = 1; y < image.size().y-1; y++) {
+ for (int x = 1; x < image.size().x-1; x++) {
+ n++;
+ if (mask[y][x] && mag[y][x] > 0) {
+ const double ygrad = ((image[y-1][x-1] + 2 * image[y-1][x] + image[y-1][x+1]) -
+ (image[y+1][x-1] + 2 * image[y+1][x] + image[y+1][x+1]));
+ const double xgrad = ((image[y-1][x-1] + 2 * image[y][x-1] + image[y+1][x-1]) -
+ (image[y-1][x+1] + 2 * image[y][x+1] + image[y+1][x+1]));
+ const double dir = std::fmod(std::atan2(xgrad, ygrad)+2*pi, pi);
+ //cout << xgrad << " " << ygrad << " " << dir;
+ TooN::Vector<2> v = TooN::makeVector((double)x, (double)y);
+ edgels.push_back(v);
+ directions.push_back(dir);
+ if (x > max_x) {
+ max_x = x;
+ }
+ if (y > max_y) {
+ max_y = y;
+ }
+ if (x < min_x) {
+ min_x = x;
+ }
+ if (y < min_y) {
+ min_y = y;
+ }
+ }
+ }
+ }
+ // Rescale to [-1,1] centered at centroid.
+ const double range_x = max_x - min_x;
+ const double range_y = max_y - min_y;
+ const double range = (range_x > range_y) ? range_x : range_y; // Preserve the aspect ratio.
+ const size_t npts(edgels.size());
+ //#if 0
+ TooN::Vector<2> minv = TooN::makeVector(min_x, min_y);
+ TooN::Vector<2> centroid = TooN::makeVector(0, 0);
+ centroid = centroid / npts;
+ centroid = (centroid - minv) / range;
+ for (size_t i = 0; i < npts; i++) {
+ edgels[i] = 2.0 * (edgels[i] - minv) / range;
+ centroid += edgels[i];
+ }
+ centroid /= npts;
+ for (size_t i = 0; i < npts; i++) {
+ edgels[i] -= centroid;
+ }
+ cout << "Centroid: " << centroid;
+ //#endif
+ }
+
+ /*
+ Computes the chamfer distance of a template edge/contour given a binary feature. A pre-computed
+ distance transform must be provided.
+
+ @param DTE The Euclidean distance transform of the feature map.
+ @param chamfer_out An image to store the chamfer distances for each feature pixel.
+ @param prop The proportion of pixels to uniformly sample from.
+ @param scale The scaling factor of the template (template edgels are [0,1] normalized).
+ */
+ template <class T>
+ void compute_chamfer_distance(const Image <T> &DTE, Image <T> &chamfer_out,
+ const double prop, const double scale) {
+ GRANTA_ASSERT(edgels.size() > 0);
+ const double nedgels = static_cast<double>(edgels.size());
+ const double nsampled_edgels = (nedgels * prop);
+ const int step_size = static_cast<int>(nedgels/(nsampled_edgels));
+ for (int y = 0; y < DTE.size().y; y++) {
+ for (int x = 0; x < DTE.size().x; x++) {
+ TooN::Vector<2> v(vec(ImageRef(x, y)));
+ double total = 0.0;
+ for (int k = 0; k < edgels.size(); y += step_size) {
+ TooN::Vector<2> scaled_edgel(edgels[k] * scale);
+ ImageRef vt_plus_v = ir_rounded(scaled_edgel + v);
+ total += DTE[vt_plus_v];
+ }
+ total /= nsampled_edgels;
+ chamfer_out[y][x] = total;
+ }
+ }
+ }
+ };
+
+ ostream &operator <<(ostream &out, const ContourTemplate &contour) {
+ out << contour.size() << endl;
+ for (size_t i = 0; i < contour.size(); i++) {
+ out << contour.get_edgel(i) << " ";
+ out << contour.get_direction(i) << " ";
+ }
+ out << endl;
+ return out;
+ }
+
+ istream &operator >>(istream &in, ContourTemplate &contour) {
+ int num_contours(0);
+ in >> num_contours;
+ for (size_t i = 0; i < contour.size(); i++) {
+ TooN::Vector<2> v;
+ double dir;
+ in >> v;
+ in >> dir;
+ contour.add_edgel(v, dir);
+ }
+ return in;
+ }
+}
+
+#endif
@@ -0,0 +1,204 @@
+#ifndef GRANTA_DISTANCE_TRANSFORM_HPP
+#define GRANTA_DISTANCE_TRANSFORM_HPP
+
+#define INF 1E20
+
+#include "debug.hpp"
+#include <cvd/image.h>
+#include <vector>
+#include <limits>
+#include <algorithm>
+#include <cmath>
+
+namespace Granta {
+
+ using namespace CVD;
+ using namespace std;
+
+
+ template <class Precision = double>
+ class DistanceTransformEuclidean {
+ public:
+ DistanceTransformEuclidean() : sz(ImageRef(-1,-1)) {}
+
+ private:
+ inline void transform_row(const int n) {
+ v[0] = 0; //std::numeric_limits<Precision>::infinity();
+ z[0] = -INF; //std::numeric_limits<Precision>::infinity();
+ z[1] = +INF; // 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] = +INF; //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(Image <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++) {
+ f[x] = DT[y][x];
+ }
+ transform_row(img_sz.x);
+ for (int x = 0; x < img_sz.x; x++) {
+ DT[y][x] = d[x];
+ }
+ }
+ //#endif
+ }
+
+ template <class T>
+ void transform_image_with_ADT(const Image <T> &feature, Image <Precision> &DT, Image <ImageRef> &ADT, T onval) {
+ const ImageRef img_sz(feature.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++) {
+ f[y] = DT[y][x];
+ }
+ transform_row(img_sz.y);
+ for (int y = 0; y < img_sz.y; y++) {
+ DT[y][x] = d[y];
+ }
+ }
+ 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];
+ 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) && feature[candA] == onval) {
+ ADT[y][x] = candA;
+ }
+ else if (DT.in_image(candB) && feature[candB] == onval) {
+ ADT[y][x] = candB;
+ }
+ else if (DT.in_image(candC) && feature[candC] == onval) {
+ ADT[y][x] = candC;
+ }
+ else if (DT.in_image(candD) && feature[candD] == onval) {
+ 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");
+ }
+ }
+ }
+ }
+
+ 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);
+ }
+ }
+
+ template <class T, class Q>
+ void find_onvals(const Image<T> &feature, const Q onval, Image <Precision> &out) {
+ for (int y = 0; y < out.size().y; y++) {
+ for (int x = 0; x < out.size().x; x++) {
+ if (feature[y][x] == onval) {
+ out[y][x] = 0;
+ }
+ else {
+ out[y][x] = INF;
+ }
+ }
+ }
+ }
+
+ public:
+ template <class T, class Q>
+ void transform_ADT(const Image <T> &feature, Image <Precision> &DT, Image <ImageRef> &ADT, const Q onval) {
+ GRANTA_ASSERT(feature.size() == DT.size());
+ GRANTA_ASSERT(feature.size() == ADT.size());
+ resize(feature.size());
+ find_onvals(feature, onval, DT);
+ transform_image_with_ADT(feature, DT, ADT, onval);
+ 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]);
+ }
+ }
+ }
+
+ template <class T, class Q>
+ void transform(const Image<T> &feature, const Q onval, Image <Precision> &out) {
+ GRANTA_ASSERT(feature.size() == out.size());
+ resize(feature.size());
+ //DistanceTransformPreProcess<T, Precision> preprocess;
+ //preprocess(feature, onval, out);
+ find_onvals(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] << ",";
+ }
+ }
+ }
+
+ private:
+ ImageRef sz;
+ std::vector <Precision> d;
+ std::vector <int> v;
+ std::vector <Precision> z;
+ std::vector <Precision> f;
+ std::vector <Precision> pos;
+ };
+
+ template <class T, class Q>
+ void euclidean_distance_transform(const Image<T> &in, const Image<Q> &out) {
+ DistanceTransformEuclidean<Q> dt;
+ dt.transform(in, out);
+ }
+};
+
+#endif

0 comments on commit 512f1ca

Please sign in to comment.