aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorFranklin Wei <franklin@rockbox.org>2020-02-12 23:05:50 -0500
committerFranklin Wei <franklin@rockbox.org>2020-02-12 23:05:50 -0500
commitac6c10a72be75bf6b1fa6260bd473cb5067c35bf (patch)
tree0fe0f04b4f0d62a0cdf45314cf9995317ac2c813
parent09acfb1ff934499e1a40aec6da496836014fe594 (diff)
downloadsloreg-ac6c10a72be75bf6b1fa6260bd473cb5067c35bf.zip
sloreg-ac6c10a72be75bf6b1fa6260bd473cb5067c35bf.tar.gz
sloreg-ac6c10a72be75bf6b1fa6260bd473cb5067c35bf.tar.bz2
sloreg-ac6c10a72be75bf6b1fa6260bd473cb5067c35bf.tar.xz
Refactor
-rw-r--r--gabor.cpp431
1 files changed, 226 insertions, 205 deletions
diff --git a/gabor.cpp b/gabor.cpp
index 3c9aa4d..24822b9 100644
--- a/gabor.cpp
+++ b/gabor.cpp
@@ -14,6 +14,7 @@
#include <iostream>
#include <opencv2/opencv.hpp>
#include <string>
+#include <vector>
using namespace cv;
using namespace std;
@@ -21,7 +22,9 @@ using namespace std;
// uncomment for headless mode -- for benchmarking
//#define imshow(a, b)
-// parameters were lifted off the internet
+// Perform a maximal Gabor filtering of in with nfilts oriented
+// filters at ksize x ksize kernel size. Parameters were lifted off
+// the internet and need fine-tuning.
Mat maximalGaborFilter(Mat in, int nfilts = 10,
int ksize = 9, // kernel size
double sig = 3, // variance -- how wide the filter is
@@ -45,8 +48,6 @@ Mat maximalGaborFilter(Mat in, int nfilts = 10,
Size sz = filtered[0].size();
- cout << "Result size: " << sz << endl;
-
Mat result = Mat(sz, CV_32F, Scalar(0));
float *rows[nfilts];
@@ -58,14 +59,29 @@ Mat maximalGaborFilter(Mat in, int nfilts = 10,
float v = -1;
for(int i = 0; i < nfilts; i++)
v = std::max(v, rows[i][x]);
- //cout << y << ", " << x << endl;
outrow[x] = v;
}
}
return result;
}
-// hist out must have size 256
+// Normalized sum of maximal gabor filtered image with varied kernel
+// size.
+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, nfilts, i);
+ result += sf * m;
+ }
+ return result;
+}
+
+// Get a list of intensity frequencies for all intensities
+// 0-255. *hist should point to a 256-element int array.
int getHist(Mat img, int *hist) {
assert(img.type() == CV_8U);
memset(hist, 0, sizeof(int) * 256);
@@ -76,6 +92,12 @@ int getHist(Mat img, int *hist) {
return ymax;
}
+void dumpHist(int *hist) {
+ for(int i = 0; i < 256; i+=1)
+ cout << setw(4) << i << " " << hist[i] << endl;
+}
+
+// Plot a histogram. highlight is a column to be highlighted in red.
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++) {
@@ -84,7 +106,7 @@ Mat plotHist(int *hist, int ymax, int highlight) {
return img;
}
-// return a 512-element integer array giving the "prominence" of each
+// 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
@@ -92,6 +114,8 @@ Mat plotHist(int *hist, int ymax, int highlight) {
//
// Sorting the array by pairs of ints will give the peaks by
// prominence.
+
+// TODO: refactor with a struct
void getPeakProminences(const int *hist, int *out) {
for(int i = 0; i < 256; i++) {
*out++ = i;
@@ -125,7 +149,8 @@ int compare_pair(const void *a, const void *b) {
#define ABS(x) (((x) < 0) ? (-(x)) : (x))
// sorted_peak_proms is a 512-element int array of (peak, prominence)
-// tuples. This function looks for numerically adjacent
+// tuples. This function looks for numerically adjacent peaks of equal
+// prominence and filters them. Might not work.
void dedupe_adjacent(int *sorted_peak_proms) {
int out[512];
int out_idx = 0;
@@ -144,34 +169,6 @@ void dedupe_adjacent(int *sorted_peak_proms) {
// go past maxLum, which is the threshold returned by Otsu's method,
// which is consistently an overestimate).
int getThreshold(const int *hist) {
-#if 0
- // "cross section" the histogram
- for(int y = 20; y < 256; y++) {
- int colorchanges = 0;
- int first_end = -1;
- bool last = hist[0] >= y;
- for(int x = 1; x < 256; x++) {
- bool curr = hist[x] >= y;
- if(curr != last)
- {
- cout << "Color changes at " << x << ", " << y << endl;
- colorchanges++;
- if(first_end < 0 && last && !curr)
- first_end = x;
- }
- last = curr;
- }
-
- cout << "y" << y << " has " << colorchanges << "changes" << endl;
-
- int sections = (colorchanges + 1) / 2;
- if(sections == 2) {
- return first_end;
- }
- }
- return -1;
-#else
- // use this version
int peak_proms[512];
getPeakProminences(hist, peak_proms);
qsort(peak_proms, 256, sizeof(int) * 2, compare_pair);
@@ -197,7 +194,7 @@ int getThreshold(const int *hist) {
idx1 = idx2;
idx2 = tmp;
}
-// assert(idx2 < idx1);
+ //assert(idx2 < idx1);
// find minimum
int min_idx = idx2;
@@ -206,22 +203,6 @@ int getThreshold(const int *hist) {
min_idx = i;
return min_idx;
-#endif
-}
-
-// normalized sum of maximal gabor filtered image with varied kernel
-// size
-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, nfilts, i);
- result += sf * m;
- }
- return result;
}
bool inBounds(int x, int y, Size sz) {
@@ -230,7 +211,11 @@ bool inBounds(int x, int y, Size sz) {
#define ARRAYLEN(x) (sizeof(x) / sizeof(x[0]))
-void floodFill(int x, int y, Mat map, Mat visited, unsigned char target, unsigned char replace) {
+// Iterative flood fill to clean up an image mask.
+void floodFill(int x, int y,
+ Mat map, Mat visited,
+ unsigned char target, unsigned char replace)
+{
queue<pair<int, int> > q;
q.push(pair<int, int>(x, y));
@@ -263,22 +248,12 @@ void floodFill(int x, int y, Mat map, Mat visited, unsigned char target, unsigne
}
}
-Mat removeSatRegion(Mat mask, Mat orig) {
- static bool firstRun = true;
-
- static Rect roi;
- if(firstRun)
- roi = selectROI("select saturation region", orig), firstRun = false;
-
- rectangle(mask, roi, Scalar(0), FILLED);
- return mask;
-}
-
-// flood fill in from the border to isolate the center region
+// Flood fill in from the border to isolate the center region.
// 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, 255, 0);
floodFill(x, h - 1, map, visited, 255, 0);
@@ -289,34 +264,43 @@ void borderFlood(Mat map, Mat visited) {
}
}
+// Clean up our raw mask, produced by thresholding the MOG filtered
+// images, by removing outside "islands" of white and "holes" of
+// black.
Mat cleanMask(Mat initial) {
- // now filter the thresholded mask to extract the region
- // we care about
int h = initial.size().height, w = initial.size().width;
+ // white = ROI
Mat mask = Mat(h, w, CV_8U, Scalar(255));
// flood fill from the four boundaries to get the region
// we want
borderFlood(initial, mask);
- //return mask;
-
- // now to get rid of disconnected blobs we flood fill from
- // the center (for now we assume the center point is in
- // the largest region -- this can be changed with a "real"
- // algorithm to find the largest connected region)
+ // Now to get rid of disconnected blobs, we flood fill from the
+ // center (this assumes the center point is in the largest region
+ // -- but this can be changed in a "real" algorithm to find the
+ // largest connected region)
Mat mask2 = Mat(h, w, CV_8U, Scalar(0));
floodFill(w / 2, h / 2, mask, mask2, 0, 255);
return mask2;
}
-void dumpHist(int *hist) {
- for(int i = 0; i < 256; i+=1)
- cout << setw(4) << i << " " << hist[i] << endl;
+// Prompt the user to select the saturation image in orig. Zeroes the
+// corresponding region in mask.
+Mat removeSatRegion(Mat mask, Mat orig) {
+ static bool firstRun = true;
+
+ static Rect roi;
+ if(firstRun)
+ roi = selectROI("select saturation region", orig), firstRun = false;
+
+ rectangle(mask, roi, Scalar(0), FILLED);
+ return mask;
}
+// Generate a red-green overlay.
Mat redGreenOverlay(Mat rPlane, Mat gPlane) {
assert(rPlane.size() == gPlane.size());
Mat planes[3] = { Mat::zeros(rPlane.size(), CV_8U),
@@ -327,168 +311,205 @@ Mat redGreenOverlay(Mat rPlane, Mat gPlane) {
return out;
}
-int main()
+// Perform automatic segmentation. Returns the masked initial frame,
+// along with auxilliary outputs.
+struct processed_frame {
+ Mat masked_frame,
+ masked_frame_filtered;
+ Mat mask;
+};
+
+processed_frame preprocessFrame(Mat frame)
{
// filtering parameters
const int nfilts = 10;
- Mat reference, reference_filtered;
- Mat mosaic;
+ resize(frame, frame, Size(500, 500));
- // alignment parameters
- //bool
+ // grayscale conversion
+ Mat frame_gray;
+ cvtColor(frame, frame_gray, COLOR_BGR2GRAY);
+ //GaussianBlur(frame_gray, frame_gray, Size(15, 15), 0, 0);
- bool homography = false;
- int ecc_iters = 100;
- double ecc_eps = 1e-5;
- int warp_mode = MOTION_TRANSLATION;
+ imshow("orig", frame_gray);
-#define NFRAMES 22
+ Mat gray_inverted = Scalar::all(255) - frame_gray;
- //while(1)
- {
- // 22 frames
- for(int i = 1; i <= NFRAMES; i++)
- {
- char buf[128];
- snprintf(buf, sizeof(buf), "SLO Data for registration/SLO001/SLO_subject001_frame%d.png", i);
- cout << "Image: " << buf << endl;
+ // convert to float
+ gray_inverted.convertTo(gray_inverted, CV_32F, 1.0, 0);
+ normalize(gray_inverted, gray_inverted, 0, 1, NORM_MINMAX);
- Mat img_c3;
- img_c3 = imread(buf);
+ // Gabor filtering
+ Mat filter_result = summedGaborFilter(gray_inverted, nfilts, 5, 13, 2);
+ filter_result.convertTo(filter_result, CV_8UC1, -255, 255); // cast back to char, uninvert while we're at it
+ imshow("filtered", filter_result);
- resize(img_c3, img_c3, Size(500, 500));
+ // calculate histogram
+ int hist[256];
+ int ymax = getHist(filter_result, hist);
+ //dumpHist(hist);
- // grayscale conversion
- Mat img;
- cvtColor(img_c3, img, COLOR_BGR2GRAY);
- //GaussianBlur(img, img, Size(15, 15), 0, 0);
+ // blur filtered filter_result
+ Mat blurred;
+ GaussianBlur(filter_result, blurred, Size(5, 5), 0, 0);
- imshow("orig", img);
+ // threshold out background
+ Mat threshed;
- Mat inverted = Scalar::all(255) - img;
+ // Otsu's method doesn't work -- too high of a threshold
+ //double t = threshold(blurred, threshed, 0, 255, THRESH_BINARY | THRESH_OTSU);
+ //imshow("thresholded", threshed);
+ //cout << "Otsu threshold: " << t << endl;
- // convert to float
- inverted.convertTo(inverted, CV_32F, 1.0, 0);
- normalize(inverted, inverted, 0, 1, NORM_MINMAX);
+ int t = getThreshold(hist);
- // Gabor filtering
- Mat result = summedGaborFilter(inverted, nfilts, 5, 13, 2);
- result.convertTo(result, CV_8UC1, -255, 255); // uninvert while we're at it
- imshow("filtered", result);
+ cout << "Our threshold: " << t << endl;
- // calculate histogram
- int hist[256];
- int ymax = getHist(result, hist);
- //dumpHist(hist);
+ // get initial mask by thresholding the filtered image
+ threshold(blurred, threshed, t, 255, THRESH_BINARY);
+ imshow("thresholded", threshed);
- // blur filtered result
- Mat blurred;
- GaussianBlur(result, blurred, Size(5, 5), 0, 0);
+ // clean up mask -- remove islands and holes
+ Mat mask = cleanMask(threshed);
- // threshold out background
- Mat threshed;
+ // remove saturation region from mask
+ mask = removeSatRegion(mask, frame_gray);
- // Otsu's method doesn't work -- too high of a threshold
- //double t = threshold(blurred, threshed, 0, 255, THRESH_BINARY | THRESH_OTSU);
- //imshow("thresholded", threshed);
- //cout << "Otsu threshold: " << t << endl;
+ imshow("mask", mask);
- int t = getThreshold(hist);
+ // if you want a magenta background
+ Mat masked_result = Mat(mask.size(), CV_8UC3, Scalar(0xff, 0, 0xff));
+ //Mat masked_result;
+ frame_gray.copyTo(masked_result, mask);
- cout << "Our threshold: " << t << endl;
+ imshow("masked result", masked_result);
- // get initial mask by thresholding the filtered image
- threshold(blurred, threshed, t, 255, THRESH_BINARY);
- imshow("thresholded", threshed);
+ Mat masked_filter_result = Mat(mask.size(), CV_8UC3, Scalar(0xff, 0, 0xff));
+ //Mat masked_filter_result;
+ filter_result.copyTo(masked_filter_result, mask);
- // clean up mask -- remove islands and holes
- Mat mask = cleanMask(threshed);
+ imshow("extracted & filtered", masked_filter_result);
- // remove saturation region from mask
- mask = removeSatRegion(mask, img);
+ // plot histogram
+ Mat histplot = plotHist(hist, ymax, t);
+ imshow("hist", histplot);
- imshow("mask", mask);
+ processed_frame ret = { masked_result, masked_filter_result, mask };
- // if you want a magenta background
- Mat masked = Mat(mask.size(), CV_8UC3, Scalar(0xff, 0, 0xff));
- //Mat masked;
- img.copyTo(masked, mask);
+ return ret;
+}
- imshow("extracted", masked);
+Mat alignFrames(Mat reference_filtered, processed_frame moving)
+{
+ int ecc_iters = 100;
+ double ecc_eps = 1e-5;
- Mat masked_filter_result = Mat(mask.size(), CV_8UC3, Scalar(0xff, 0, 0xff));
- //Mat masked_filter_result;
- result.copyTo(masked_filter_result, mask);
- imshow("extracted & filtered", masked_filter_result);
+ // Translation works best. Higher DoF modes (affine, homography)
+ // are slow and introduce errors.
+ int warp_mode = MOTION_TRANSLATION;
- // plot histogram
- Mat histplot = plotHist(hist, ymax, t);
- imshow("hist", histplot);
+ Mat warp_matrix;
+ if(warp_mode == MOTION_HOMOGRAPHY)
+ warp_matrix = Mat::eye(3, 3, CV_32F);
+ else
+ warp_matrix = Mat::eye(2, 3, CV_32F);
- // try alignment
- if(reference.empty())
- {
- reference = masked;
- reference_filtered = masked_filter_result;
+ double coeff;
- mosaic = Mat(reference.rows, reference.cols, CV_32F);
- mosaic = reference / NFRAMES;
- }
- else
- {
- Mat warp_matrix;
- if(warp_mode == MOTION_HOMOGRAPHY)
- warp_matrix = Mat::eye(3, 3, CV_32F);
- else
- warp_matrix = Mat::eye(2, 3, CV_32F);
-
- double coeff;
-
- try
- {
- coeff = findTransformECC(/*reference,*/ reference_filtered,
- /*masked,*/ masked_filter_result,
- warp_matrix,
- warp_mode,
- TermCriteria (TermCriteria::COUNT+TermCriteria::EPS,
- ecc_iters, ecc_eps),
- mask);
-
- Mat warped_image = Mat(reference.rows, reference.cols, CV_32FC1);
-
- if (warp_mode != MOTION_HOMOGRAPHY)
- warpAffine(masked, warped_image, warp_matrix, warped_image.size(),
- INTER_LINEAR + WARP_INVERSE_MAP);
- else
- warpPerspective(masked, warped_image, warp_matrix, warped_image.size(),
- INTER_LINEAR + WARP_INVERSE_MAP);
-
- imshow("aligned", warped_image);
-
- cout << "ECC coef: " << coeff << endl;
- cout << "Sizes: " << warped_image.size() << ", " << reference.size() << endl;
-
- // create red/green overlay
- Mat overlay = redGreenOverlay(warped_image, reference);
- imshow("overlay", overlay);
-
- // assemble mosaic
- mosaic += warped_image / NFRAMES;
- }
- catch(cv::Exception e) {
- cerr << "Alignment failed: " << e.what() << endl;
- }
- }
+ try
+ {
+ coeff = findTransformECC(reference_filtered,
+ moving.masked_frame_filtered,
+ warp_matrix,
+ warp_mode,
+ TermCriteria (TermCriteria::COUNT+TermCriteria::EPS,
+ ecc_iters, ecc_eps),
+ moving.mask);
+
+ Mat warped_image = Mat(reference_filtered.rows, reference_filtered.cols, CV_32FC1);
+
+ if (warp_mode != MOTION_HOMOGRAPHY)
+ warpAffine(moving.masked_frame, warped_image, warp_matrix, warped_image.size(),
+ INTER_LINEAR + WARP_INVERSE_MAP);
+ else
+ warpPerspective(moving.masked_frame, warped_image, warp_matrix, warped_image.size(),
+ INTER_LINEAR + WARP_INVERSE_MAP);
- imshow("mosaic", mosaic);
+ imshow("aligned", warped_image);
- waitKey(0);
- }
- reference = mosaic;
- mosaic = Mat(reference.rows, reference.cols, CV_32F);
+ cout << "ECC coef: " << coeff << endl;
- cout << "Registration complete." << endl;
- waitKey(0);
+ return warped_image;
}
+ catch(cv::Exception e) {
+ cerr << "Alignment failed: " << e.what() << endl;
+ }
+
+ return Mat();
+}
+
+Mat averageStack(vector<Mat> aligned_frames)
+{
+ double sf = 1. / aligned_frames.size();
+
+ Mat result = Mat(aligned_frames[0].size(), CV_32F, Scalar(0));
+ for(Mat i : aligned_frames)
+ {
+ Mat scaled;
+ i.convertTo(scaled, CV_32F, sf / 255.0);
+ result += scaled;
+ }
+
+ return result;
+}
+
+Mat produceMosaic(vector<Mat> frames)
+{
+ vector<processed_frame> preproc(frames.size());
+
+ for(int i = 0; i < frames.size(); i++) {
+ preproc[i] = preprocessFrame(frames[i]);
+ }
+
+ // align all to first frame and mosaic
+ vector<Mat> aligned(frames.size()), overlays(frames.size());
+
+ for(int i = 0; i < frames.size(); i++) {
+ aligned[i] = alignFrames(preproc[0].masked_frame, preproc[i]);
+ overlays[i] = redGreenOverlay(aligned[i], preproc[0].masked_frame);
+ imshow("overlay", overlays[i]);
+ }
+
+ // now average
+ Mat mosaic = averageStack(aligned);
+
+ imshow("mosaic", mosaic);
+
+ waitKey(0);
+
+ return mosaic;
+}
+
+#define NFRAMES 22
+
+vector<Mat> loadData()
+{
+ vector<Mat> stack(NFRAMES);
+
+ for(int i = 1; i <= NFRAMES; i++)
+ {
+ char buf[128];
+ snprintf(buf, sizeof(buf), "SLO Data for registration/SLO001/SLO_subject001_frame%d.png", i);
+
+ stack[i - 1] = imread(buf);
+ }
+
+ return stack;
+}
+
+int main()
+{
+ vector<Mat> stack = loadData();
+
+ produceMosaic(stack);
}