From ea830704b82d49ff57dc386d4a1bf9376a82cadb Mon Sep 17 00:00:00 2001
From: gr5 <gr@gr5.org>
Date: Sat, 8 Nov 2025 18:57:48 -0500
Subject: [PATCH 1/7] Wrote a "better" gaussian blur function but it's barely
 better. And it's 200ms slower.  Will do more experiments.

---
 gaussianblur.cpp | 154 +++++++++++++++++++++++++++++++++++++++++++++++
 gaussianblur.h   |  10 +++
 2 files changed, 164 insertions(+)
 create mode 100644 gaussianblur.cpp
 create mode 100644 gaussianblur.h

diff --git a/gaussianblur.cpp b/gaussianblur.cpp
new file mode 100644
index 00000000..30bc3823
--- /dev/null
+++ b/gaussianblur.cpp
@@ -0,0 +1,154 @@
+// written by gr@gr5.org
+// gaussian blur but pays attention to borders - specifically the outer radius
+
+#include <opencv2/opencv.hpp>
+
+/*
+    cx = wf->m_outside.m_center.rx();
+    cy = wf->m_outside.m_center.ry();
+    double rad = wf->m_outside.m_radius-2;
+    */
+
+
+void CropGaussianBlur(cv::Mat in, cv::Mat out, int kernelSize, int centerx, int centery, int radius)
+{
+    // it's okay if in and out are the same materix as we do calculations in a single row buffer
+
+    // make sure kernelSize is positive and odd
+    if (kernelSize < 1) kernelSize=1;
+    if (kernelSize %2 == 0) kernelSize++;
+
+    int rad2 = (radius)*(radius);
+    int kernel_radius = (kernelSize-1)/2;
+
+    cv::Mat kernel = cv::getGaussianKernel(kernelSize,0);
+
+    // we will be using this 1 dimensional kernal matrix to do a 1 dimensional gaussian that blurrs horizontal pixels 
+    // together but ignores vertical neighbors.  Then we will start over and do the vertical blurring.  Gaussian blur lets
+    // you do this (separate vertical and horizontal blurring steps) which saves a LOT of time even with 3x3 blurring kernels
+
+     int size = in.cols;
+    double * onerow = new double[size];
+    bool * mask = new bool[size];
+    bool some_points_masked=false;
+
+    for( int y=0; y<in.rows; y++){
+        int y2 = (y-centery)*(y-centery);
+        for(int x=0; x<size; x++){
+            // first make the mask - true means inside the circle
+            int x2 = (x-centerx)*(x-centerx);
+            mask[x] = ((x2+y2) < rad2);
+        }
+        for( int x=0; x<size; x++){
+            if (mask[x] == false) {
+                // outside radius
+                onerow[x] = in.at<double>(y,x);
+                continue;
+            }
+
+            // check to see if all the points needed to create the gaussian blur for this pixel are available
+            some_points_masked=false;
+            int i = x-kernel_radius;
+            if (i<0 || i>= size || mask[i]==false)
+                some_points_masked=true;
+            else {
+                i = x+kernel_radius;
+                if (i<0 || i>= size || mask[i]==false)
+                    some_points_masked=true;
+            }
+
+            if (some_points_masked) {
+                double kernel_sum=0;
+                double sum=0;
+                int kernel_index=0;
+                for(int i = x-kernel_radius; i<=x+kernel_radius; i++) {
+                    if (i>=0 && i< size && mask[i]) {
+                        sum+= kernel.at<double>(kernel_index)*in.at<double>(y,i);
+                        kernel_sum += kernel.at<double>(kernel_index);
+                    }
+                    kernel_index++;
+                }
+                onerow[x] = sum/kernel_sum;
+
+            } else {
+                // all points available
+                double avg=0;
+                int kernel_index=0;
+                for(int i = x-kernel_radius; i<=x+kernel_radius; i++)
+                    avg+= kernel.at<double>(kernel_index++)*in.at<double>(y,i);
+                onerow[x] = avg;
+            }
+        }
+
+        // copy points to out matrix
+
+        for( int x=0; x<size; x++)
+            out.at<double>(y,x) = onerow[x];
+    }
+    delete[] onerow;
+    delete[] mask;
+
+    // now we do the vertical gaussian
+
+    size = in.rows;
+    double * onecol = new double[size];
+    mask = new bool[size];
+
+    for( int x=0; x<out.cols; x++){
+        int x2 = (x-centerx)*(x-centerx);
+        for( int y=0; y<size; y++){
+            // first make the mask - true means inside the circle
+            int y2 = (y-centery)*(y-centery);
+            mask[y] = ((x2+y2) < rad2);
+        }
+        for( int y=0; y<size; y++){
+            if (mask[y] == false) {
+                // outside radius
+                onecol[y] = out.at<double>(y,x);
+                continue;
+            }
+
+            // check to see if all the points needed to create the gaussian blur for this pixel are available
+            some_points_masked=false;
+            int i = y-kernel_radius;
+            if (i<0 || i>= size || mask[i]==false)
+                some_points_masked=true;
+            else {
+                i = y+kernel_radius;
+                if (i<0 || i>= size || mask[i]==false)
+                    some_points_masked=true;
+            }
+
+            if (some_points_masked) {
+                double kernel_sum=0;
+                double sum=0;
+                int kernel_index=0;
+                for(int i = y-kernel_radius; i<=y+kernel_radius; i++) {
+                    if (i>=0 && i< size && mask[i]) {
+                        sum+= kernel.at<double>(kernel_index) * out.at<double>(i,x);
+                        kernel_sum += kernel.at<double>(kernel_index);
+                    }
+                    kernel_index++;
+                }
+                onecol[y] = sum/kernel_sum;
+
+            } else {
+                // all points available
+                double avg=0;
+                int kernel_index=0;
+                for(int i = y-kernel_radius; i<=y+kernel_radius; i++)
+                    avg+= kernel.at<double>(kernel_index++) * out.at<double>(i,x);
+                onecol[y] = avg;
+            }
+        }
+
+        // copy points to out matrix
+
+        for(int y=0; y<size; y++)
+            out.at<double>(y,x) = onecol[y];
+    }
+
+    delete[] onecol;
+    delete[] mask;
+
+}
diff --git a/gaussianblur.h b/gaussianblur.h
new file mode 100644
index 00000000..cf0ecc1a
--- /dev/null
+++ b/gaussianblur.h
@@ -0,0 +1,10 @@
+#include <opencv2/opencv.hpp>
+
+/*
+    cx = wf->m_outside.m_center.rx();
+    cy = wf->m_outside.m_center.ry();
+    double rad = wf->m_outside.m_radius-2;
+    */
+
+
+void CropGaussianBlur(cv::Mat in, cv::Mat out, int kernelSize, int centerx, int centery, int radius);

From d2049d0caab5d6e7c462688a94bec5917004fef8 Mon Sep 17 00:00:00 2001
From: gr5 <gr@gr5.org>
Date: Thu, 13 Nov 2025 14:54:52 -0500
Subject: [PATCH 2/7] fixes ronchi issues with gaussianBlur.  gaussianBlur
 completely rewritten to handle outer and inner outlines with no turned edge.
 Now it works better for ronchi than zernike smoothing (in my opinion).

Also commented out the unused function, processSmoothing()
---
 gaussianblur.cpp    | 154 ------------------------
 gaussianblur.h      |  10 --
 simulationsview.cpp |   8 +-
 surfacemanager.cpp  |  23 ++--
 surfacemanager.h    |   2 +-
 utils.cpp           | 286 ++++++++++++++++++++++++++++++++++++++++++++
 utils.h             |   3 +
 7 files changed, 312 insertions(+), 174 deletions(-)
 delete mode 100644 gaussianblur.cpp
 delete mode 100644 gaussianblur.h

diff --git a/gaussianblur.cpp b/gaussianblur.cpp
deleted file mode 100644
index 30bc3823..00000000
--- a/gaussianblur.cpp
+++ /dev/null
@@ -1,154 +0,0 @@
-// written by gr@gr5.org
-// gaussian blur but pays attention to borders - specifically the outer radius
-
-#include <opencv2/opencv.hpp>
-
-/*
-    cx = wf->m_outside.m_center.rx();
-    cy = wf->m_outside.m_center.ry();
-    double rad = wf->m_outside.m_radius-2;
-    */
-
-
-void CropGaussianBlur(cv::Mat in, cv::Mat out, int kernelSize, int centerx, int centery, int radius)
-{
-    // it's okay if in and out are the same materix as we do calculations in a single row buffer
-
-    // make sure kernelSize is positive and odd
-    if (kernelSize < 1) kernelSize=1;
-    if (kernelSize %2 == 0) kernelSize++;
-
-    int rad2 = (radius)*(radius);
-    int kernel_radius = (kernelSize-1)/2;
-
-    cv::Mat kernel = cv::getGaussianKernel(kernelSize,0);
-
-    // we will be using this 1 dimensional kernal matrix to do a 1 dimensional gaussian that blurrs horizontal pixels 
-    // together but ignores vertical neighbors.  Then we will start over and do the vertical blurring.  Gaussian blur lets
-    // you do this (separate vertical and horizontal blurring steps) which saves a LOT of time even with 3x3 blurring kernels
-
-     int size = in.cols;
-    double * onerow = new double[size];
-    bool * mask = new bool[size];
-    bool some_points_masked=false;
-
-    for( int y=0; y<in.rows; y++){
-        int y2 = (y-centery)*(y-centery);
-        for(int x=0; x<size; x++){
-            // first make the mask - true means inside the circle
-            int x2 = (x-centerx)*(x-centerx);
-            mask[x] = ((x2+y2) < rad2);
-        }
-        for( int x=0; x<size; x++){
-            if (mask[x] == false) {
-                // outside radius
-                onerow[x] = in.at<double>(y,x);
-                continue;
-            }
-
-            // check to see if all the points needed to create the gaussian blur for this pixel are available
-            some_points_masked=false;
-            int i = x-kernel_radius;
-            if (i<0 || i>= size || mask[i]==false)
-                some_points_masked=true;
-            else {
-                i = x+kernel_radius;
-                if (i<0 || i>= size || mask[i]==false)
-                    some_points_masked=true;
-            }
-
-            if (some_points_masked) {
-                double kernel_sum=0;
-                double sum=0;
-                int kernel_index=0;
-                for(int i = x-kernel_radius; i<=x+kernel_radius; i++) {
-                    if (i>=0 && i< size && mask[i]) {
-                        sum+= kernel.at<double>(kernel_index)*in.at<double>(y,i);
-                        kernel_sum += kernel.at<double>(kernel_index);
-                    }
-                    kernel_index++;
-                }
-                onerow[x] = sum/kernel_sum;
-
-            } else {
-                // all points available
-                double avg=0;
-                int kernel_index=0;
-                for(int i = x-kernel_radius; i<=x+kernel_radius; i++)
-                    avg+= kernel.at<double>(kernel_index++)*in.at<double>(y,i);
-                onerow[x] = avg;
-            }
-        }
-
-        // copy points to out matrix
-
-        for( int x=0; x<size; x++)
-            out.at<double>(y,x) = onerow[x];
-    }
-    delete[] onerow;
-    delete[] mask;
-
-    // now we do the vertical gaussian
-
-    size = in.rows;
-    double * onecol = new double[size];
-    mask = new bool[size];
-
-    for( int x=0; x<out.cols; x++){
-        int x2 = (x-centerx)*(x-centerx);
-        for( int y=0; y<size; y++){
-            // first make the mask - true means inside the circle
-            int y2 = (y-centery)*(y-centery);
-            mask[y] = ((x2+y2) < rad2);
-        }
-        for( int y=0; y<size; y++){
-            if (mask[y] == false) {
-                // outside radius
-                onecol[y] = out.at<double>(y,x);
-                continue;
-            }
-
-            // check to see if all the points needed to create the gaussian blur for this pixel are available
-            some_points_masked=false;
-            int i = y-kernel_radius;
-            if (i<0 || i>= size || mask[i]==false)
-                some_points_masked=true;
-            else {
-                i = y+kernel_radius;
-                if (i<0 || i>= size || mask[i]==false)
-                    some_points_masked=true;
-            }
-
-            if (some_points_masked) {
-                double kernel_sum=0;
-                double sum=0;
-                int kernel_index=0;
-                for(int i = y-kernel_radius; i<=y+kernel_radius; i++) {
-                    if (i>=0 && i< size && mask[i]) {
-                        sum+= kernel.at<double>(kernel_index) * out.at<double>(i,x);
-                        kernel_sum += kernel.at<double>(kernel_index);
-                    }
-                    kernel_index++;
-                }
-                onecol[y] = sum/kernel_sum;
-
-            } else {
-                // all points available
-                double avg=0;
-                int kernel_index=0;
-                for(int i = y-kernel_radius; i<=y+kernel_radius; i++)
-                    avg+= kernel.at<double>(kernel_index++) * out.at<double>(i,x);
-                onecol[y] = avg;
-            }
-        }
-
-        // copy points to out matrix
-
-        for(int y=0; y<size; y++)
-            out.at<double>(y,x) = onecol[y];
-    }
-
-    delete[] onecol;
-    delete[] mask;
-
-}
diff --git a/gaussianblur.h b/gaussianblur.h
deleted file mode 100644
index cf0ecc1a..00000000
--- a/gaussianblur.h
+++ /dev/null
@@ -1,10 +0,0 @@
-#include <opencv2/opencv.hpp>
-
-/*
-    cx = wf->m_outside.m_center.rx();
-    cy = wf->m_outside.m_center.ry();
-    double rad = wf->m_outside.m_radius-2;
-    */
-
-
-void CropGaussianBlur(cv::Mat in, cv::Mat out, int kernelSize, int centerx, int centery, int radius);
diff --git a/simulationsview.cpp b/simulationsview.cpp
index 6228ba40..03971c8f 100644
--- a/simulationsview.cpp
+++ b/simulationsview.cpp
@@ -35,6 +35,8 @@
 #include <QTextDocument>
 #include <QtMath>
 #include <opencv2/core/core_c.h>
+#include "utils.h"
+
 double M2PI = M_PI * 2.;
 SimulationsView *SimulationsView::m_Instance = 0;
 class arcSecScaleDraw: public QwtScaleDraw
@@ -205,13 +207,15 @@ cv::Mat SimulationsView::nulledSurface(double defocus){
     zernEnables[3] = false;
     cv::Mat nulled_surface = zp.null_unwrapped( *(m_Instance->m_wf), newZerns, zernEnables);
     zernEnables[3] = saved_defocus_enable;
+
     if (GB_enabled){
         double gbValue = settings.value("GBValue", 21).toInt();
-        int blurRad = .01 * gbValue * md->diameter;
+        int blurRad = .01 * gbValue * m_wf->m_outside.m_radius * 2;
         blurRad &= 0xfffffffe;
         ++blurRad;
-        cv::GaussianBlur( nulled_surface, nulled_surface , cv::Size( blurRad, blurRad ),0,0);
+        CropGaussianBlur(nulled_surface.clone(), nulled_surface, blurRad, m_wf->m_outside, m_wf->m_inside);
     }
+
     nulled_surface  *= M2PI * md->lambda/outputLambda;
     return nulled_surface;
 }
diff --git a/surfacemanager.cpp b/surfacemanager.cpp
index 0bd7dd3c..a1028ff1 100644
--- a/surfacemanager.cpp
+++ b/surfacemanager.cpp
@@ -74,7 +74,8 @@
 #include "oglrendered.h"
 #include "ui_oglrendered.h"
 #include "astigpolargraph.h"
-
+#include <chrono>
+#include "utils.h"
 
 cv::Mat theMask;
 cv::Mat deb;
@@ -359,7 +360,7 @@ void SurfaceManager::generateSurfacefromWavefront(wavefront * wf){
                 gaussianRad &= 0xfffffffe;
 
                 ++gaussianRad;
-                    cv::GaussianBlur( wf->nulledData.clone(), wf->workData,
+                cv::GaussianBlur( wf->nulledData.clone(), wf->workData,
                                       cv::Size( gaussianRad, gaussianRad ),0,0,cv::BORDER_REFLECT);
             }
             else {
@@ -429,14 +430,22 @@ void SurfaceManager::generateSurfacefromWavefront(wavefront * wf){
 
 
     if (m_GB_enabled){
-            expandBorder(wf);
             // compute blur radius
             int gaussianRad = 2 * wf->m_outside.m_radius * m_gbValue * .01;
 
             gaussianRad &= 0xfffffffe;
             ++gaussianRad;
-            cv::GaussianBlur( wf->nulledData.clone(), wf->workData,
-                              cv::Size( gaussianRad, gaussianRad ),0,0,cv::BORDER_REFLECT);
+
+            //std::chrono::milliseconds ms = std::chrono::duration_cast< std::chrono::milliseconds >(
+            //    std::chrono::system_clock::now().time_since_epoch());
+
+            CropGaussianBlur(wf->nulledData.clone(), wf->workData, gaussianRad, wf->m_outside, wf->m_inside);
+
+            //std::chrono::milliseconds ms2 = std::chrono::duration_cast< std::chrono::milliseconds >(
+            //    std::chrono::system_clock::now().time_since_epoch());
+
+            //int duration = ms2.count() - ms.count();
+            //spdlog::get("logger")->trace("guassian blur time in ms: {}", duration);
     }
 
     wf->nulledData.release();
@@ -1366,7 +1375,7 @@ void SurfaceManager::deleteCurrent(){
 
     emit currentNdxChanged(m_currentNdx);
 }
-
+/*
 void SurfaceManager::processSmoothing(){
     if (m_wavefronts.size() == 0)
         return;
@@ -1387,7 +1396,7 @@ void SurfaceManager::processSmoothing(){
     }
 
     sendSurface(wf);
-}
+}*/
 
 void SurfaceManager::next(){
     if (m_wavefronts.size() == 0)
diff --git a/surfacemanager.h b/surfacemanager.h
index 5f6f1052..de8108f3 100644
--- a/surfacemanager.h
+++ b/surfacemanager.h
@@ -68,7 +68,7 @@ class SurfaceManager : public QObject
     void previous();
     void next();
     void deleteCurrent();
-    void processSmoothing();
+    //void processSmoothing();
     void saveAllWaveFrontStats();
     void SaveWavefronts(bool saveNulled);
     void writeWavefront(const QString &fname, wavefront *wf, bool saveNulled);
diff --git a/utils.cpp b/utils.cpp
index 0a2fa633..35dacb7c 100644
--- a/utils.cpp
+++ b/utils.cpp
@@ -1,4 +1,6 @@
 #include <QDebug>
+#include <opencv2/opencv.hpp>
+#include "Circleoutline.h"
 
 #if defined(_WIN32)
 #include <Windows.h>
@@ -30,3 +32,287 @@ long showmem(const QString & /*title*/){
 #else
 int showmem(const QString&){return 1000;}
 #endif
+
+
+
+/*
+void debugshow(cv::Mat m) {
+    // PLEASE DON'T DELETE
+    // outputs do debug stream several key areas of a matrix - mainly helps figure out where the edge of the data is
+    int height = m.rows;
+    int width = m.cols;
+
+    int y=height/2;
+    int count_zeros=0;
+    int data_index=0;
+    double data[6];
+    for(int x=0;x<width;x++) {
+        double val = m.at<double>(y,x);
+        if (val == 0.0) {
+            count_zeros++;
+            continue;
+        }
+        data[data_index++]=val;
+        if (data_index>=5)
+            break;
+    }
+
+    spdlog::get("logger")->debug("middle zeros: {} data {} {} {} {} {}", count_zeros, data[0], data[1], data[2], data[3], data[4]  );
+
+    y=height/3;
+    count_zeros=0;
+    data_index=0;
+    for(int x=0;x<width;x++) {
+        double val = m.at<double>(y,x);
+        if (val == 0.0) {
+            count_zeros++;
+            continue;
+        }
+        data[data_index++]=val;
+        if (data_index>=5)
+            break;
+    }
+
+    y--;
+    count_zeros=0;
+    data_index=0;
+    for(int x=0;x<width;x++) {
+        double val = m.at<double>(y,x);
+        if (val == 0.0) {
+            count_zeros++;
+            continue;
+        }
+        data[data_index++]=val;
+        if (data_index>=5)
+            break;
+    }
+
+    spdlog::get("logger")->debug("third0 zeros: {} data {} {} {} {} {}", count_zeros, data[0], data[1], data[2], data[3], data[4]  );
+
+    y--;
+    count_zeros=0;
+    data_index=0;
+    for(int x=0;x<width;x++) {
+        double val = m.at<double>(y,x);
+        if (val == 0.0) {
+            count_zeros++;
+            continue;
+        }
+        data[data_index++]=val;
+        if (data_index>=5)
+            break;
+    }
+
+    spdlog::get("logger")->debug("third1 zeros: {} data {} {} {} {} {}", count_zeros, data[0], data[1], data[2], data[3], data[4]  );
+
+    //
+    // right side
+    //
+
+    y = height/3;
+    count_zeros=0;
+    for(int x=width;x>0;x--) {
+        double val = m.at<double>(y,x);
+        if (val == 0.0) {
+            count_zeros++;
+            continue;
+        }
+        else break;
+    }
+
+    spdlog::get("logger")->debug("third0 right side zeros: {} ", count_zeros);
+
+    y--;
+    count_zeros=0;
+    for(int x=width;x>0;x--) {
+        double val = m.at<double>(y,x);
+        if (val == 0.0) {
+            count_zeros++;
+            continue;
+        }
+        else break;
+    }
+
+    spdlog::get("logger")->debug("third1 right side zeros: {} ", count_zeros);
+
+    y--;
+    count_zeros=0;
+    for(int x=width;x>0;x--) {
+        double val = m.at<double>(y,x);
+        if (val == 0.0) {
+            count_zeros++;
+            continue;
+        }
+        else break;
+    }
+
+    spdlog::get("logger")->debug("third2 right side zeros: {} ", count_zeros);
+
+    //
+    // inside circle
+    //
+
+    y=height/2;
+    int x=width/2;
+    int min,max;
+    while(x>0 && m.at<double>(y,x) == 0.0)x--;
+    min = x;
+    x=width/2;
+    while(x<width && m.at<double>(y,x) == 0.0)x++;
+    max = x;
+    spdlog::get("logger")->debug("middle left/right side zeros: {} {} ", min, max);
+
+    y=height/3;
+    x=width/2;
+    while(x>0 && m.at<double>(y,x) == 0.0)x--;
+    min = x;
+    x=width/2;
+    while(x<width && m.at<double>(y,x) == 0.0)x++;
+    max = x;
+    spdlog::get("logger")->debug("third0 left/right side zeros: {} {} ", min, max);
+
+    y--;
+    x=width/2;
+    while(x>0 && m.at<double>(y,x) == 0.0)x--;
+    min = x;
+    x=width/2;
+    while(x<width && m.at<double>(y,x) == 0.0)x++;
+    max = x;
+    spdlog::get("logger")->debug("third1 left/right side zeros: {} {} ", min, max);
+
+    y--;
+    x=width/2;
+    while(x>0 && m.at<double>(y,x) == 0.0)x--;
+    min = x;
+    x=width/2;
+    while(x<width && m.at<double>(y,x) == 0.0)x++;
+    max = x;
+    spdlog::get("logger")->debug("third2 left/right side zeros: {} {} ", min, max);
+
+
+}
+
+*/
+
+
+void CropGaussianBlur(cv::Mat in, cv::Mat out, int kernelSize, const CircleOutline &outside, const CircleOutline &center)
+{
+    // Warning - this function is extremely sensitive to getting the mask off by even just a few pixels out of all 360k or so.  The mask
+    // must be perfect to every last pixel.  If the outer mask zeros some valid pixels we get a turned edge.  Very visible in the 3d view.
+    // if it's off the other way you will see some pixels at level 0 outside of the wavefront.  Either error messes up foucault and
+    // ronchi simulations
+
+    // Right now we shrink outer outline by 2 pixels and enlarge inner outline by 1 pixel. As far as I can tell, this is an arbitrary
+    // hack.  These values should at a minimum be in some constant somewhere or better, part of the outline class or maybe we should
+    // get rid of these alltogether?  For now, know that if you mess with it one one place in DFTFringe, you have to mess with it
+    // in other places as well.
+
+
+    // first we make a mask - with ones where the data is and zeros where there is no data.  So if 0 is black it will be a white donut
+    // using outside and center outlines (often center outline has radius zero so then a disk instead of a donut)
+    // we also zero out the input data which probably is already zeroed but we can't take that chance particularly if there were wavefronts averaged
+    // because the average feature fills in the center and outside the outer outline.
+
+    int height = in.rows;
+    int width = in.cols;
+    cv::Mat mask = cv::Mat::ones(height,width,CV_64F);
+
+    //
+    // outer circle of mask
+    //
+
+    double cx = outside.m_center.x();
+    double cy = outside.m_center.y();
+    double radius2 = (outside.m_radius-2)*(outside.m_radius-2);
+    for (int x=0; x<width; x++) {
+        double dx = (cx-x)*(cx-x);
+        for (int y=0; y<width; y++) {
+            double dy = (cy-y)*(cy-y);
+            if ((dx+dy) >= radius2) {
+                mask.at<double>(y,x) = 1e-5; // can't use zero or we get divide by zero issues later (in cells we don't care about)
+                in.at<double>(y,x) = 0;
+            }
+        }
+    }
+
+    //
+    // inner circle of mask
+    //
+
+    cx = center.m_center.x();
+    cy = center.m_center.y();
+    radius2 = (center.m_radius+1)*(center.m_radius+1);
+    if (center.m_radius > 0) {
+        for (int x=0; x<width; x++) {
+            double dx = (cx-x)*(cx-x);
+            for (int y=0; y<width; y++) {
+                double dy = (cy-y)*(cy-y);
+                if ((dx+dy) < radius2) {
+                    mask.at<double>(y,x) = 1e-5; // can't use zero or we get divide by zero issues later (in cells we don't care about)
+                    in.at<double>(y,x) = 0;
+                }
+            }
+        }
+    }
+
+
+
+    // these next 3 lines of code are brilliant and from julien asking chatgpt but I get it completely - if someone wants an explanation
+    // I'll have to draw some pictures...  Although I try to explain with words farther below
+
+    cv::GaussianBlur(in,out,cv::Size(kernelSize,kernelSize),0);
+    cv::GaussianBlur(mask,mask,cv::Size(kernelSize,kernelSize),0);
+    out = out / mask; // corrects for the zeroed pixels that are outside the circle (were set to zero and were in the outer mask area or center area for annular wavefronts)
+
+    //
+    // now zero out points outside of mask
+    //
+
+    //
+    // outer circle of mask
+    //
+
+
+    cx = outside.m_center.x();
+    cy = outside.m_center.y();
+    radius2 = (outside.m_radius-2)*(outside.m_radius-2);
+
+    for (int x=0; x<width; x++) {
+        double dx2 = (cx-x)*(cx-x);
+        for (int y=0; y<width; y++) {
+            double dy2 = (cy-y)*(cy-y);
+            if ((dx2+dy2) >= radius2) {
+                out.at<double>(y,x) = 0;
+            }
+        }
+    }
+
+    //
+    // inner circle of mask
+    //
+
+    cx = center.m_center.x();
+    cy = center.m_center.y();
+    radius2 = (center.m_radius+1)*(center.m_radius+1);
+    if (center.m_radius > 0) {
+        for (int x=0; x<width; x++) {
+            double dx = (cx-x)*(cx-x);
+            for (int y=0; y<width; y++) {
+                double dy = (cy-y)*(cy-y);
+                if ((dx+dy) < radius2) {
+                    out.at<double>(y,x) = 0;
+                }
+            }
+        }
+    }
+
+    // here's an attempt at an explanation of the 3 "magic" lines of code.  When you are getting the gaussian blur of a pixel near the edge you are adding
+    // up a bunch of points convolved (multiplied) by the kernel shape.  But the pixels outside the circle are all zero and because the kernel is
+    // normalized you normally don't have to divide by the sum of the kernel values. but you don't want to use the kernel values for points outside
+    // the circle (and you sort of haven't used them because they are all multiplied by zero) but in this case you need to divide by the kernel
+    // values that fall inside the circle.  So we make a mask and blur it the same amount and this gives us the proper weighting.  Points far from
+    // the edge will all be set to 1 so it won't affect the weighting.  Points near the edge will have a smaller weighting value.  aka normalizing
+    // value aka sum of kernel points.
+
+
+}
diff --git a/utils.h b/utils.h
index bb6c2555..131e34a1 100644
--- a/utils.h
+++ b/utils.h
@@ -2,4 +2,7 @@
 #define UTILS_H
 int showmem(const QString &title = "");
 extern double outputLambda; //TODO nothing to do here. But there are many other occurences
+#include <opencv2/opencv.hpp>
+void CropGaussianBlur(cv::Mat in, cv::Mat out, int kernelSize, const CircleOutline &outside, const CircleOutline &center);
+//void debugshow(cv::Mat m);
 #endif // UTILS_H

From 251ca52f4fddb3070f382b5eb62775739e041b9f Mon Sep 17 00:00:00 2001
From: gr5 <gr5@users.noreply.github.com>
Date: Fri, 14 Nov 2025 09:52:48 -0500
Subject: [PATCH 3/7] Update utils.h

Co-authored-by: Julien Staub <atsju2@yahoo.fr>
---
 utils.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/utils.h b/utils.h
index 131e34a1..3272eb76 100644
--- a/utils.h
+++ b/utils.h
@@ -3,6 +3,6 @@
 int showmem(const QString &title = "");
 extern double outputLambda; //TODO nothing to do here. But there are many other occurences
 #include <opencv2/opencv.hpp>
-void CropGaussianBlur(cv::Mat in, cv::Mat out, int kernelSize, const CircleOutline &outside, const CircleOutline &center);
+void CropGaussianBlur(const cv::Mat &in, cv::Mat &out, int kernelSize, const CircleOutline &outside, const CircleOutline &center);
 //void debugshow(cv::Mat m);
 #endif // UTILS_H

From 59b5238b405b6f5783bf6e0ef37583d6326ab0cf Mon Sep 17 00:00:00 2001
From: gr5 <gr5@users.noreply.github.com>
Date: Fri, 14 Nov 2025 09:53:03 -0500
Subject: [PATCH 4/7] Update surfacemanager.cpp

---
 surfacemanager.cpp | 1 -
 1 file changed, 1 deletion(-)

diff --git a/surfacemanager.cpp b/surfacemanager.cpp
index 06cff821..d3528be1 100644
--- a/surfacemanager.cpp
+++ b/surfacemanager.cpp
@@ -74,7 +74,6 @@
 #include "oglrendered.h"
 #include "ui_oglrendered.h"
 #include "astigpolargraph.h"
-#include <chrono>
 #include "utils.h"
 
 cv::Mat theMask;

From 43c3cc8231cced14e35ee76866ce7eda4a12be80 Mon Sep 17 00:00:00 2001
From: gr5 <gr5@users.noreply.github.com>
Date: Fri, 14 Nov 2025 09:53:14 -0500
Subject: [PATCH 5/7] Update utils.cpp

Co-authored-by: Julien Staub <atsju2@yahoo.fr>
---
 utils.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/utils.cpp b/utils.cpp
index 35dacb7c..7e41875c 100644
--- a/utils.cpp
+++ b/utils.cpp
@@ -213,8 +213,8 @@ void CropGaussianBlur(cv::Mat in, cv::Mat out, int kernelSize, const CircleOutli
     // we also zero out the input data which probably is already zeroed but we can't take that chance particularly if there were wavefronts averaged
     // because the average feature fills in the center and outside the outer outline.
 
-    int height = in.rows;
-    int width = in.cols;
+    const int height = in.rows;
+    const int width = in.cols;
     cv::Mat mask = cv::Mat::ones(height,width,CV_64F);
 
     //

From cf9343219614c180388c9bfb8411ab7c23ec709d Mon Sep 17 00:00:00 2001
From: gr5 <gr5@users.noreply.github.com>
Date: Fri, 14 Nov 2025 09:53:27 -0500
Subject: [PATCH 6/7] Update utils.cpp

Co-authored-by: Julien Staub <atsju2@yahoo.fr>
---
 utils.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/utils.cpp b/utils.cpp
index 7e41875c..38cc72b3 100644
--- a/utils.cpp
+++ b/utils.cpp
@@ -36,7 +36,7 @@ int showmem(const QString&){return 1000;}
 
 
 /*
-void debugshow(cv::Mat m) {
+void debugshow(const cv::Mat &m) {
     // PLEASE DON'T DELETE
     // outputs do debug stream several key areas of a matrix - mainly helps figure out where the edge of the data is
     int height = m.rows;

From 243dc9ad8ba7adfcfbb617e631c41b0719224752 Mon Sep 17 00:00:00 2001
From: gr5 <gr@gr5.org>
Date: Fri, 14 Nov 2025 13:16:15 -0500
Subject: [PATCH 7/7] Fixed comment about how gaussianBlur works.  Made
 CropGaussianBlur() first parameter be a const and pass by reference and
 second parameter is also passed by reference.

---
 simulationsview.cpp |  2 +-
 surfacemanager.cpp  |  2 +-
 utils.cpp           | 26 ++++++++++++++++----------
 3 files changed, 18 insertions(+), 12 deletions(-)

diff --git a/simulationsview.cpp b/simulationsview.cpp
index 7c975b96..12fff464 100644
--- a/simulationsview.cpp
+++ b/simulationsview.cpp
@@ -218,7 +218,7 @@ cv::Mat SimulationsView::nulledSurface(double defocus){
         int blurRad = .01 * gbValue * m_wf->m_outside.m_radius * 2;
         blurRad &= 0xfffffffe;
         ++blurRad;
-        CropGaussianBlur(nulled_surface.clone(), nulled_surface, blurRad, m_wf->m_outside, m_wf->m_inside);
+        CropGaussianBlur(nulled_surface, nulled_surface, blurRad, m_wf->m_outside, m_wf->m_inside);
     }
 
     nulled_surface  *= M2PI * md->lambda/outputLambda;
diff --git a/surfacemanager.cpp b/surfacemanager.cpp
index d3528be1..d594de2c 100644
--- a/surfacemanager.cpp
+++ b/surfacemanager.cpp
@@ -440,7 +440,7 @@ void SurfaceManager::generateSurfacefromWavefront(wavefront * wf){
             //std::chrono::milliseconds ms = std::chrono::duration_cast< std::chrono::milliseconds >(
             //    std::chrono::system_clock::now().time_since_epoch());
 
-            CropGaussianBlur(wf->nulledData.clone(), wf->workData, gaussianRad, wf->m_outside, wf->m_inside);
+            CropGaussianBlur(wf->nulledData, wf->workData, gaussianRad, wf->m_outside, wf->m_inside);
 
             //std::chrono::milliseconds ms2 = std::chrono::duration_cast< std::chrono::milliseconds >(
             //    std::chrono::system_clock::now().time_since_epoch());
diff --git a/utils.cpp b/utils.cpp
index 38cc72b3..11b8dc30 100644
--- a/utils.cpp
+++ b/utils.cpp
@@ -195,7 +195,7 @@ void debugshow(const cv::Mat &m) {
 */
 
 
-void CropGaussianBlur(cv::Mat in, cv::Mat out, int kernelSize, const CircleOutline &outside, const CircleOutline &center)
+void CropGaussianBlur(const cv::Mat &in_, cv::Mat &out, int kernelSize, const CircleOutline &outside, const CircleOutline &center)
 {
     // Warning - this function is extremely sensitive to getting the mask off by even just a few pixels out of all 360k or so.  The mask
     // must be perfect to every last pixel.  If the outer mask zeros some valid pixels we get a turned edge.  Very visible in the 3d view.
@@ -213,6 +213,7 @@ void CropGaussianBlur(cv::Mat in, cv::Mat out, int kernelSize, const CircleOutli
     // we also zero out the input data which probably is already zeroed but we can't take that chance particularly if there were wavefronts averaged
     // because the average feature fills in the center and outside the outer outline.
 
+    cv::Mat in = in_.clone();
     const int height = in.rows;
     const int width = in.cols;
     cv::Mat mask = cv::Mat::ones(height,width,CV_64F);
@@ -257,8 +258,20 @@ void CropGaussianBlur(cv::Mat in, cv::Mat out, int kernelSize, const CircleOutli
 
 
 
-    // these next 3 lines of code are brilliant and from julien asking chatgpt but I get it completely - if someone wants an explanation
-    // I'll have to draw some pictures...  Although I try to explain with words farther below
+    // these next 3 lines of code are brilliant and from chatgpt but I get it completely.  But it's very hard to explain without pictures.
+    // here's an attempt: When you are getting the gaussian blur of a pixel near the edge (near the edge of the mirror or other mask edge) you are adding
+    // up a bunch of points convolved (multiplied) by the kernel shape.  But the pixels outside the mask should be ignored. and because the kernel is
+    // normalized (adds up to 1) you normally don't have to divide by the sum of the kernel values. but you don't want to use the kernel values for points outside
+    // the mirror edge so we set those to zero in the "in" matrix in code above (and you sort of haven't used them because they are all multiplied by zero) but
+    // in this case you need to divide by the kernel values that fall inside the mirror outline.  So we make a mask (with 1s for valid data) and blur it
+    // the same amount and this gives us the proper weighting. In the blurred mask, points far from the edge will all be set to 1 so it won't affect the weighting.
+    // Points just inside the edge of the mirror will have a smaller weighting value. These points near the edge will have a value of the sum of kernel points.  Just
+    // what we need to divide the points by.
+    // This is "mask" in the division below.  To be clear, "mask" will have values less than 1 (greater than 0) near the edge of the mask and is the perfect amount
+    // to divide into "out" to get the properly weighted final result.
+    //
+    // Note the 3 lines of code below can have a very complicated mask with many masked regions and it all still just works and in fact cropGaussianBlur() supports
+    // the inner outline as well as part of the mask.
 
     cv::GaussianBlur(in,out,cv::Size(kernelSize,kernelSize),0);
     cv::GaussianBlur(mask,mask,cv::Size(kernelSize,kernelSize),0);
@@ -306,13 +319,6 @@ void CropGaussianBlur(cv::Mat in, cv::Mat out, int kernelSize, const CircleOutli
         }
     }
 
-    // here's an attempt at an explanation of the 3 "magic" lines of code.  When you are getting the gaussian blur of a pixel near the edge you are adding
-    // up a bunch of points convolved (multiplied) by the kernel shape.  But the pixels outside the circle are all zero and because the kernel is
-    // normalized you normally don't have to divide by the sum of the kernel values. but you don't want to use the kernel values for points outside
-    // the circle (and you sort of haven't used them because they are all multiplied by zero) but in this case you need to divide by the kernel
-    // values that fall inside the circle.  So we make a mask and blur it the same amount and this gives us the proper weighting.  Points far from
-    // the edge will all be set to 1 so it won't affect the weighting.  Points near the edge will have a smaller weighting value.  aka normalizing
-    // value aka sum of kernel points.
 
 
 }