WIP: Bayes estimation

.vscode/c_cpp_properties.json

@@ -6,7 +6,8 @@
                 "${workspaceFolder}/**",
                 "${workspaceFolder}/src/vstarstack/library/fine_movement/libimagedeform/include",
                 "/usr/include/python3.11",
-                "/usr/include/python3.13"
+                "/usr/include/python3.13",
+                "/home/vlad/Astronomy/software/vstarstack/lib/python3.13/site-packages/numpy/_core/include"
             ],
             "defines": [],
             "compilerPath": "/usr/bin/gcc",

.vscode/settings.json

@@ -43,7 +43,8 @@
         "span": "cpp",
         "stdexcept": "cpp",
         "typeinfo": "cpp",
-        "imagegrid.h": "c"
+        "map": "c",
+        "format": "c"
     },
     "editor.tabSize": 4,
     "editor.insertSpaces": true,

doc/deepsky.tex

@@ -0,0 +1,303 @@
+\documentclass{article}
+\usepackage{mathtools}
+
+\begin{document}
+
+\section{Task}
+
+We have frames $F(x,y)$:
+
+\begin{eqnarray}
+    F(x,y) \sim Poisson\left( \lambda(x,y) \right) \\
+    \lambda(x,y) = f_{dark}(x,y) + v(x,y) \cdot \left( f_{sky}(x,y) + f_{true}(x,y) \right)
+\end{eqnarray}
+
+Where 
+\begin{itemize}
+    \item $f_{dark}(x,y)$ --- dark signal
+    \item $f_{sky}(x,y)$ --- sky signal (light pollution)
+    \item $f_{true}(x,y)$ --- true signal
+    \item $v(x,y)$ --- flat relative calibration function (vignetting, particles on sensor, etc)
+\end{itemize}
+
+What does "relative" means? We have telescope lens or mirror, sensor with some QE, exposure, sensor gain. We have Poisson--distributed amount
+of photons and dark electrons on sensor, which is multiplied by x,y--independent factor of QE, lens size, exposure, and x,y--dependent factor
+of vignetting, particles on sensor, etc. We want to compensate x,y--dependent factors. So we can assume $max(v) = 1$, because it shows
+relative coefficient to x,y--maximal value. If $max(v)$ would be, for example, $0.5$, it would be same as twice less QE with $max(v) = 1$.
+But we know that it still will be Poisson distribution.
+
+We want to estimate true value of signal: $f_{true}(x,y)$.
+
+\section{Calibrations}
+
+\subsection{Dark signal}
+
+We can find $f_{dark}(x,y)$ by taking enough dark frames $D_k(x,y)$ and find mean of them.
+In later research we assume that we know $f_{dark}(x,y)$ with absolute precision.
+
+\begin{equation}
+    D_k(x,y) \sim Poisson(f_{dark}(x,y))
+\end{equation}
+
+So we find
+
+\begin{equation}
+    f_{dark}(x,y) = \frac{\sum_{k=1}^K D_k(x,y)}{K}
+\end{equation}
+
+\subsection{Flats}
+
+We need to calculate $v(x,y)$. We can either use flatbox or capture enough images of sky in different positions.
+Building flat is not a subject of current article. In later we assume that we know $v(x,y)$ with absolute precision.
+
+\section{Signal obtaining}
+
+As input we have:
+\begin{itemize}
+    \item $d(x,y)$ --- constant signal (dark or dark + sky, depending on case)
+    \item $v(x,y)$ --- normalized flats function
+    \item $F_j(x,y)$ --- captured frames
+\end{itemize}
+
+We know that $F_i$ have Poisson distribution:
+\begin{eqnarray}
+    F_i(x,y) = Poisson\left( d(x,y) + v(x,y) \cdot f(x,y) \right)
+\end{eqnarray}
+
+How $d$ is related with $f_{dark}$, $f_{sky}$? It depends on approach.
+In most basic approach we do 2--step calculation:
+\begin{itemize}
+    \item Capture dark frames
+    \item Capture flat frames
+    \item Capture target frames
+    \item Considering $d = f_{dark}$ find image signal $f_{full}$. It is sum of $f_{sky}$ and $f$
+    \item Analyze $f_{full}$ and find $f_{sky}$ from it
+    \item Considering $d = f_{dark} + v \cdot f_{sky}$ find target image signal $f$
+\end{itemize}
+
+In case of high light pollution we even don't need bayesian calculation to obtain $f_{sky}$
+
+\subsection{First approach}
+
+We act the same way as for flats calculation:
+\begin{eqnarray}
+    <F_i(x,y)> = d(x,y) + v(x,y) \cdot f(x,y) \\
+    f(x,y) = \frac{<F_i(x,y)> - d(x,y)}{v(x,y)} = \frac{<F_i(x,y) - d(x,y)>}{v(x,y)}
+\end{eqnarray}
+
+Here we use Sigma-Clipping:
+
+\begin{eqnarray}
+    f(x,y) = \frac{SC(\left\{ F_i(x,y) - d(x,y)\right\})}{v(x,y)}
+\end{eqnarray}
+
+This implements the most basic approach: substract darks, then apply flats, then do Sigma-Clipping.
+This method is simple, but can generate negative values, especially for low SNR case.
+
+\subsection{Bayesian analysis}
+
+We have set of frames $F_i(x,y)$. We want to find estimation for $f(x,y)$:
+
+\begin{eqnarray}
+    E\left[ f(x,y) \right] = \int df(x,y) \cdot f(x,y) \cdot p(f(x,y) | \left\{ F_i(x,y) \right\})
+\end{eqnarray}
+
+Let's find $p(f(x,y) | \left\{ F_i(x,y) \right\})$. We process independently for each $x,y$, so they are just parameters,
+not integration varables. So we won't write them for readability:
+
+\begin{eqnarray}
+    E\left[ f \right] = \int df \cdot f \cdot p(f | \left\{ F_i \right\})
+\end{eqnarray}
+
+According to Bayes' theorem
+
+\begin{eqnarray}
+    p(f | \left\{ F_i \right\}) = \frac{p(\left\{ F_i \right\} | f) \cdot p(f)}{p(\left\{ F_i \right\})} =
+    \frac{p(\left\{ F_i \right\} | f) \cdot p(f)}{\int p(\left\{ F_i \right\} | f') \cdot p(f') \cdot df'} \label{eq:bayes}
+\end{eqnarray}
+
+where $p(f)$ --- apriori probability of $f$. We should estimate it before starting Bayesian analysis.
+
+It's very important that
+\begin{equation}
+    p(\left\{ F_i \right\}) \ne \prod_i p(F_i)
+\end{equation}
+
+It happens because $F_i$ are not marginally independent: if $P(AB|C_i) = P(A|C_i) \cdot P(B|C_i)$ for each $i$, it doesn't mean that $P(AB) = P(A)\cdot P(B)$
+
+But $F_i|f$ are conditionally independent, so we can write
+\begin{equation}
+    p(\left\{ F_i \right\} | f) = \prod_i p(F_i | f)
+\end{equation}
+
+So we can continue (\ref{eq:bayes}) and write it as
+\begin{eqnarray}
+    p(f | \left\{ F_i \right\}) = \frac{\prod_i p(F_i|f) \cdot p(f)}{\int \prod_i p(F_i|f') \cdot p(f') df'} \label{eq:bayes2}
+\end{eqnarray}
+
+Let's remember what $p(F_i | f)$ is:
+\begin{equation}
+\begin{array}{r@{}l}
+    p(F_i | f) = \frac{\lambda^{F_i} \exp(-\lambda)} {F_i!} \\
+    \lambda = d(x,y) + v(x,y) \cdot f(x,y)
+\end{array}
+\end{equation}
+
+Let's substutute it into (\ref{eq:bayes2}):
+\begin{equation}
+\begin{array}{r@{}l}
+    p(f | \left\{ F_i \right\}) = \frac{ \prod_i \frac{\lambda(f)^{F_i} \cdot \exp(-\lambda(f))}{F_i!} \cdot p(f)  }
+    {\int \prod_i \frac{\lambda(f)^{F_i} \cdot \exp(-\lambda(f))}{F_i!} p(f') df'}
+\end{array}
+\end{equation}
+
+After simplyfying it we have:
+\begin{equation}
+    p(f | \left\{ F_i \right\}) = \frac{\lambda(f)^{\sum F_i} \cdot \exp(-N\lambda(f)) \cdot p(f)}
+    {\int \lambda(f')^{\sum F_i} \cdot \exp(-N\lambda(f')) \cdot p(f') df'}
+\end{equation}
+
+where $N$ --- amount of $F_i$. Factorials are gone. And for case where $p(f) \ne 0$, $\lambda(f) \ne 0$ we can simplify:
+
+\begin{equation}
+    p(f | \left\{ F_i \right\}) = \left( \int \left( \frac{\lambda(f')}{\lambda(f)} \right)^{\sum F_i} \cdot
+                                        \exp\left\{-N \cdot (\lambda(f')-\lambda(f)) \right\} \cdot \left( \frac{p(f')}{p(f)} \right) df'\right)^{-1}
+\end{equation}
+
+But there can be some practical case: what if different frames have different $f_{dark}$ or $f_{sky}$ or $\nu$? In that case we have series
+of $\lambda_i(f)$ and formula slightly changes:
+\begin{equation}
+    p(f | \left\{ F_i \right\}) = \left( \int \prod_i \left( \frac{\lambda_i(f')}{\lambda_i(f)} \right)^{F_i} \cdot 
+    \exp\left\{-\sum_i \left(\lambda_i(f')-\lambda_i(f)\right)\right\} \cdot \left( \frac{p(f')}{p(f)} \right) df'\right)^{-1}
+\end{equation}
+
+In this formula we have 2 components $\prod_i \left( \frac{\lambda_i(f')}{\lambda_i(f)} \right)^{F_i}$ and
+$\exp\left\{-\sum_i \left(\lambda_i(f')-\lambda_i(f)\right) \right\}$. For big $f'$ first would be big, but it will be compensated by exponent.
+But it can cause overflow. So we can write:
+
+\begin{eqnarray}
+    p(f | \left\{ F_i \right\}) = \left( \int \exp\left\{\sum_i \left( F_i ln \lambda_i(f') - \lambda_i(f')\right) 
+            - \left( F_i ln \lambda_i(f) - \lambda_i(f)\right) \right\} \cdot \frac{p(f')}{p(f)} df' \right)^{-1} \label{eq:bayes3}
+\end{eqnarray}
+
+Or, introducing $\Lambda_i$:
+
+\begin{eqnarray}
+    \Lambda_i(f) = F_i ln \lambda_i(f) - \lambda_i(f) \\[9pt]
+    p(f | \left\{ F_i \right\}) = \left( \int \exp\left\{\sum_i \Lambda_i(f') - \Lambda_i(f) \right\}\cdot \frac{p(f')}{p(f)} df' \right)^{-1}
+\end{eqnarray}
+
+Also we can use that formula in case when $f$ is a vector. But we have to adjust $\lambda(f)$ in that case:
+
+\begin{eqnarray}
+    \lambda_i (f_{\alpha}) = d_{i} +  v_{i} \cdot \sum_{\beta} K_{i, \alpha} \cdot f_{\alpha}
+\end{eqnarray}
+
+Here:
+\begin{itemize}
+\item $i$ --- index for frames
+\item $\alpha$ --- index for $f$ dimensions
+\item $K$ --- matrix for contribution of each component into common signal
+\end{itemize}
+
+Challenges with such approach:
+\begin{itemize}
+\item we need some a-priori knowledge of probability $p(f)$
+\item integral calculation over $R_{+}$ is pretty long procedure. So we can calculate it only in some area around maxima of $p(f')$
+\end{itemize}
+
+Good point is that for some priors, integral can be taken analytically.
+
+\subsection{Pipeline}
+So the pipeline is following:
+
+\begin{itemize}
+    \item Perform a "naive" estimation to obtain initial $f$
+    \item Using a sliding window, build local histograms of $f$ to estimate the local prior $p(f)$. Normalize and optionally smooth these histograms.
+    \item Use the Bayesian estimator for $E\left[ f \right]$ to refine $f$, based on the current prior
+    \item Iterate the prior estimation and refinement steps until convergence or acceptable error level
+\end{itemize}
+
+\section{True signal obtain}
+
+During previous step we studied how to obtain signal $f$. But it contains sky light pollution:
+\begin{eqnarray}
+    f = f_{sky} + f_{signal} \\
+    \lambda(f_{signal}) = f_{dark} + \nu \left( f_{sky} + f_{signal} \right)
+\end{eqnarray}
+
+So, we need to estimate $f_{sky}$.
+
+\subsection{Sky estimation}
+Firstly, let's using formula $\lambda(f) = f_{dark} + \nu f$ find whole captured signal $f$.
+
+We can estimate $f_{sky}$ from $f$ using regions on the side of image, where target signal is expected to be absent.
+But we should be careful, not to overestimate.
+
+\subsection{Bayesian analysis}
+
+So, now we use $\lambda(f_{signal}) = f_{dark} + \nu \left( f_{sky} + f_{signal} \right)$ to find $f_{signal}$.
+
+\subsection{Pipeline}
+So the pipeline is following:
+
+\begin{itemize}
+    \item obtain whole signal $f(x,y)$
+    \item find $f_{sky}$ from $f(x,y)$
+    \item find expected $f_{signal}$
+\end{itemize}
+
+\section{Narrowband signal extraction}
+
+In this case we have 2 signal images, captured with narrow (1) and wide (2) filters. Whole signal has 2 components: $f_{c,*}$ --- continuum,
+and $f_{n,*}$ --- narrow signal:
+
+\begin{eqnarray}
+    f_{signal, 1} = f_{c, 1} + f_{n} \\
+    f_{signal, 2} = f_{c, 2} + f_{n}
+\end{eqnarray}
+
+And we assume that continuum for both cases is proportional. We can assume that if wide filter not too wide.
+\begin{equation}
+    f_{c,2} = k \cdot f_{c,1} = k \cdot f_c
+\end{equation}, where $k > 1$. Value $k$ we obtain from calibration by stars --- they emit only continuum.
+
+So
+\begin{eqnarray}
+    f_{signal, 1} = f_c + f_n \\
+    f_{signal, 2} = k \cdot f_c + f_n
+\end{eqnarray}
+
+So we can find $f_n$ as
+\begin{eqnarray}
+    (1-1/k) \cdot f_n = f_{signal, 1} - \frac{f_{signal, 2}}{k} \\
+    f_n = \frac{k \cdot f_{signal, 1} - f_{signal, 2}}{k-1} \\
+    f_c = f_{signal,1} - f_n = \frac{f_{signal, 2} - f_{signal, 1}}{k-1}
+\end{eqnarray}
+
+But as in previous parts, it's also a first approach.
+
+\subsection{Bayesian analysis}
+
+\begin{eqnarray}
+    F_{i,1} = Poisson\left( \lambda_{i,1}\left({f_n \atop f_c}\right) \right) \\
+    F_{j,2} = Poisson\left( \lambda_{j,2}\left({f_n \atop f_c}\right) \right)
+\end{eqnarray}
+
+\begin{eqnarray}
+    \lambda_{i,1}\left({f_n \atop f_c}\right) = f_{dark, 1, i} + \nu_{1,i} \cdot \left( f_{sky, 1, i} + f_c + f_n \right) \\
+    \lambda_{j,2}\left({f_n \atop f_c}\right) = f_{dark, 2, j} + \nu_{2,j} \cdot \left( f_{sky, 2, j} + k \cdot f_c + f_n \right)
+\end{eqnarray}
+
+We can concatinate $\left\{F_{i,1}, F_{j,2}\right\}$ frames to single frames list and
+$\left\{ \lambda_{i,1}, \lambda_{j,2} \right\}$ to single $\lambda$ list, and like in (\ref{eq:bayes3}) we get:
+
+\begin{eqnarray}
+p\left({f_n \atop f_c} | \left\{ F_{i,1}, F_{j,2} \right\}\right) =
+\left( \int e^{
+    -\sum_i \left[\Lambda_{i,1}\left({f'_n \atop f'_c}\right) - \Lambda_{i,1}\left({f_n \atop f_c}\right) \right]
+    -\sum_j \left[\Lambda_{j,2}\left({f'_n \atop f'_c}\right) - \Lambda_{j,2}\left({f_n \atop f_c}\right) \right]}
+\cdot \frac{p\left({f'_n \atop f'_c}\right)}{p\left({f_n \atop f_c}\right)} df'_n df'_c \right)^{-1}
+\end{eqnarray}
+
+\end{document}

pyproject.toml

@@ -4,7 +4,11 @@ build-backend = "py_build_cmake.build"
 
 [project]
 name = "vstarstack"
+<<<<<<< HEAD
 version = "0.3.7.11"
+=======
+version = "0.3.8"
+>>>>>>> 4006f6e (Test bayes)
 requires-python = ">=3.10"
 authors = [{name="Vladislav Tsendrovskii", email="vtcendrovskii@gmail.com"}]
 description = 'Stacking astrophotos'
@@ -36,10 +40,11 @@ directory = "src/"
 name = "vstarstack"
 
 [tool.py-build-cmake.cmake]
+build_type = "RelWithDebInfo"
 source_path = "src/c_modules" 
 install_components = [
     "python_modules"
 ]
 
 [tool.py-build-cmake.sdist]
-include = ["src/c_modules/*"]
+include = ["src/c_modules/*"]

src/c_modules/CMakeLists.txt

@@ -1,6 +1,11 @@
 cmake_minimum_required(VERSION 3.10)
 project(vstarstack)
 
+if (BUILD_TESTS)
+    enable_testing()
+    include(CTest)
+endif()
+
 find_package(Python3 REQUIRED COMPONENTS Interpreter Development.Module)
 
 message("Python: '${Python3_EXECUTABLE}'")
@@ -12,10 +17,10 @@ if(NUMPY_NOT_FOUND)
     message(FATAL_ERROR "NumPy headers not found: ${NUMPY_NOT_FOUND}")
 endif()
 
-
 message("Numpy headers: '${NUMPY_INCLUDE_DIR}'")
 
 add_subdirectory(projections)
 add_subdirectory(movements)
 add_subdirectory(fine_movement)
 add_subdirectory(clusterization)
+add_subdirectory(bayes)

src/c_modules/bayes/CMakeLists.txt

@@ -0,0 +1,13 @@
+cmake_minimum_required(VERSION 3.10)
+
+add_subdirectory(libbayes)
+
+Python3_add_library(bayes_py MODULE module.c WITH_SOABI)
+
+target_include_directories(bayes_py PUBLIC lib "${NUMPY_INCLUDE_DIR}")
+target_link_libraries(bayes_py PUBLIC bayes)
+set_target_properties(bayes_py PROPERTIES OUTPUT_NAME "bayes")
+
+install(TARGETS bayes_py
+        COMPONENT python_modules
+        DESTINATION ${PY_BUILD_CMAKE_IMPORT_NAME}/library/bayes)

src/c_modules/bayes/libbayes/CMakeLists.txt

@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 3.10)
+project(bayes)
+
+add_library(bayes STATIC bayes.c)
+target_include_directories(bayes INTERFACE .)
+
+if (BUILD_TESTS)
+    enable_testing()
+    include(CTest)
+    add_subdirectory(tests)
+endif()