UPDATE
The code from this answer seems to not work as expected (at least I didn't manage to make it work satisfactorily).
So I run the code on which that implementation was based, that you can find here
The (more reasonble) result is:
I'll leave the code below for future references.
I adapted the mex implementation from here. The result on your image with the below code is:
I'll let you evaluate if this is helpful to you. The code is below.
swt.h
#include <opencv2\opencv.hpp> #include <vector> #include <map> #include <set> #include <algorithm> using namespace std; namespace sw { #define PI 3.14159265 struct Point2d { int x; int y; float SWT; }; struct Point2dFloat { float x; float y; }; struct Ray { Point2d p; Point2d q; std::vector<Point2d> points; }; void strokeWidthTransform(const float * edgeImage, const float * gradientX, const float * gradientY, bool dark_on_light, float * SWTImage, int h, int w, std::vector<Ray> & rays) { // First pass float prec = .05f; for (int row = 0; row < h; row++){ const float* ptr = edgeImage + row*w; for (int col = 0; col < w; col++){ if (*ptr > 0) { Ray r; Point2d p; p.x = col; p.y = row; r.p = p; std::vector<Point2d> points; points.push_back(p); float curX = (float)col + 0.5f; float curY = (float)row + 0.5f; int curPixX = col; int curPixY = row; float G_x = gradientX[col + row*w]; float G_y = gradientY[col + row*w]; // normalize gradient float mag = sqrt((G_x * G_x) + (G_y * G_y)); if (dark_on_light){ G_x = -G_x / mag; G_y = -G_y / mag; } else { G_x = G_x / mag; G_y = G_y / mag; } while (true) { curX += G_x*prec; curY += G_y*prec; if ((int)(floor(curX)) != curPixX || (int)(floor(curY)) != curPixY) { curPixX = (int)(floor(curX)); curPixY = (int)(floor(curY)); // check if pixel is outside boundary of image if (curPixX < 0 || (curPixX >= w) || curPixY < 0 || (curPixY >= h)) { break; } Point2d pnew; pnew.x = curPixX; pnew.y = curPixY; points.push_back(pnew); if (edgeImage[curPixY*w + curPixX] > 0) { r.q = pnew; // dot product float G_xt = gradientX[curPixY*w + curPixX]; float G_yt = gradientY[curPixY*w + curPixX]; mag = sqrt((G_xt * G_xt) + (G_yt * G_yt)); if (dark_on_light){ G_xt = -G_xt / mag; G_yt = -G_yt / mag; } else { G_xt = G_xt / mag; G_yt = G_yt / mag; } if (acos(G_x * -G_xt + G_y * -G_yt) < PI / 2.0) { float length = sqrt(((float)r.q.x - (float)r.p.x)*((float)r.q.x - (float)r.p.x) + ((float)r.q.y - (float)r.p.y)*((float)r.q.y - (float)r.p.y)); for (std::vector<Point2d>::iterator pit = points.begin(); pit != points.end(); pit++) { float* pSWT = SWTImage + w * pit->y + pit->x; if (*pSWT < 0) { *pSWT = length; } else { *pSWT = std::min(length, *pSWT); } } r.points = points; rays.push_back(r); } break; } } } } ptr++; } } } bool Point2dSort(const Point2d &lhs, const Point2d &rhs) { return lhs.SWT < rhs.SWT; } void SWTMedianFilter(float * SWTImage, int h, int w, std::vector<Ray> & rays, float maxWidth = -1) { for (std::vector<Ray>::iterator rit = rays.begin(); rit != rays.end(); rit++) { for (std::vector<Point2d>::iterator pit = rit->points.begin(); pit != rit->points.end(); pit++) { pit->SWT = SWTImage[w*pit->y + pit->x]; } std::sort(rit->points.begin(), rit->points.end(), &Point2dSort); //std::nth_element( rit->points.begin(), rit->points.end(), rit->points.size()/2, &Point2dSort ); float median = (rit->points[rit->points.size() / 2]).SWT; if (maxWidth > 0 && median >= maxWidth) { median = -1; } for (std::vector<Point2d>::iterator pit = rit->points.begin(); pit != rit->points.end(); pit++) { SWTImage[w*pit->y + pit->x] = std::min(pit->SWT, median); } } } typedef std::vector< std::set<int> > graph_t; // graph as a list of neighbors per node void connComp(const graph_t& g, std::vector<int>& c, int i, int l) { // starting from node i labe this conn-comp with label l if (i < 0 || i > g.size()) { return; } std::vector< int > stack; // push i stack.push_back(i); c[i] = l; while (!stack.empty()) { // pop i = stack.back(); stack.pop_back(); // go over all nieghbors for (std::set<int>::const_iterator it = g[i].begin(); it != g[i].end(); it++) { if (c[*it] < 0) { stack.push_back(*it); c[*it] = l; } } } } int findNextToLabel(const graph_t& g, const vector<int>& c) { for (int i = 0; i < c.size(); i++) { if (c[i] < 0) { return i; } } return c.size(); } int connected_components(const graph_t& g, vector<int>& c) { // check for empty graph! if (g.empty()) { return 0; } int i = 0; int num_conn = 0; do { connComp(g, c, i, num_conn); num_conn++; i = findNextToLabel(g, c); } while (i < g.size()); return num_conn; } std::vector< std::vector<Point2d> > findLegallyConnectedComponents(const float* SWTImage, int h, int w, std::vector<Ray> & rays) { std::map<int, int> Map; std::map<int, Point2d> revmap; std::vector<std::vector<Point2d> > components; // empty int num_vertices = 0, idx = 0; graph_t g; // Number vertices for graph. Associate each point with number for (int row = 0; row < h; row++){ for (int col = 0; col < w; col++){ idx = col + w * row; if (SWTImage[idx] > 0) { Map[idx] = num_vertices; Point2d p; p.x = col; p.y = row; revmap[num_vertices] = p; num_vertices++; std::set<int> empty; g.push_back(empty); } } } if (g.empty()) { return components; // nothing to do with an empty graph... } for (int row = 0; row < h; row++){ for (int col = 0; col < w; col++){ idx = col + w * row; if (SWTImage[idx] > 0) { // check pixel to the right, right-down, down, left-down int this_pixel = Map[idx]; float thisVal = SWTImage[idx]; if (col + 1 < w) { float right = SWTImage[w*row + col + 1]; if (right > 0 && (thisVal / right <= 3.0 || right / thisVal <= 3.0)) { g[this_pixel].insert(Map[w*row + col + 1]); g[Map[w*row + col + 1]].insert(this_pixel); //boost::add_edge(this_pixel, map.at(row * SWTImage->width + col + 1), g); } } if (row + 1 < h) { if (col + 1 < w) { float right_down = SWTImage[w*(row + 1) + col + 1]; if (right_down > 0 && (thisVal / right_down <= 3.0 || right_down / thisVal <= 3.0)) { g[this_pixel].insert(Map[w*(row + 1) + col + 1]); g[Map[w*(row + 1) + col + 1]].insert(this_pixel); // boost::add_edge(this_pixel, map.at((row+1) * SWTImage->width + col + 1), g); } } float down = SWTImage[w*(row + 1) + col]; if (down > 0 && (thisVal / down <= 3.0 || down / thisVal <= 3.0)) { g[this_pixel].insert(Map[w*(row + 1) + col]); g[Map[w*(row + 1) + col]].insert(this_pixel); //boost::add_edge(this_pixel, map.at((row+1) * SWTImage->width + col), g); } if (col - 1 >= 0) { float left_down = SWTImage[w*(row + 1) + col - 1]; if (left_down > 0 && (thisVal / left_down <= 3.0 || left_down / thisVal <= 3.0)) { g[this_pixel].insert(Map[w*(row + 1) + col - 1]); g[Map[w*(row + 1) + col - 1]].insert(this_pixel); //boost::add_edge(this_pixel, map.at((row+1) * SWTImage->width + col - 1), g); } } } } } } std::vector<int> c(num_vertices, -1); int num_comp = connected_components(g, c); components.reserve(num_comp); //std::cout << "Before filtering, " << num_comp << " components and " << num_vertices << " vertices" << std::endl; for (int j = 0; j < num_comp; j++) { std::vector<Point2d> tmp; components.push_back(tmp); } for (int j = 0; j < num_vertices; j++) { Point2d p = revmap[j]; (components[c[j]]).push_back(p); } return components; } enum { EIN = 0, GXIN, GYIN, DOLFIN, MAXWIN, NIN }; void swt_mex(const float* edgeImage, const float* gradientX, const float* gradientY, float* SWTImage, float* pComp, int* nstrokes, int w, int h, bool dark_on_light) { float maxWidth = w; std::vector<Ray> rays; strokeWidthTransform(edgeImage, gradientX, gradientY, dark_on_light, SWTImage, h, w, rays); SWTMedianFilter(SWTImage, h, w, rays, maxWidth); std::vector<std::vector<Point2d> > components = findLegallyConnectedComponents(SWTImage, h, w, rays); *nstrokes = components.size(); for (int ci = 0; ci < components.size(); ci++) { for (std::vector<Point2d>::iterator it = components[ci].begin(); it != components[ci].end(); it++) { pComp[w * it->y + it->x] = ci + 1; } } } void swt(const cv::Mat1b& img, cv::Mat1f& strokes, int* nstrokes, bool dark_on_light = true) { cv::Mat1b edgeMap; cv::Canny(img, edgeMap, 400, 200); cv::Mat1f floatEdgeMap; edgeMap.convertTo(floatEdgeMap, CV_32F); cv::Mat1b blurred; cv::GaussianBlur(img, blurred, cv::Size(5, 5), 0.3*(2.5 - 1) + .8); cv::Mat1f gx, gy; cv::Sobel(blurred, gx, CV_32F, 1, 0); cv::Sobel(blurred, gy, CV_32F, 0, 1); cv::medianBlur(gx, gx, 3); cv::medianBlur(gy, gy, 3); cv::Mat1f swtimg(img.rows, img.cols, -1.f); strokes = cv::Mat1f(img.rows, img.cols, 0.f); swt_mex((float*)floatEdgeMap.data, (float*)gx.data, (float*)gy.data, (float*)swtimg.data, (float*)strokes.data, nstrokes, img.cols, img.rows, dark_on_light); } }
main
#include <opencv2/opencv.hpp> #include "swt.h" using namespace cv; int main(int, char** argv) { Mat1b img = cv::imread("path_to_image", IMREAD_GRAYSCALE); // Compute SWT Mat1f strokes; int nstrokes; sw::swt(img, strokes, &nstrokes); // Create color table vector<Vec3b> colors(nstrokes+1); colors[0] = Vec3b(0, 0, 0); RNG rng; for (int i = 0; i < nstrokes; ++i) { colors[i + 1] = Vec3b(rng.uniform(0, 255), rng.uniform(0, 255), rng.uniform(0, 255)); } // Colors strokes Mat3b coloredStrokes(strokes.size(), Vec3b(0,0,0)); for (int r = 0; r < strokes.rows; ++r) { for (int c = 0; c < strokes.cols; ++c) { coloredStrokes(r, c) = colors[strokes(r,c)]; } } imshow("Strokes", coloredStrokes); waitKey(); return 0; }
