ということでEkosのソースを読みます。まずは星像の大きさを求めるところ。星像の大きさを求めるステップは
1)星の位置を検出する
2)星像の大きさを計算する
3)安定した星像の大きさを選択する

の3段階からなります。

●まずは1)の星の検出から。
Ekosでは星の検出にいくつかの方法が用意されています。
EkosのAFの画面で、Process->Detectionを選択すると、
Gradient, Centroid, Threshold, SEP, Bahtinov の5つから選択することができます。

Screenshot from 2021-09-19 13-14-06

ソースを見てみると、この星の検出やこの次の星像の大きさをHFR(後述)の計算はすべて、FITS画像の処理として
kstars/kstars/fitsviewer/
にまとめられています。

名前 デフォルト 説明 ソース
Gradient Sobelフィルタでエッジ検出を行ったあとエッジのまとまりごとに中心を計算。中心はSobelの微分値の大きさを重みとして重心を取っている。その中で微分値の和が最も大きいエッジのまとまりについてHFRを計算する。
微分値の大きさで重心を取る意味はなんだろう???
fitsgradientdetector.cpp
Centroid   画素値のスレッショルドでエッジを求め、エッジのまとまりごとに中心を計算する。スレッショルドは最初大きな値からはじめ、指定の数のエッジのまとまりが検出されるまで下げていく。求まったエッジのまとまりの大きさが、4σより大きいものはエラーとして排除する fitscentroiddetector.cpp
Threshold   スレッショルド以上の領域について画素値で重心をもとめ、星の位置とする。指定した領域に複数の星が存在した場合は正しい重心が計算できない fitsthresholddetector.cpp
SEP  SourceExtractor による星の検出 fitssepdetector.cpp
Bahtinov   ハフ変換で三本の光芒の位置を見つけ、その位置のズレを計算する fitsbahtinovdetector.cpp

これらの計算はFITSData::findStars()メソッドで呼び出されます。FITSDataクラスはFITS画像データの保持及びその画像処理を行うクラスです。どのdetectorもすべてFITSStarDetectorクラスを継承していて、同じ呼び出しかたで星の位置を算出することができます。

せっかく同じクラスを継承しているのに、switch文でdetectorごとに処理を分けているのが今ひとつ美しくないですが。。
QFuture FITSData::findStars(StarAlgorithm algorithm, const QRect &trackingBox)
{
    if (m_StarFindFuture.isRunning())
        m_StarFindFuture.waitForFinished();

    starAlgorithm = algorithm;
    qDeleteAll(starCenters);
    starCenters.clear();
    starsSearched = true;

    switch (algorithm)
    {
        case ALGORITHM_SEP:
        {
            QPointer detector = new FITSSEPDetector(this);
            detector->setSettings(m_SourceExtractorSettings);
            if (m_Mode == FITS_NORMAL && trackingBox.isNull())
            {
                if (Options::quickHFR())
                {
                    //Just finds stars in the center 25% of the image.
                    const int w = getStatistics().width;
                    const int h = getStatistics().height;
                    QRect middle(static_cast(w * 0.25), static_cast(h * 0.25), w / 2, h / 2);
                    return detector->findSources(middle);
                }
            }
            m_StarFindFuture = detector->findSources(trackingBox);
            return m_StarFindFuture;
        }

        case ALGORITHM_GRADIENT:
        default://defaultはGradient
        {
            QPointer detector = new FITSGradientDetector(this);
            detector->setSettings(m_SourceExtractorSettings);
            m_StarFindFuture = detector->findSources(trackingBox);
            return m_StarFindFuture;
        }

        case ALGORITHM_CENTROID:
        {
#ifndef KSTARS_LITE
            QPointer detector = new FITSCentroidDetector(this);
            detector->setSettings(m_SourceExtractorSettings);
            // We need JMIndex calculated from histogram
            if (!isHistogramConstructed())
                constructHistogram();
            detector->configure("JMINDEX", m_JMIndex);
            m_StarFindFuture = detector->findSources(trackingBox);
            return m_StarFindFuture;
        }
#else
            {
                QPointer detector = new FITSCentroidDetector(this);
                return detector->findSources(starCenters, trackingBox);
            }
#endif

        case ALGORITHM_THRESHOLD:
        {
            QPointer detector = new FITSThresholdDetector(this);
            detector->setSettings(m_SourceExtractorSettings);
            detector->configure("THRESHOLD_PERCENTAGE", Options::focusThreshold());
            m_StarFindFuture =  detector->findSources(trackingBox);
            return m_StarFindFuture;
        }

        case ALGORITHM_BAHTINOV:
        {
            QPointer detector = new FITSBahtinovDetector(this);
            detector->setSettings(m_SourceExtractorSettings);
            detector->configure("NUMBER_OF_AVERAGE_ROWS", Options::focusMultiRowAverage());
            m_StarFindFuture = detector->findSources(trackingBox);
            return m_StarFindFuture;
        }
    }
}

どの方法が良いのかやってみないとわかりませんが、GradientとCentroidはどっちもどっちな感じがします。複数の星の位置が同時に欲しい場合はCentroidになります。またCentroidでは変に大きな領域がある場合はその領域を省く処理が入っているので、星以外のものが写り込んでいるような場合にも有効そうです。

また星の大体の場所がわかっていて、指定領域内に星が一つしかないことが保証されているようなシチュエーションだと、Threthioldで十分な気がします。

●では次、2)星像の大きさの計算です。
これは定番のHFD(EkosではHFR)を計算します。HFDは Half Flux Diameterの略。Diameterなので星の直径を表します。EkosはRadius(半径)なのでHFR。HFDと意味は同じです。
HFDについてはこちらが詳しいです。



HFDの計算方法は、中心からの距離を重みとした輝度の重み付き平均を、単なる輝度の平均で割ったもの、です。式は下の通り。Viは各画素の輝度、diは各画素の重心からの距離です。
halffluxdiameter_html_m3052dbde

その意味は、星の領域とその外の領域とが、輝度の重み付き和が等しくなる直径。なので輝度が一点に集中して明るい方が、値が小さくなります。また星が検出されない場合は、その値は計算範囲の1/√2 となります。

HFD(HFR)の計算は、(1)の星を検出したときにそれぞれのdetectorのfindsources()の中の終わりの方で計算されています。Bahtinovを選んだ場合は、光芒の三本線の中央線のズレをHFRの値として保存しています。ですのでBahtinovではHFRの意味が異なりますが、フォーカスが合っているほど値が小さくなる、という意味では同じで、このあとのフォーカスシフトによるAFの評価値として、星像の大きさと同じように扱うことができます。ただBahtinovマスクの長所は中央の光芒がどちらにずれているかでピントのズレ方向がわかるのが特長ですが、それを使っていないのはなんかもったいない気がします。


●さて最後は安定した星像の大きさの選択です。星像の大きさはシンチレーションや風、振動などの影響を受けます。そのためEkosではHFRの値が±2σに入る値が指定の数蓄積されるまでキャプチャを続ける、ということをやっています。

やっている場所はFocus::appendHFR()。focus.cppの中にあります。

bool Focus::appendHFR(double newHFR)
{
    // Add new HFR to existing values, even if invalid
    HFRFrames.append(newHFR);

    // Prepare a work vector with valid HFR values
    QVector  samples(HFRFrames);
    samples.erase(std::remove_if(samples.begin(), samples.end(), [](const double HFR)
    {
        return HFR == FocusAlgorithmInterface::IGNORED_HFR;
    }), samples.end());

    // Perform simple sigma clipping if more than a few samples
    if (samples.count() > 3)
    {
        // Sort all HFRs and extract the median
        std::sort(samples.begin(), samples.end());
        const auto median =	//中央値を計算
            ((samples.size() % 2) ?
             samples[samples.size() / 2] :
             (static_cast(samples[samples.size() / 2 - 1]) + samples[samples.size() / 2]) * .5);

        // Extract the mean
        const auto mean = std::accumulate(samples.begin(), samples.end(), .0) / samples.size();

        // Extract the variance
        double variance = 0;
        foreach (auto val, samples)
            variance += (val - mean) * (val - mean);

        // Deduce the standard deviation
        const double stddev = sqrt(variance / samples.size());//標準偏差を計算

        // Reject those 2 sigma away from median//中央値から2σ外れている値は除外
        const double sigmaHigh = median + stddev * 2;
        const double sigmaLow  = median - stddev * 2;

        // FIXME: why is the first value not considered?
        // FIXME: what if there are less than 3 samples after clipping?
        QMutableVectorIterator i(samples);
        while (i.hasNext())
        {
            auto val = i.next();
            if (val > sigmaHigh || val < sigmaLow)
                i.remove();
        }
    }

    // Consolidate the average HFR
    currentHFR = samples.isEmpty() ? -1 : std::accumulate(samples.begin(), samples.end(), .0) / samples.size();

    // Return whether we need more frame based on user requirement
    return HFRFrames.count() < focusFramesSpin->value();
}

ここで一定数2σ内にHFRの値が入るようになったあと、実際に使うのは得られたHFRの平均とかではなく最新の値だけを使っています。

やっていることは単純ですが、これでうまく行っているようです。