File indexing completed on 2024-05-12 03:44:33

 /*
     SPDX-FileCopyrightText: 2023 Joseph McGee <joseph.mcgee@sbcglobal.net>
 
     SPDX-License-Identifier: GPL-2.0-or-later
 */
 
 #include <QLoggingCategory>
 #include <cmath>
 #include "optimalsubexposurecalculator.h"
 #include "imagingcameradata.h"
 #include "calculatedgainsubexposuretime.h"
 #include "cameraexposureenvelope.h"
 #include "optimalexposuredetail.h"
 #include <ekos_capture_debug.h>
 
 namespace OptimalExposure
 {
 
 /*
  *  For UI presentation provide two calculation methods.
  *  1) A method for the overall calculation for
  *  exposure times over the range of gain read-noise pairs.
  *      This method would be called when changes are made to:
  *          Camera selection, Focal Ratio, Sky Quality, Filter Compensation, or Noise Tolerance
  *      This method will return an obect (CameraExposureEnvelope) containing (for ui presentation)
  *          lightPollutionElectronBaseRate, lightPollutionForOpticFocalRatio,
  *          a vector of gain sub-exposure pairs, and int values the max and min exposure time.
  *      This method is intented to be used to present the graphical display of exposures for the camera,
  *      based on a specifc SQM.  The SQM can be adjusted in the UI, and an this method will be used to refresh ui.
  *
  *  2) A method for a calculation of an exposure time (OptimalExposureDetail) at a specific gain.
  *      This method is intented to be used to present a specific exposure calculation, and
  *      resulting shot and stack noise information.  This method will interpolate for a
  *      read-noise value (between camera gain read-noise pairs), and will be called when
  *      a change is made to selected Gain on the UI.
  *      This method will return an object that conains:
  *          Precise Exposure Time in (interpolated), Shot pollution electrons, Exposure shot noise, Exposure total noise
  *          a Vector for stacks of 1 to n hours, which includes Exposure Count, Stack Total Time, Stack total noise
  *
  */
 
 /*
  *
  * More work is needed for the calculation of effect of a Filter on exposure noise and time.
  * Original effort (in Java) used an estimage of the spectrum bandwidth passed by a filter. For
  * example: on broadband filters for a one shot color camera:
  *      Optolong l-Pro apprears to pass about 165 nanometers.
  *      Optolong l-Enhance apprears to pass only about 33 nanometers.
  *
  * In this code the filter compensation has is applied as a reducer of the calulation of the
  * light pollution rate for the optic. For example, the filter compensation to the light pollution
  * would be 165 / 300 for an l-Pro filter.
  *
  * But this filter compensatoin approach is imprecise & and does not reflect reality. It might
  * be better to analyze a filter for its ability to block the distinct emmission lines found
  * in light pollution:
  *
  *      Hg  435.8, 546.1, 577, 578.1
  *      Na  598, 589.6, 615.4, 616.1
  *
  * And finally, the use of 300 nanometers (the bandwidth of the visible light spectrum) as the
  * basis for filter compensation may not be accurate. An unfiltered camera sensor may be able
  * to "see" a spectrum that is much more broad.  Is an unfiltered sensor collecting more noise
  * from beyond the range of the visible spectrum?  But other discussions on web forums regarding
  * Dr. Glovers calculations appear to use 300 as the basis for filter compensation.
  *
  */
 
 OptimalSubExposureCalculator::OptimalSubExposureCalculator(
     double aNoiseTolerance, double aSkyQuality, double aFocalRatio, double aFilterCompensation,
     ImagingCameraData &anImagingCameraData) :
     aNoiseTolerance(aNoiseTolerance),
     aSkyQuality(aSkyQuality),
     aFocalRatio(aFocalRatio),
     aFilterCompensation(aFilterCompensation),
     anImagingCameraData(anImagingCameraData) {}
 
 
 
 double OptimalSubExposureCalculator::getASkyQuality()
 {
     return aSkyQuality;
 }
 
 void OptimalSubExposureCalculator::setASkyQuality(double newASkyQuality)
 {
     aSkyQuality = newASkyQuality;
 }
 
 double OptimalSubExposureCalculator::getANoiseTolerance()
 {
     return aNoiseTolerance;
 }
 
 void OptimalSubExposureCalculator::setANoiseTolerance(double newANoiseTolerance)
 {
     aNoiseTolerance = newANoiseTolerance;
 }
 
 double OptimalSubExposureCalculator::getAFocalRatio()
 {
     return aFocalRatio;
 }
 
 void OptimalSubExposureCalculator::setAFocalRatio(double newAFocalRatio)
 {
     aFocalRatio = newAFocalRatio;
 }
 
 double OptimalSubExposureCalculator::getAFilterCompensation()
 {
     return aFilterCompensation;
 }
 
 void OptimalSubExposureCalculator::setAFilterCompensation(double newAFilterCompensation)
 {
     aFilterCompensation = newAFilterCompensation;
 }
 
 ImagingCameraData &OptimalSubExposureCalculator::getImagingCameraData()
 {
     return anImagingCameraData;
 }
 
 void OptimalSubExposureCalculator::setImagingCameraData(ImagingCameraData &newanImagingCameraData)
 {
     anImagingCameraData = newanImagingCameraData;
 }
 
 int OptimalSubExposureCalculator::getASelectedGain()
 {
     return aSelectedGain;
 }
 
 void OptimalSubExposureCalculator::setASelectedGain(int newASelectedGain)
 {
     aSelectedGain = newASelectedGain;
 }
 
 int OptimalSubExposureCalculator::getASelectedCameraReadMode() const
 {
     return aSelectedCameraReadMode;
 }
 
 void OptimalSubExposureCalculator::setASelectedCameraReadMode(int newASelectedCameraReadMode)
 {
     aSelectedCameraReadMode = newASelectedCameraReadMode;
 }
 
 
 
 QVector<CalculatedGainSubExposureTime> OptimalSubExposureCalculator::calculateGainSubExposureVector(double cFactor,
         double lightPollutionForOpticFocalRatio)
 {
     // qCInfo(KSTARS_EKOS_CAPTURE) << "Calculating gain sub-exposure vector: ";
     QVector<CalculatedGainSubExposureTime> aCalculatedGainSubExposureTimeVector;
 
     OptimalExposure::CameraGainReadMode aSelectedReadMode =
         anImagingCameraData.getCameraGainReadModeVector()[aSelectedCameraReadMode];
     // qCInfo(KSTARS_EKOS_CAPTURE) << "\t with camera read mode: "
     //  << aSelectedReadMode.getCameraGainReadModeName();
 
     QVector<OptimalExposure::CameraGainReadNoise> aCameraGainReadNoiseVector
         = aSelectedReadMode.getCameraGainReadNoiseVector();
 
     for(QVector<OptimalExposure::CameraGainReadNoise>::iterator rn = aCameraGainReadNoiseVector.begin();
             rn != aCameraGainReadNoiseVector.end(); ++rn)
     {
         int aGain = rn->getGain();
         // qCInfo(KSTARS_EKOS_CAPTURE) << "At a gain of: " << aGain;
         double aReadNoise = rn->getReadNoise();
         // qCInfo(KSTARS_EKOS_CAPTURE) << "\t camera read-noise is: " << aReadNoise;
 
         // initial sub exposure.
         double anOptimalSubExposure = 0.0;
 
 
         switch(anImagingCameraData.getSensorType())
         {
             case SENSORTYPE_MONOCHROME:
                 anOptimalSubExposure =  cFactor * pow(aReadNoise, 2) / lightPollutionForOpticFocalRatio;
                 break;
 
             case SENSORTYPE_COLOR:
                 anOptimalSubExposure = (double)3.0 * cFactor * pow(aReadNoise, 2) / lightPollutionForOpticFocalRatio;
                 break;
         }
 
         // qCInfo(KSTARS_EKOS_CAPTURE) << "anOptimalSubExposure is: " << anOptimalSubExposure;
 
         CalculatedGainSubExposureTime *aSubExposureTime = new CalculatedGainSubExposureTime(aGain, anOptimalSubExposure);
 
         aCalculatedGainSubExposureTimeVector.append(*aSubExposureTime);
 
     }
 
     return aCalculatedGainSubExposureTimeVector;
 }
 /*
     double OptimalSubExposureCalculator::calculateLightPolutionForOpticFocalRatio(double lightPollutionElectronBaseRate, double aFocalRatio) {
     // double lightPollutionRateForOptic =  lightPollutionElectronBaseRate * java.lang.Math.pow(anOptic.getFocalRatio(),-2);
         return((double) (lightPollutionElectronBaseRate * pow(aFocalRatio,-2)));
     }
 */
 
 double OptimalSubExposureCalculator::calculateLightPolutionForOpticFocalRatio(double lightPollutionElectronBaseRate,
         double aFocalRatio, double aFilterCompensation)
 {
     // double lightPollutionRateForOptic =
     //  lightPollutionElectronBaseRate
     //  * java.lang.Math.pow(anOptic.getFocalRatio(),-2);
     return(((double) (lightPollutionElectronBaseRate * pow(aFocalRatio, -2))) * aFilterCompensation);
 }
 
 double OptimalSubExposureCalculator::calculateCFactor(double aNoiseTolerance)
 {
     // cFactor = 1/( (((100+allowableNoiseIncreasePercent)/100)^2) -1)
 
     return((double) (1 / ( pow((((double)100 + (double) aNoiseTolerance) / (double)100), (double)2) - 1)));
 
 }
 
 double OptimalSubExposureCalculator::calculateLightPollutionElectronBaseRate(double skyQuality)
 {
     // Conversion curve fitted from Dr Glover data
     double base = std::stod("1.25286030612621E+27");
     double power = (double) -19.3234809465887;
 
     // New version of Dr Glover's function calculates the Electron rate at thet pixel level
     // and requires pixel size in microns and QE of sensor.
     // double baseV2 = 1009110388.7838
     // Calculation of an initial electron rate will
     // use a new fitted curve based upon a 1.0 micon pixel and 100% QE
     // curve appears to be f(sqm) = 1009110388.7838 exp (-0.921471189594521 * sqm)
     //
 
 
     return(base * pow(skyQuality, power));
 }
 
 OptimalExposure::CameraExposureEnvelope OptimalSubExposureCalculator::calculateCameraExposureEnvelope()
 {
     /*
         This method calculates the exposures for each gain setting found in the camera gain readnoise table.
         It is used to refresh the ui presentation, in prparation for a calculation of sub-exposure data with a
         specific (probably interpolated) read-noise value.
     */
 
     double lightPollutionElectronBaseRate = OptimalSubExposureCalculator::calculateLightPollutionElectronBaseRate(aSkyQuality);
     /*
     qCInfo(KSTARS_EKOS_CAPTURE) << "Calculating CameraExposureEnvelope..."
                                 << "Using Sky Quality: " << aSkyQuality
                                 << "\nConverted to lightPollutionElectronBaseRate: " << lightPollutionElectronBaseRate
                                 << "Using an Optical Focal Ratio: " << aFocalRatio;
     */
     double lightPollutionForOpticFocalRatio = calculateLightPolutionForOpticFocalRatio(lightPollutionElectronBaseRate,
             aFocalRatio, aFilterCompensation);
     /*
     qCInfo(KSTARS_EKOS_CAPTURE)
             << "\tResulting in an Light Pollution Rate for the Optic of: "
             << lightPollutionForOpticFocalRatio;
     */
     // qCDebug(KSTARS_EKOS_CAPTURE) << "Using a camera Id: " << anImagingCameraData.getCameraId();
     OptimalExposure::CameraGainReadMode aSelectedReadMode =
         anImagingCameraData.getCameraGainReadModeVector()[aSelectedCameraReadMode];
     /*
     qCInfo(KSTARS_EKOS_CAPTURE) << "\t with camera read mode: "
                                 << aSelectedReadMode.getCameraGainReadModeName()
                                 << "\tWith a read-noise table of : "
                                 << aSelectedReadMode.getCameraGainReadNoiseVector().size() << " values";
     */
     // qCInfo(KSTARS_EKOS_CAPTURE) << "Using a Noise Tolerance of: " << aNoiseTolerance;
 
     double cFactor = calculateCFactor(aNoiseTolerance);
     // qCInfo(KSTARS_EKOS_CAPTURE) << "Calculated CFactor is: " << cFactor;
 
     QVector<CalculatedGainSubExposureTime> aSubExposureTimeVector = calculateGainSubExposureVector(cFactor,
             lightPollutionForOpticFocalRatio);
 
     double exposureTimeMin = 99999999999.9;
     double exposureTimeMax = -1.0;
     // Iterate the times to capture and store the min and max exposure times.
     // (But the QCustomPlot can rescale, so this may be unnecessary
     for(QVector<OptimalExposure::CalculatedGainSubExposureTime>::iterator oe = aSubExposureTimeVector.begin();
             oe != aSubExposureTimeVector.end(); ++oe)
     {
         if(exposureTimeMin > oe->getSubExposureTime()) exposureTimeMin = oe->getSubExposureTime();
         if(exposureTimeMax < oe->getSubExposureTime()) exposureTimeMax = oe->getSubExposureTime();
     }
 
     CameraExposureEnvelope aCameraExposureEnvelope = CameraExposureEnvelope(
                 lightPollutionElectronBaseRate,
                 lightPollutionForOpticFocalRatio,
                 aSubExposureTimeVector,
                 exposureTimeMin,
                 exposureTimeMax);
 
     return(aCameraExposureEnvelope);
 }
 
 // OptimalExposure::OptimalExposureDetail OptimalSubExposureCalculator::calculateSubExposureDetail(int aSelectedGainValue)
 OptimalExposure::OptimalExposureDetail OptimalSubExposureCalculator::calculateSubExposureDetail()
 {
     /*
         This method calculates some details for a sub-exposure. It will interpolate between
         points in the camera read-noise table, to find a reasonable appoximation of the read-noise
         for a non-listed gain value.
     */
 
     // qCDebug(KSTARS_EKOS_CAPTURE) << "aSelectedGainValue: " << aSelectedGainValue;
 
     int aSelectedGain = getASelectedGain();
 
     // Look for a matching gain from the camera gain read-noise table, or identify a bracket for interpolation.
     // Interpolation may result in slight errors when the read-noise data is curve. (there's probably a better way to code this)
     double aReadNoise = 100.0;  // defaulted high just to make an error in the interpolation obvious
     bool matched = false;
     int lowerReadNoiseIndex = 0;
 
     OptimalExposure::CameraGainReadMode aSelectedReadMode =
         anImagingCameraData.getCameraGainReadModeVector()[aSelectedCameraReadMode];
 
     QVector<OptimalExposure::CameraGainReadNoise> aCameraGainReadNoiseVector =
         aSelectedReadMode.getCameraGainReadNoiseVector();
 
 
 
     for(int readNoiseIndex = 0; readNoiseIndex < aSelectedReadMode.getCameraGainReadNoiseVector().size(); readNoiseIndex++)
     {
         if(!matched)
         {
             CameraGainReadNoise aCameraGainReadNoise = aSelectedReadMode.getCameraGainReadNoiseVector()[readNoiseIndex];
             if(aCameraGainReadNoise.getGain() == aSelectedGain)
             {
                 matched = true;
                 // qCInfo(KSTARS_EKOS_CAPTURE) << " matched a camera gain ";
                 aReadNoise = aCameraGainReadNoise.getReadNoise();
                 break;
             }
             if(aCameraGainReadNoise.getGain() < aSelectedGain)
             {
                 lowerReadNoiseIndex = readNoiseIndex;
             }
         }
     }
 
     if(!matched)
     {
         int upperReadNoiseIndex = lowerReadNoiseIndex + 1;
         // qCInfo(KSTARS_EKOS_CAPTURE)
         //  << "Interpolating read noise gain bracket: " << lowerReadNoiseIndex
         //  << " and " << upperReadNoiseIndex;
 
         if (lowerReadNoiseIndex < aSelectedReadMode.getCameraGainReadNoiseVector().size() - 1)
         {
 
             CameraGainReadNoise aLowerIndexCameraReadNoise = aSelectedReadMode.getCameraGainReadNoiseVector()[lowerReadNoiseIndex];
             CameraGainReadNoise anUpperIndexCameraReadNoise = aSelectedReadMode.getCameraGainReadNoiseVector()[upperReadNoiseIndex];
 
             // interpolate a read-noise value
             double m = (anUpperIndexCameraReadNoise.getReadNoise() - aLowerIndexCameraReadNoise.getReadNoise()) /
                        (anUpperIndexCameraReadNoise.getGain() - aLowerIndexCameraReadNoise.getGain());
             aReadNoise = aLowerIndexCameraReadNoise.getReadNoise() + (m * (aSelectedGain - aLowerIndexCameraReadNoise.getGain()));
             // qCInfo(KSTARS_EKOS_CAPTURE)
             //        << "The camera read-noise for the selected gain value is interpolated to: " << aReadNoise;
         }
         else
         {
             // qCInfo(KSTARS_EKOS_CAPTURE) << " using max camera gain ";
             int maxIndex = aSelectedReadMode.getCameraGainReadNoiseVector().size() - 1;
             CameraGainReadNoise aMaxIndexCameraReadNoise = aSelectedReadMode.getCameraGainReadNoiseVector()[maxIndex];
             aReadNoise = aMaxIndexCameraReadNoise.getReadNoise();
         }
     }
     else
     {
         // qCInfo(KSTARS_EKOS_CAPTURE)
         // << "The camera read-noise for the selected gain value was matched: " << aReadNoise;
     }
 
     double anOptimalSubExposure = 0.0;
 
     double lightPollutionElectronBaseRate = OptimalSubExposureCalculator::calculateLightPollutionElectronBaseRate(aSkyQuality);
     double lightPollutionForOpticFocalRatio = calculateLightPolutionForOpticFocalRatio(lightPollutionElectronBaseRate,
             aFocalRatio, aFilterCompensation);
     double cFactor = calculateCFactor(aNoiseTolerance);
 
     switch(anImagingCameraData.getSensorType())
     {
         case SENSORTYPE_MONOCHROME:
             anOptimalSubExposure =  cFactor * pow(aReadNoise, 2) / lightPollutionForOpticFocalRatio;
             break;
 
         case SENSORTYPE_COLOR:
             anOptimalSubExposure = (double)3.0 * cFactor * pow(aReadNoise, 2) / lightPollutionForOpticFocalRatio;
             break;
     }
 
     // qCInfo(KSTARS_EKOS_CAPTURE) << "an Optimal Sub Exposure at gain: "
     //  << aSelectedGainValue << " is " << anOptimalSubExposure << " seconds";
 
     // Add the single exposure noise calcs to the response
     double anExposurePollutionElectrons = lightPollutionForOpticFocalRatio * anOptimalSubExposure;
     // qCInfo(KSTARS_EKOS_CAPTURE) << "Exposure Pollution Electrons: " << anExposurePollutionElectrons;
 
     double anExposureShotNoise = sqrt(lightPollutionForOpticFocalRatio * anOptimalSubExposure);
     // qCInfo(KSTARS_EKOS_CAPTURE) << "Exposure Shot Noise: " << anExposureShotNoise;
 
     double anExposureTotalNoise = sqrt( pow(aReadNoise, 2)
                                         + lightPollutionForOpticFocalRatio * anOptimalSubExposure);
     // qCInfo(KSTARS_EKOS_CAPTURE) << "Exposure Total Noise: " << anExposureTotalNoise;
 
     // Need to build and populate the stack estimations
     QVector<OptimalExposureStack> aStackSummary;
 
     // Loop through planned sessions for hours of stacking
     for(int sessionHours = 1; sessionHours < 41; sessionHours++)
     {
         // rounding up to ensure that exposure count is at least "sessionHours" of data.
         int anExposureCount = (int) ceil( (double)(sessionHours * 3600) /
                                           anOptimalSubExposure );
 
         int aStackTime = anExposureCount * anOptimalSubExposure;
         double aStackTotalNoise = sqrt( (double)anExposureCount * pow(aReadNoise, 2)
                                         + lightPollutionForOpticFocalRatio * (anExposureCount * anOptimalSubExposure));
 
         OptimalExposureStack *anOptimalExposureStack
             = new OptimalExposureStack(sessionHours, anExposureCount, aStackTime, aStackTotalNoise);
 
         /*
                 qCInfo(KSTARS_EKOS_CAPTURE) << sessionHours << "\t" << anExposureCount
                                             << "\t" << aStackTime << "\t" << aStackTotalNoise
                                             << "\t" << aStackTime / aStackTotalNoise;
 
                 //  << aSelectedGainValue << " is " << anOptimalSubExposure << " seconds";
         */
         aStackSummary.push_back(*anOptimalExposureStack);
 
 
     }
 
     OptimalExposureDetail anOptimalExposureDetail = OptimalExposureDetail(getASelectedGain(), anOptimalSubExposure,
             anExposurePollutionElectrons, anExposureShotNoise,  anExposureTotalNoise, aStackSummary);
 
     return(anOptimalExposureDetail);
 }
 
 
 }