diff options
author | Franklin Wei <franklin@rockbox.org> | 2019-12-11 14:39:02 -0500 |
---|---|---|
committer | Franklin Wei <franklin@rockbox.org> | 2019-12-11 14:39:02 -0500 |
commit | 699805f2bf28a8602a2b906c71d9f9fa0757fa91 (patch) | |
tree | 13caa370ddbdffd7e217af3e7ffbd6793b8f9dce | |
parent | 98224b1be453c99fa07457848e9c6b3adce5f314 (diff) | |
download | sloreg-699805f2bf28a8602a2b906c71d9f9fa0757fa91.zip sloreg-699805f2bf28a8602a2b906c71d9f9fa0757fa91.tar.gz sloreg-699805f2bf28a8602a2b906c71d9f9fa0757fa91.tar.bz2 sloreg-699805f2bf28a8602a2b906c71d9f9fa0757fa91.tar.xz |
Add background segmentation
-rw-r--r-- | Makefile | 8 | ||||
-rw-r--r-- | gabor.cpp | 221 |
2 files changed, 212 insertions, 17 deletions
@@ -1,5 +1,7 @@ -all: gabor Makefile +all: gabor test Makefile CXXFLAGS+=$(shell pkg-config --cflags --libs opencv4) -CXXFLAGS+=-O2 -gabor: gabor.cpp +CXXFLAGS+=-O3 +CXXFLAGS+=-g +clean: + rm -f gabor test @@ -2,6 +2,7 @@ #include <cassert> #include <cmath> #include <cstdio> +#include <iomanip> #include <iostream> #include <opencv2/opencv.hpp> #include <string> @@ -46,67 +47,259 @@ Mat maximalGaborFilter(Mat in, int nfilts = 10, Mat result = Mat(sz, CV_32F, Scalar(0)); + float *rows[nfilts]; for(int y = 0; y < sz.height; y++) { + for(int i = 0; i < nfilts; i++) + rows[i] = (float*)filtered[i].ptr<float>(y); + float *outrow = (float*)result.ptr(y); for(int x = 0; x < sz.width; x++) { float v = -1; for(int i = 0; i < nfilts; i++) - v = std::max(v, filtered[i].at<float>(y, x)); + v = std::max(v, rows[i][x]); //cout << y << ", " << x << endl; - result.at<float>(y, x) = v; - //outrow[x] = v; + outrow[x] = v; } } return result; } +// hist out must have size 256 +int getHist(Mat img, int *hist) { + assert(img.type() == CV_8U); + memset(hist, 0, sizeof(int) * 256); + unsigned char *ptr = img.ptr<unsigned char>(0); + int ymax = 0; + for(int i = 0; i < img.size().height * img.size().width; i++) + ymax = std::max(++hist[*ptr++], ymax); + return ymax; +} + +Mat plotHist(int *hist, int ymax, int highlight) { + Mat img = Mat(Size(256, 256), CV_8UC3, Scalar(0xff,0xff,0xff)); + for(int x = 0; x < 256; x++) { + line(img, Point(x, 255), Point(x, 255 - hist[x] / (double)ymax * 255), x == highlight ? Scalar(0, 0, 0xff) : Scalar(0,0,0)); + } + return img; +} + +// return a 512-element integer array giving the "prominence" of each +// element of hist with its preceding index -- that is, the minimum +// distance in either direction that one needs to travel to find an +// element of equal or greater value (for the maximum value in the +// array, it will have prominence 256). +// +// Sorting the array by pairs of ints will give the peaks by +// prominence. +void getPeakProminences(const int *hist, int *out) { + for(int i = 0; i < 256; i++) { + *out++ = i; + int elem = hist[i]; + bool done = false; + for(int j = 1; !done && (i - j >= 0 || i + j < 256); j++) + { + if(i - j >= 0 && hist[i - j] >= elem) + *out = j, done = true; + else if(i + j < 256 && hist[i + j] >= elem) + *out = j, done = true; + } + if(!done) + *out = 999; + out++; + } +} + +int compare_pair(const void *a, const void *b) { + const int *l = (const int*)a + 1, *r = (const int*)b + 1; + if(*l < *r) return 1; + if(*l > *r) return -1; + return 0; +} + +// get the threshold to filter out the periphial blind spots (will not +// go past maxLum, which is the threshold returned by Otsu's method, +// which is consistently an overestimate). +int getThreshold(const int *hist) { + // "cross section" the histogram +#if 0 + for(int y = 0; y < 256; y++) { + int sections = 0; + + } +#endif + + int peak_proms[512]; + getPeakProminences(hist, peak_proms); + qsort(peak_proms, 256, sizeof(int) * 2, compare_pair); + int idx1 = peak_proms[0], idx2 = peak_proms[2]; + cout << "Minimum of interest lies between " << idx2 << ", " << idx1 << endl; + +// assert(idx2 < idx1); + + // find minimum + int min_idx = idx2; + for(int i = idx2 + 1; i < idx1; i++) + if(hist[i] < hist[min_idx]) + min_idx = i; + + return min_idx; +} + // normalized sum of maximal gabor filtered image with varied kernel // size -Mat summedGaborFilter(Mat in, int start, int stop, int step) { +Mat summedGaborFilter(Mat in, int nfilts, int start, int stop, int step) { assert(!(step & 1) && (start & 1) && (stop & 1)); int n = (stop - start) / step + 2; double sf = 1. / n; Mat result = Mat(in.size(), CV_32F, Scalar(0)); for(int i = start; i <= stop; i += step) { - Mat m = maximalGaborFilter(in, 10, i); + Mat m = maximalGaborFilter(in, nfilts, i); result += sf * m; } return result; } +bool inBounds(int x, int y, Size sz) { + return (0 <= x && x < sz.width) && (0 <= y && y < sz.height); +} + +#define ARRAYLEN(x) (sizeof(x) / sizeof(x[0])) + +void floodFill(int x, int y, Mat map, Mat visited) { + queue<pair<int, int> > q; + + q.push(pair<int, int>(x, y)); + + while(q.size() > 0) { + pair<int, int> p = q.front(); + q.pop(); + x = p.first; + y = p.second; + if(!inBounds(x, y, map.size())) + continue; + if(!visited.at<char>(y, x) || map.at<char>(y, x)) + continue; + visited.at<char>(y, x) = 0; + + int delts[] = { -1, 0, 1, 0, 0, -1, 0, 1 }; + for(int i = 0; i < ARRAYLEN(delts); i += 2) { + int xp = x + delts[i + 0], + yp = y + delts[i + 1]; + q.push(pair<int, int>(xp, yp)); + } + } + + /* + if(!inBounds(x, y, map.size()) || !visited.at<unsigned char>(y, x) || map.at<unsigned char>(y, x)) + return; + visited.at<unsigned char>(y, x) = 0; + floodFill(x + 1, y, map, visited); + floodFill(x - 1, y, map, visited); + floodFill(x, y + 1, map, visited); + floodFill(x, y - 1, map, visited); + */ + + //imshow("progress", visited); + //waitKey(1); +} + +// in visited, != 0 means not visited, 0 means visited +void borderFlood(Mat map, Mat visited) { + Size sz = map.size(); + int w = sz.width, h = sz.height; + for(int x = 0; x < w; x++) { + floodFill(x, 0, map, visited); + floodFill(x, h - 1, map, visited); + } + for(int y = 1; y < h - 1; y++) { + floodFill(0, y, map, visited); + floodFill(w - 1, y, map, visited); + } +} + int main() { + const int nfilts = 10; while(1) { - for(int i = 1; i < 23; i++) + // 22 frames + for(int i = 1; i <= 22; i++) { char buf[64]; snprintf(buf, sizeof(buf), "SLO Data for registration/SLO001/SLO_subject001_frame%d.png", i); - Mat img; - img = imread(buf); - cvtColor(img, img, COLOR_BGR2GRAY); + Mat img_c3; + img_c3 = imread(buf); + + resize(img_c3, img_c3, Size(512, 512)); - resize(img, img, Size(512, 512)); + // grayscale + Mat img; + cvtColor(img_c3, img, COLOR_BGR2GRAY); imshow("orig", img); - Mat img_invert = Scalar::all(255) - img; + Mat inverted = Scalar::all(255) - img; - img_invert.convertTo(img_invert, CV_8UC1, 1, 0); + // alpha and beta values boost contrast + inverted.convertTo(inverted, CV_32F, 1.0, 0); + normalize(inverted, inverted, 0, 1, NORM_MINMAX); /* for(int i = 5; i < 15; i += 2) { - Mat result = maximalGaborFilter(img_invert, 10, i); + Mat result = maximalGaborFilter(inverted, 10, i); imshow(string("filtered") + to_string(i), result); } */ - Mat result = summedGaborFilter(img_invert, 3, 13, 2); + Mat result = summedGaborFilter(inverted, nfilts, 5, 13, 2); + result.convertTo(result, CV_8UC1, -255, 255); // uninvert while we're at it imshow("filtered", result); + // calculate and display histogram + int hist[256]; + int ymax = getHist(result, hist); + for(int i = 0; i < 256; i+=1) + cout << setw(4) << i << " " << hist[i] << endl; + + Mat blurred; + GaussianBlur(result, blurred, Size(5, 5), 0, 0); + Mat threshed; +// double t = threshold(blurred, threshed, 0, 255, THRESH_BINARY | THRESH_OTSU); +// imshow("thresholded", threshed); + //cout << "Otsu threshold: " << t << endl; + + int t = getThreshold(hist); + + cout << "Our threshold: " << t << endl; + + threshold(blurred, threshed, t, 255, THRESH_BINARY); + imshow("thresholded", threshed); + + // now filter the thresholded mask to extract the region + // we care about + Mat mask = Mat(threshed.size().height, threshed.size().width, CV_8U, Scalar(0xff)); + + // flood fill from the four boundaries to get the region + // we want + borderFlood(threshed, mask); + + imshow("mask", mask); + + //Mat masked = Mat(mask.size(), CV_8UC3, Scalar(0xff, 0, 0xff)); + Mat masked; + img.copyTo(masked, mask); + + imshow("extracted", masked); + + Mat histplot = plotHist(hist, ymax, t); + imshow("hist", histplot); + + + //calcHist(&result, 1, &chans, Mat(), hist, 2, &histSize, &histRange, true, false); + waitKey(0); } } |