Flesh out driver program a little. Changes optimization flags.HEADmaster

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;
 }