aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--Makefile2
-rw-r--r--gabor.cpp193
2 files changed, 137 insertions, 58 deletions
diff --git a/Makefile b/Makefile
index 72a8bf2..fde6f50 100644
--- a/Makefile
+++ b/Makefile
@@ -1,7 +1,7 @@
all: gabor test Makefile
CXXFLAGS+=$(shell pkg-config --cflags --libs opencv4)
-CXXFLAGS+=-O3
+CXXFLAGS+=-O1
CXXFLAGS+=-g
clean:
rm -f gabor test
diff --git a/gabor.cpp b/gabor.cpp
index 24822b9..d07f5ec 100644
--- a/gabor.cpp
+++ b/gabor.cpp
@@ -1,8 +1,8 @@
/***************************************************************************
*
- * Real-time SLO image registration
+ * Real-time retinal image registration
*
- * Copyright (C) 2019 Franklin Wei
+ * Copyright (C) 2019-2020 Franklin Wei
*
****************************************************************************/
@@ -14,8 +14,15 @@
#include <iostream>
#include <opencv2/opencv.hpp>
#include <string>
+#include <sys/stat.h>
#include <vector>
+//#define NULL_ALIGNMENT
+
+#ifdef NULL_ALIGNMENT
+#warning Null alignment - testing only
+#endif
+
using namespace cv;
using namespace std;
@@ -25,9 +32,9 @@ using namespace std;
// 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,
+Mat maximalGaborFilter(Mat in, int nfilts = 12,
int ksize = 9, // kernel size
- double sig = 3, // variance -- how wide the filter is
+ double sig = 2.5, // variance -- how wide the filter is
double lm = 8, // wavelength ?
double gm = 0.02, // aspect ratio (smaller = longer vessels), 1 : square, >1 : bad
double ps = 0) { // phase shift -- don't change?
@@ -165,10 +172,21 @@ void dedupe_adjacent(int *sorted_peak_proms) {
memcpy(sorted_peak_proms, out, sizeof(out));
}
-// 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) {
+#define HISTOGRAM_NOISE_FLOOR .05
+
+void zeroNoise(int *hist, int thres) {
+ for(int i = 0; i < 256; i++)
+ if(hist[i] < thres) hist[i] = 0;
+}
+
+// Get the threshold to filter out the peripheral noise.
+// Does a preprocessing step followed by a peak-finding.
+int getThreshold(const int *hist_orig, int ymax) {
+ int hist[256];
+ memcpy(hist, hist_orig, sizeof(hist));
+
+ zeroNoise(hist, ymax * HISTOGRAM_NOISE_FLOOR);
+
int peak_proms[512];
getPeakProminences(hist, peak_proms);
qsort(peak_proms, 256, sizeof(int) * 2, compare_pair);
@@ -301,22 +319,30 @@ Mat removeSatRegion(Mat mask, Mat orig) {
}
// Generate a red-green overlay.
-Mat redGreenOverlay(Mat rPlane, Mat gPlane) {
+Mat redGreenOverlay(Mat rPlane, Mat gPlane, double rWeight = 1, double gWeight = 1) {
assert(rPlane.size() == gPlane.size());
Mat planes[3] = { Mat::zeros(rPlane.size(), CV_8U),
- gPlane * .6,
- rPlane };
+ gPlane * gWeight,
+ rPlane * rWeight};
Mat out;
merge(planes, 3, out);
return out;
}
+// A Gabor-filtered frame and its unfiltered counterpart
+struct filtered_frame {
+ Mat filtered, orig;
+};
+
+// Holds an image and
// Perform automatic segmentation. Returns the masked initial frame,
// along with auxilliary outputs.
struct processed_frame {
+ Mat orig_frame;
Mat masked_frame,
masked_frame_filtered;
Mat mask;
+ int histogram[256];
};
processed_frame preprocessFrame(Mat frame)
@@ -331,7 +357,7 @@ processed_frame preprocessFrame(Mat frame)
cvtColor(frame, frame_gray, COLOR_BGR2GRAY);
//GaussianBlur(frame_gray, frame_gray, Size(15, 15), 0, 0);
- imshow("orig", frame_gray);
+ //imshow("orig", frame_gray);
Mat gray_inverted = Scalar::all(255) - frame_gray;
@@ -342,7 +368,7 @@ processed_frame preprocessFrame(Mat frame)
// 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);
+ //imshow("filtered", filter_result);
// calculate histogram
int hist[256];
@@ -361,45 +387,48 @@ processed_frame preprocessFrame(Mat frame)
//imshow("thresholded", threshed);
//cout << "Otsu threshold: " << t << endl;
- int t = getThreshold(hist);
+ int t = getThreshold(hist, ymax);
cout << "Our threshold: " << t << endl;
// get initial mask by thresholding the filtered image
threshold(blurred, threshed, t, 255, THRESH_BINARY);
- imshow("thresholded", threshed);
+ //imshow("thresholded", threshed);
// clean up mask -- remove islands and holes
Mat mask = cleanMask(threshed);
// remove saturation region from mask
- mask = removeSatRegion(mask, frame_gray);
+ // disabled
+ //mask = removeSatRegion(mask, frame_gray);
- imshow("mask", mask);
+ //imshow("mask", mask);
// 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);
- imshow("masked result", masked_result);
+ //imshow("masked result", masked_result);
Mat masked_filter_result = Mat(mask.size(), CV_8UC3, Scalar(0xff, 0, 0xff));
//Mat masked_filter_result;
filter_result.copyTo(masked_filter_result, mask);
- imshow("extracted & filtered", masked_filter_result);
+ //imshow("extracted & filtered", masked_filter_result);
// plot histogram
Mat histplot = plotHist(hist, ymax, t);
- imshow("hist", histplot);
+ //imshow("hist", histplot);
- processed_frame ret = { masked_result, masked_filter_result, mask };
+ processed_frame ret = { frame, masked_result, masked_filter_result, mask };
+
+ memcpy(ret.histogram, hist, sizeof(ret.histogram));
return ret;
}
-Mat alignFrames(Mat reference_filtered, processed_frame moving)
+Mat alignECC(Mat reference_filtered, processed_frame moving)
{
int ecc_iters = 100;
double ecc_eps = 1e-5;
@@ -416,36 +445,57 @@ Mat alignFrames(Mat reference_filtered, processed_frame moving)
double coeff;
- 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);
+ coeff = findTransformECC(reference_filtered,
+ moving.masked_frame_filtered,
+ warp_matrix,
+ warp_mode,
+ TermCriteria (TermCriteria::COUNT+TermCriteria::EPS,
+ ecc_iters, ecc_eps),
+ moving.mask);
- imshow("aligned", warped_image);
+ Mat globally_aligned = Mat(reference_filtered.rows, reference_filtered.cols, CV_32FC1);
- cout << "ECC coef: " << coeff << endl;
+ if (warp_mode != MOTION_HOMOGRAPHY)
+ warpAffine(moving.masked_frame, globally_aligned, warp_matrix, globally_aligned.size(),
+ INTER_LINEAR + WARP_INVERSE_MAP);
+ else
+ warpPerspective(moving.masked_frame, globally_aligned, warp_matrix, globally_aligned.size(),
+ INTER_LINEAR + WARP_INVERSE_MAP);
- return warped_image;
- }
- catch(cv::Exception e) {
- cerr << "Alignment failed: " << e.what() << endl;
- }
+ imshow("aligned", globally_aligned);
+
+ cout << "ECC coef: " << coeff << endl;
+
+ // do local registration
+
+
+ return globally_aligned;
+}
+
+Mat alignFrames(Mat reference_filtered, processed_frame moving)
+{
+#ifdef NULL_ALIGNMENT
+ return moving.masked_frame;
+#endif
+
+ Mat globally_aligned = alignECC(reference_filtered, moving);
- return Mat();
+ // do local registration
+
+ // subdivide into nxn
+ const int n = 5;
+
+ int region_w = globally_aligned.width / n,
+ region_h = globally_aligned.height / n;
+
+ for(int i = 0; i < n; i++) {
+ for(int j = 0; j < n; j++) {
+ Mat sub = globally_aligned(Rect(region_w * i, region_h * j,
+ region_w, region_h)); // pointer semantics - careful!
+
+
+ }
+ }
}
Mat averageStack(vector<Mat> aligned_frames)
@@ -465,6 +515,8 @@ Mat averageStack(vector<Mat> aligned_frames)
Mat produceMosaic(vector<Mat> frames)
{
+ mkdir("output", 0755);
+
vector<processed_frame> preproc(frames.size());
for(int i = 0; i < frames.size(); i++) {
@@ -475,9 +527,25 @@ Mat produceMosaic(vector<Mat> frames)
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]);
+ try {
+ aligned[i] = alignFrames(preproc[0].masked_frame, preproc[i]);
+ } catch (cv::Exception e) {
+ aligned[i] = preproc[i].masked_frame; // null alignment if failed
+ cerr << "Alignment failed: " << e.what() << endl;
+ }
+
overlays[i] = redGreenOverlay(aligned[i], preproc[0].masked_frame);
+
imshow("overlay", overlays[i]);
+ imshow("filtered", preproc[i].masked_frame_filtered);
+
+ char out_overlay[128], out_filtered[128];
+ snprintf(out_overlay, sizeof out_overlay, "output/overlay_%d.png", i + 1);
+ snprintf(out_filtered, sizeof out_filtered, "output/filtered_%d.png", i + 1);
+ imwrite(out_overlay, overlays[i]);
+ imwrite(out_filtered, preproc[i].masked_frame_filtered);
+
+ waitKey(0);
}
// now average
@@ -492,24 +560,35 @@ Mat produceMosaic(vector<Mat> frames)
#define NFRAMES 22
-vector<Mat> loadData()
+vector<Mat> loadData(int argc, char *argv[])
{
- vector<Mat> stack(NFRAMES);
+ if(argc <= 1)
+ throw "Usage: gabor [IMAGE...]";
+ vector<Mat> stack(argc - 1);
- for(int i = 1; i <= NFRAMES; i++)
+ for(int i = 1; i < argc; i++)
{
- char buf[128];
- snprintf(buf, sizeof(buf), "SLO Data for registration/SLO001/SLO_subject001_frame%d.png", i);
+ stack[i - 1] = imread(argv[i]);
- stack[i - 1] = imread(buf);
+ if(!stack[i - 1].data)
+ throw "File not found";
}
return stack;
}
-int main()
+int main(int argc, char *argv[])
{
- vector<Mat> stack = loadData();
+ vector<Mat> stack;
+
+ try {
+ stack = loadData(argc, argv);
+ } catch(const char *e) {
+ cerr << e << endl;
+ return 1;
+ }
produceMosaic(stack);
+
+ return 0;
}