diff --git a/simulationsview.cpp b/simulationsview.cpp index 31247701..12fff464 100644 --- a/simulationsview.cpp +++ b/simulationsview.cpp @@ -35,6 +35,7 @@ #include #include #include +#include "utils.h" #include #include #include @@ -211,17 +212,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 = 2 *(m_Instance->m_wf)->m_outside.m_radius * gbValue * .01; + if (GB_enabled){ + double gbValue = settings.value("GBValue", 21).toInt(); + int blurRad = .01 * gbValue * m_wf->m_outside.m_radius * 2; blurRad &= 0xfffffffe; ++blurRad; - - - qDebug() << "Blurr" << blurRad; - cv::GaussianBlur( nulled_surface, nulled_surface , cv::Size( blurRad, blurRad ),0,0); + CropGaussianBlur(nulled_surface, 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 4686f7a2..d594de2c 100644 --- a/surfacemanager.cpp +++ b/surfacemanager.cpp @@ -74,7 +74,7 @@ #include "oglrendered.h" #include "ui_oglrendered.h" #include "astigpolargraph.h" - +#include "utils.h" cv::Mat theMask; cv::Mat deb; @@ -359,8 +359,9 @@ void SurfaceManager::generateSurfacefromWavefront(wavefront * wf){ gaussianRad &= 0xfffffffe; ++gaussianRad; - qDebug() << "Blurr" << gaussianRad; - cv::GaussianBlur( wf->nulledData.clone(), wf->workData, + + qDebug() << "Blurr" << gaussianRad; + cv::GaussianBlur( wf->nulledData.clone(), wf->workData, cv::Size( gaussianRad, gaussianRad ),0,0,cv::BORDER_REFLECT); } else { @@ -430,15 +431,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, 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(); @@ -1368,7 +1376,7 @@ void SurfaceManager::deleteCurrent(){ emit currentNdxChanged(m_currentNdx); } - +/* void SurfaceManager::processSmoothing(){ if (m_wavefronts.size() == 0) return; @@ -1389,7 +1397,7 @@ void SurfaceManager::processSmoothing(){ } sendSurface(wf); -} +}*/ void SurfaceManager::next(){ if (m_wavefronts.size() == 0) @@ -2440,7 +2448,7 @@ textres SurfaceManager::Phase2(QList list, QList inp void SurfaceManager::computeStandAstig(define_input *wizPage, QList list){ // check for pairs QVector lookat = list.toVector(); - spdlog::get("logger")->trace("computeStandAstig()"); + spdlog::get("logger")->trace("computeStandAstig()"); while (lookat.size()){ for (int i = 0; i < lookat.size(); ++i){ double angle1 = wrapAngle(lookat[i]->angle); 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..11b8dc30 100644 --- a/utils.cpp +++ b/utils.cpp @@ -1,4 +1,6 @@ #include +#include +#include "Circleoutline.h" #if defined(_WIN32) #include @@ -30,3 +32,293 @@ long showmem(const QString & /*title*/){ #else int showmem(const QString&){return 1000;} #endif + + + +/* +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; + int width = m.cols; + + int y=height/2; + int count_zeros=0; + int data_index=0; + double data[6]; + for(int x=0;x(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(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(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(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(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(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(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(y,x) == 0.0)x--; + min = x; + x=width/2; + while(x(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(y,x) == 0.0)x--; + min = x; + x=width/2; + while(x(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(y,x) == 0.0)x--; + min = x; + x=width/2; + while(x(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(y,x) == 0.0)x--; + min = x; + x=width/2; + while(x(y,x) == 0.0)x++; + max = x; + spdlog::get("logger")->debug("third2 left/right side zeros: {} {} ", min, max); + + +} + +*/ + + +void CropGaussianBlur(const cv::Mat &in_, cv::Mat &out, int kernelSize, const CircleOutline &outside, const CircleOutline ¢er) +{ + // 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. + + 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); + + // + // 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= radius2) { + mask.at(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(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(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(y,x) = 0; + } + } + } + } + + + + // 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); + 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= radius2) { + out.at(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(y,x) = 0; + } + } + } + } + + + +} diff --git a/utils.h b/utils.h index bb6c2555..3272eb76 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 +void CropGaussianBlur(const cv::Mat &in, cv::Mat &out, int kernelSize, const CircleOutline &outside, const CircleOutline ¢er); +//void debugshow(cv::Mat m); #endif // UTILS_H