aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorFranklin Wei <franklin@rockbox.org>2019-12-11 14:39:02 -0500
committerFranklin Wei <franklin@rockbox.org>2019-12-11 14:39:02 -0500
commit699805f2bf28a8602a2b906c71d9f9fa0757fa91 (patch)
tree13caa370ddbdffd7e217af3e7ffbd6793b8f9dce
parent98224b1be453c99fa07457848e9c6b3adce5f314 (diff)
downloadsloreg-699805f2bf28a8602a2b906c71d9f9fa0757fa91.zip
sloreg-699805f2bf28a8602a2b906c71d9f9fa0757fa91.tar.gz
sloreg-699805f2bf28a8602a2b906c71d9f9fa0757fa91.tar.bz2
sloreg-699805f2bf28a8602a2b906c71d9f9fa0757fa91.tar.xz
Add background segmentation
-rw-r--r--Makefile8
-rw-r--r--gabor.cpp221
2 files changed, 212 insertions, 17 deletions
diff --git a/Makefile b/Makefile
index 0198ce4..72a8bf2 100644
--- a/Makefile
+++ b/Makefile
@@ -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
diff --git a/gabor.cpp b/gabor.cpp
index f465912..e2c9029 100644
--- a/gabor.cpp
+++ b/gabor.cpp
@@ -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);
}
}