star-gs

Literate program for a geometric sensitivity calculation
git clone git://git.meso-star.fr/star-gs.git
Log | Files | Refs | README | LICENSE

commit 2edab611144a9c9d1524b0ce07c54098dd48f13e
parent 1800cfa83842618f00890aeaf5a0b5088ca2fdfa
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed, 12 Apr 2023 22:04:55 +0200

Merge remote-tracking branch 'origin/feature_trans_sensib' into feature_trans_sensib

Diffstat:
Msrc/sgs_compute_sensitivity_translation.nw | 154++++++++++++++++++++++++++++++++++++++++++++++---------------------------------
1 file changed, 89 insertions(+), 65 deletions(-)

diff --git a/src/sgs_compute_sensitivity_translation.nw b/src/sgs_compute_sensitivity_translation.nw @@ -770,6 +770,14 @@ struct sgs_hit hit1; \section{Résultats} \label{sec:results} +Au delà du simple calcul de sensibilité nous proposons ici une vérification des +résultats par le calcul de différences finies. La différentiation est effectuée +à partir des estimations du flux $\varphi(\PI)$ reçu par le capteur pour +différentes valeurs du paramètre géométrique $\PI$, soit pour différentes +hauteurs de la paroi spéculaire. Le calcul du poids associé au calcul du flux +est décrit en annexe \ref{flux}. La figure \ref{fig:resultats} présente les +estimations de la sensibilité du flux et des différences finies correspondantes. + \begin{figure}[h!] \centering \begin{tikzpicture} @@ -818,21 +826,19 @@ struct sgs_hit hit1; \end{tikzpicture} \flushleft - \caption{Sensibilité du flux reçu par le capteur de surface $A_r$ en fonction - de l'ouverture du parallélépipède paramétrée par $\PI$. La sensibilité du flux - $\partial_{\PI} \hat{\varphi} = \frac{\partial_{\PI} - \varphi}{\varphi_{max}} h$ avec $\varphi_{max} = \varphi(\PI=0)$.} + \caption{Estimations des sensibilités du flux par Monte Carlo et comparaison + avec les différences finies du flux. Les résultats sont adimensionnalisés et + affichés en fonction de l'ouverture du parallélépipède, qui est paramétrée + par $\PI$. La sensibilité du flux $\partial_{\PI} \hat{\varphi} = + \frac{\partial_{\PI} \varphi}{\varphi_{max}} h$ avec $\varphi_{max} = + \varphi_{spec}(\PI=0)$. Dans cette expression, $\varphi_{spec}$ correspond + uniquement à la partie du flux qui arrive au récepteur après avoir été + réfléchie sur la paroi spéculaire (voir annexe \ref{flux}). Le nombre + d'échantillonnage du poids de Monte Carlo nécessaire à la reproduction de ce + résultat est $10^{8}$.} + \label{fig:resultats} \end{figure} -\paragraph{TODO} À écrire qd on aura le programme fonctionnel. - -\paragraph{temp} Au delà du simple calcul de sensibilité nous proposons ici une -vérification des résultats par le calcul de différences finies. La -différentiation est effectuée à partir des estimations du flux reçu par le -capteur pour différente valeur du paramètre géométrique $\PI$, soit pour -différentes hauteurs de la paroi spéculaire. Le calcul du poids associé au -calcul du flux est décrit en annexe \ref{flux}. - \appendix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -844,8 +850,8 @@ calcul du flux est décrit en annexe \ref{flux}. Nous récupérons ici la condition à la limite pour une paroi spéculaire donnée dans \cite{papier_sensib} $s = s(\vec{x},\vec{\omega},\PI)$: \begin{equation} -s = \mathcal{C}_b \lbrack s \rbrack + S_{b,\PI} \lbrack I, \partial_{1,\vec{u}}I, -\partial_{1,\vec{\chi}} I, \partial_{2,\vec{\gamma}_t}I \rbrack \quad \quad \quad +s = \mathcal{C}_b \lbrack s \rbrack + S_{b,\PI} \lbrack L, \partial_{1,\vec{u}}L, +\partial_{1,\vec{\chi}} L, \partial_{2,\vec{\gamma}_t}L \rbrack \quad \quad \quad \vec{x} \in \partial G(\PI) ; \vec{\omega} \cdot \vec{n} > 0 \label{eq:cl_sensib_gen} \end{equation} @@ -865,11 +871,17 @@ p_{\Omega'}(-\vec{\omega}'|\vec{x},-\vec{\omega})d\vec{\omega}' L \\ \end{aligned} \end{equation} -Dans cette équation $\mathcal{C}$ est l'opérateur collisionnel du milieu décrit -dans l'équation \ref{eq:C_operator}, il est ici appliqué à la luminance. La -source $S$ est la source radiative du milieu. On trouve également -$\mathcal{C}_b$ l'opérateur collisionnel de la surface. Appliqué à la -luminance et dans le cas d'une paroi spéculaire il s'écrit: +Dans cette équation $\mathcal{C}$ est l'opérateur collisionnel du milieu: +\begin{equation} +\mathcal{C}[L] = -k_a L(\vec{x},\vec{\omega},\PI) - k_s +L(\vec{x},\vec{\omega},\PI) + k_s \int_{\mathcal{S}} +p_{\Omega'}(\vec{\omega}'|\vec{x},\vec{\omega}) d\vec{\omega}' +L(\vec{x},\vec{\omega}',\PI) +\end{equation} +La source $S$ est la source radiative du milieu (émission thermique $k_a +L^{eq}(T)$). On trouve également $\mathcal{C}_b$ l'opérateur collisionnel de la +surface. Appliqué à la luminance et dans le cas d'une paroi spéculaire il +s'écrit: \begin{equation} \mathcal{C}_b[L] = \rho(\vec{x},-\vec{\omega}) \int_{H'} \delta(\vec{\omega}' - \vec{\omega}_{spec}) L(\vec{x},\vec{\omega}',\PI) d\vec{\omega}' @@ -878,22 +890,24 @@ avec $\vec{\omega}_{spec} = \vec{\omega} - 2(\vec{\omega} \cdot \vec{n})\vec{n}$ \paragraph{Condition à la limite de notre exemple} -Pour commencer seule la paroi spéculaire est source de sensibilité géométrique. -Nous voyons que dans l'équation \ref{eq:cl_sensib_gen} le terme collisionnel +Pour commencer, seule la paroi spéculaire est source de sensibilité géométrique. +Nous voyons dans l'équation \ref{eq:cl_sensib_gen} que le terme collisionnel $\mathcal{C}_b[s]$ traduit la réflexion de la sensibilité incidente à la paroi spéculaire. Ce terme est obligatoirement nul puisque il n'y a aucune autre source qui pourrait émettre une sensibilité géométrique et aucune autre paroi réfléchissante qui pourrait réfléchir la sensibilité émise par la paroi spéculaire. -Dans notre exemple le milieu est transparent, les termes $\mathcal{C}[L]$ et -$S$ sont donc nuls. La paroi spéculaire est froide, la source surfacique $S_b$ -qui dans cet exemple correspondrait à l'émission thermique de la paroi est donc -aussi nulle. L'opérateur collisionnel de la surface $\mathcal{C}_b$ est -indépendant de $\PI$, la dérivée $\partial_{\PI} \mathcal{C}_b$ est donc nulle. -Pour finir la déformation géométrique de la paroi spéculaire est une -translation, l'axe de rotation $\vec{\gamma}$ est donc nul et toutes les -dérivées angulaires n'ont plus lieux d'être dans la condition aux limites. +Dans notre exemple le milieu est transparent, les termes +$\mathcal{C}[L]$ et $S$ sont donc nuls. La paroi spéculaire est +froide, la source surfacique $S_b$ qui dans cet exemple +correspondrait à l'émission thermique de la paroi est donc aussi +nulle. L'opérateur collisionnel de la surface $\mathcal{C}_b$ est +indépendant de $\PI$, la dérivée $\partial_{\PI} \mathcal{C}_b$ est +donc nulle. Pour finir la déformation géométrique de la paroi +spéculaire est une translation, l'axe de rotation $\vec{\gamma}$ est +donc nul et toutes les dérivées angulaires n'ont plus lieux d'être +dans la condition à la limite. La condition à la limite de sensibilité de la paroi spéculaire de la boite devient donc: @@ -932,17 +946,19 @@ L(\vec{x},\vec{\omega}_{spec},\PI) \\ \section{Décomposition du vecteur de déformation} \label{ann:proj} -Dans le modèle de sensibilité la déformation est caractérisée par le vecteur de -déformation $\vec{\chi}$. La condition à la limite de sensibilité dépend alors -de la dérivée spatiale $\partial_{1,\vec{\chi}}I$ (équation \ref{eq:clsensib}). -Le champs de luminance n'étant pas connu il n'existe pas de solution analytique -à cette dérivée. À la frontière nous avons donc choisi de décomposer la -direction $\vec{\chi}$ ([[chi]]) en deux directions distinctes, une direction -tangente à la frontière $\vec{u}$ ([[u]]) et la direction de transport -$\vec{\omega}$ ([[omega]]). Ainsi la dérivée de $I$ selon $\vec{\chi}$ devient -une composition de dérivées de $I$ le long de la paroi (sur laquelle la luminance -est connue) et le long de la direction de transport (retrouvant ainsi le terme -de transport de l'ETR). +Dans le modèle de sensibilité la déformation est caractérisée par le +vecteur de déformation $\vec{\chi}$. La condition à la limite de +sensibilité dépend alors de la dérivée spatiale +$\partial_{1,\vec{\chi}}L$ (équation \ref{eq:clsensib}). Le champs +de luminance n'étant pas connu il n'existe pas de solution +analytique à cette dérivée. À la frontière nous avons donc choisi de +décomposer la direction $\vec{\chi}$ ([[chi]]) en deux directions +distinctes, une direction tangente à la frontière $\vec{u}$ ([[u]]) +et la direction de transport $\vec{\omega}$ ([[omega]]). Ainsi la +dérivée de $L$ selon $\vec{\chi}$ devient une composition de +dérivées de $L$ le long de la paroi (sur laquelle la luminance est +connue) et le long de la direction de transport (retrouvant ainsi le +terme de transport de l'ETR). La décomposition de $\vec{\chi}$ s'écrit: \begin{equation} @@ -956,20 +972,24 @@ nous renvoyons le lecteur vers \cite{papier_sensib}. Nous nous contentons de donner ici les résultats: \begin{equation} \alpha = \frac{\vec{\chi}.\vec{n}}{\vec{\omega}.\vec{n}} \ \; \quad \beta = -\|\vec{chi} - \alpha \vec{\omega} \| \ \; \quad \vec{u} = \frac{\vec{chi} - +\|\vec{\chi} - \alpha \vec{\omega} \| \ \; \quad \vec{u} = \frac{\vec{\chi} - \alpha \vec{\omega}}{\beta} \end{equation} Ce qui permet d'obtenir: \begin{equation} -\partial_{1,\vec{\chi}}I = \alpha \partial_{1,\vec{\omega}} I + \beta -\partial_{1,\vec{u}} I +\partial_{1,\vec{\chi}}L = \alpha \partial_{1,\vec{\omega}} L + \beta +\partial_{1,\vec{u}} L \end{equation} -et de résoudre la dérivée en estimant $\partial_{1,\vec{u}}I$ le long de la -surface et en utilisant l'ETR pour résoudre $\partial{1,\vec{\omega}}I$: +et de résoudre la dérivée en estimant $\partial_{1,\vec{u}}L$ le long de la +surface et en utilisant l'ETR pour résoudre $\partial_{1,\vec{\omega}}L$: \begin{equation} -\partial_{1,\vec{\omega}}I = \mathcal{C} \lbrack I \rbrack +\partial_{1,\vec{\omega}}L = \mathcal{C} \lbrack L \rbrack \end{equation} +La fonction [[decomposition]] réalise cette décomposition et les résultats +([[alpha]], [[beta]] et [[u]]) sont répertoriés dans la structure +[[projection]]. + <<fonctions utilitaires>>= struct projection { double alpha; @@ -1002,32 +1022,31 @@ decomposition \section{Calcul de la contribution de la luminance qui dépend de $\PI$} \label{flux} -Le flux reçu par le capteur est décrit par l'équation \ref{eq:flux}. Étant -donné la configuration radiative de notre problème, la luminance -$L(\vec{x},\vec{\omega},\PI)$ incidente au récepteur ne dépend que de deux -types de chemins de transports de la source radiative qui est émise par la -paroi de droite: +Le flux reçu par le capteur est décrit par l'équation \ref{eq:flux}. Dans notre +problème, la luminance $L(\vec{x},\vec{\omega},\PI)$ incidente au récepteur +dépend de deux sortes de chemins radiatifs qui traduisent le transport de la +source émise par la paroi de droite : \begin{itemize} \item soit la source est transportée directement depuis la paroi de droite jusqu'au récepteur; - \item soit la source est transportée en direction de la paroi spéculaire et - est ensuite réfléchie jusqu'à atteindre le récepteur. + \item soit la source est transportée en direction de la paroi spéculaire puis + réfléchie jusqu'au récepteur. \end{itemize} Dans notre configuration le paramètre géométrique $\PI$ n'a d'influence que sur -la hauteur de la paroi spéculaire. Ainsi, la contribution radiative de la -source émise par la paroi de droite, transportée directement vers le -récepteur, reste identique pour toutes valeurs de $\PI$ (c'est à dire quelle que -soit la hauteur de la paroi spéculaire). Du point de vue de la sensibilité +la hauteur de la paroi spéculaire. Ainsi, la contribution radiative de la source +émise par la paroi de droite, transportée directement vers le récepteur, reste +identique pour toutes valeurs de $\PI$ (c'est à dire quelle que soit la hauteur +de la paroi spéculaire). Du point de vue de la sensibilité cette contribution n'aura donc aucune influence. En pratique cela signifie que, en vue d'une validation via un calcul des différences finies, le calcul par Monte Carlo du flux au récepteur n'est pas -entièrement nécessaire. Nous pouvons nous contenter d'estimer l'unique partie -du flux qui sera perturbée par une variation de $\PI$, soit les contributions -portées par les chemins réfléchis par la paroi spéculaire. En terme -d'algorithme nous pouvons donc réutiliser les chemins déjà échantillonnés pour -le problème couplé puisque leur statistique correspond exactement à celle d'un -chemin émis par la paroi de droite et réfléchi par la paroi du haut. +entièrement nécessaire. Nous pouvons nous contenter d'estimer l'unique partie du +flux qui sera perturbée par une variation de $\PI$, soit les contributions +portées par les chemins réfléchis par la paroi spéculaire. En terme d'algorithme +nous pouvons donc réutiliser les chemins déjà échantillonnés pour le problème +couplé puisque leur statistique correspond exactement à celle d'un chemin émis +par la paroi de droite et réfléchi par la paroi du haut. Le poids ([[weight_flux_part_spec]]) correspondant à la partie du flux venant de la réflection sur la paroi spéculaire est donc calculé au même moment que @@ -1037,6 +1056,11 @@ celui de la sensibilité. double weight_flux_part_spec; @ +Le poids de cette contribution correspond à la source radiative $S_b$ (equation +\ref{eq:cl_rad}) représentée par la variable [[Sb]] multipliées par l'angle +$\pi$ ([[PI]]), la surface $A_h$ et le coefficient de réflection $\rho$ +([[rho]]). + <<calcul du poids>>= weight_flux_part_spec = Sb * rho * PI * get_Sr(scene); @ @@ -1429,7 +1453,7 @@ indépendante de ses dimensions réelles. Dans notre script ces dimensions sont fixées via deux variables qui définissent le coin inférieur ([[lower]]) et le coin supérieur ([[upper]]) d'un parallélépipède aligné aux axes, deux variables passées en arguments -du programme via l'option [[-b]]. +du programme via l'option [[-b]]. <<lancer un calcul Monte Carlo>>= out=$(./sgs \ -n "${nrealisations}" \