Recent Changes - Search:

Menu

editer le Menu


Francis Dutil
Thomas Goulet-Soucy
David Kabanza
Maxence Mayrand
Sciences de la nature
Groupe 9140




Optimisation d'ILLUMINA





Rapport de laboratoire présenté à
M. Martin AUBÉ
dans le cadre du cours
Initiation à la recherche
Département de Physique





16 décembre 2011

Cégep de Sherbrooke





Introduction

Résumé

La pollution lumineuse nous affecte tous d'une façon ou d'une autre. Mais réussir à bien la calculer, et à bien la modéliser, n’est pas simple. Toutefois, avec ILLUMINA, un nouveau modèle prenant en considération bon nombre de variables, nous en sommes capables. Cela demande énormément de calcul. En effet, même sur un superordinateur comme le Mammouth série II, il nous est impossible de calculer plus qu'un ordre de diffusion pour des larges domaines géographiques de l'ordre de 10000 km2 et plus (résolution de 1 km2). Notre projet consiste donc à optimiser ILLUMINA en diminuant son temps de calcul pour ainsi pouvoir effectuer des calculs plus poussés (comme insérer la deuxième diffusion) à l'intérieur d'un temps raisonnable.

Étant donné la complexité du modèle, nous ne modifierons que les cartes d’entrées et de sorties. Pour ce faire, nous allons simplifier les cartes de luminosité, en attribuant à un certain point la luminosité totale de tous les points avoisinant, mettant ainsi leur valeur à zéro, ce qui allégera considérablement le temps de calcul, puisqu’ILLUMINA ne calcule pas les luminosités égales à zéro. Le nombre de points non-nuls diminuera ainsi en fonction de la distance au point d’observation. Nous donnerons toutefois plus d’importance aux endroits ayant une grande luminosité, même si cette zone se trouve plus éloignée. Puis une fois les calculs d’ILLUMINA effectués, nous redistribuons la somme des résultats aux points avoisinants selon leur proportion de la luminosité initiale.

Les résultats sont encourageant, ILLUMINA prenant nettement moins de temps de calculs, mais les paramètres du programme optimisant le produit du temps de calcul et de l'erreur commise par rapport au calcul complet restent à déterminer.

Dans un deuxième ordre d’idée, nous avons mis à jour les données intrants élaborés par Andréane D’Arcy-Lepage, Maude Fontaine et Caroline Pfister en 2010-2011 pour modéliser la réserve de ciel étoilé du Mont-Mégantic., qui avait grandement sous-estimé la présence de lanterne de ferme dans les milieux ruraux. Nous avons donc défini une nouvelle zone, la zone rurale, fait une étude sur ce type de lumière, et pris un échantillon de lampadaires dans la région de Stoke, pour en venir à la conclusion que le 2/3 des lampadaires ruraux sont des lanternes de ferme, De plus, selon les observations que Martin Aubé nous a fourni, 1/3 de ces lanternes avaient des ampoules au sodium haute pression et le 2/3 des ampoules au mercure. Il ne reste donc plus qu’à intégrer ces intrants aux autres déjà présentes, et de refaire une autre modélisation.

Abstract

Light pollution affects us all one way or another. But to succeed in calculating and modeling the quantity of light pollution present in a certain area is not simple since very few models of calculation are able to do so thoroughly. However, with ILLUMINA, a new model taking into account many variables, it is possible to calculate light pollution. However, this requires a lot of calculations. Even on a supercomputer such as Mammouth series II, it is impossible to calculate more than an order of distribution for large geographical areas of about 10,000 km2 and more. Our project is therefore to optimize ILLUMINA by reducing the computation time and power for more advanced calculations (like inserting the second distribution) within a reasonable time.

Given the complexity of the model, we will only change the card input and output. To do so, we will simplify the maps of brightness, by assigning to a certain point the total brightness of all points nearby, thus bringing the value of the nearby points to zero, which considerably reduces the computation time, does not calculate since ILLUMINA brightness equal to zero. The number of non-zero points will decrease depending on the distance between the point and the observation point. We will, however, place more emphasis on areas with a lot of light, even though some these areas can be farther away from the observation point. Then once ILLUMINA’s calculations are performed, we redistribute the sum of the results to neighboring points according to their proportion of contribution to the initial brightness.

The results are encouraging; ILLUMINA taking significantly less computing time, but the parameters of the optimizing program the results of the computation time and the error margin from the full calculation program need to be determined. We have also updated the data inputs developed by Andréane D'Arcy-Lepage, Maude Fontaine and Caroline Pfister in 2010-2011 to model the Dark Sky Preserve Mount Mégantic. They had greatly underestimated the presence of lantern farm in rural areas. We have defined a new zone, the rural area, and studied on its type of light, and took a sample of streetlights in the area of Stoke, to come to the conclusion that the 2 / 3 of rural street lights are farm lanterns. In addition, according to the observations provided to us by Martin Aubé, 1 / 3 of the lanterns had high pressure sodium bulbs and 2 / 3 of mercury bulbs. The only task left is to integrate these inputs to the other inputs already present, and make ILLUMINA execute new calculations.

Introduction

Lors des dernières années, Aubé et al. ont mis au point une méthode et un logiciel (ILLUMINA) capable de calculer la pollution lumineuse perceptible pour un site d'observation et une ligne de visée quelconque. Cette méthode a été appliquée au cas du Mont-Mégantic, la première réserve de ciel étoilé au monde, par Andréane D’Arcy-Lepage, Maude Fontaine et Caroline Pfister en 2010-2011. Toutefois, ce logiciel, prenant en considération plusieurs paramètres, comme la réflectance du sol, l'altitude, la hauteur des lampadaires, des obstacles, les différents types de lampadaires et possédant une très grande résolution, prend un temps énorme pour faire tout ses calculs, l'empêchant donc d’être utilisé à son plein potentiel. Le temps de calcul pour une seule diffusion est raisonable (2h environ sur le Mammouth serie II), toutefois elle est loin de l’être pour une deuxième diffusion puisque le temps de calcul nécessaire pour faire la modélisation dépasse notre temps alloué sur le superordinateur.

Nous voulons donc optimiser ILLUMINA pour pouvoir l'utiliser à son plein potentiel, et ainsi avoir des résultats qui tiennent comte d’un ordre de diffusion supérieur à un. Nous saurons si nous avons réussi si le temps gagné avec ces optimisations est notable, et si la perte de précisions vis-à-vis les données est en dessous d'un seuil critique que vous avons estimé à 1%.

Nous avons aussi peaufiné les données entrantes ayant servies à la modélisation de la pollution lumineuse au Mont-Mégantic en 2010 par Andréane D’Arcy-Lepage, Maude Fontaine et Caroline Pfister en 2010-2011. En effet, ces derniers n'avaient pas considéré la prédominance des lanternes de ferme dans les régions rurales, bien que ces lanternes soient parmi les plus polluantes.

Cadre théorique

Modèle ILLUMINA

Lorsqu’une lumière est dirigée vers le haut, que ce soit de façon directe ou par réflexion, elle risque, à un moment ou à un autre, de frapper un aérosol (petite particule de matière en suspension dans l’air) ou une molécule et d’être diffusée vers le sol. Ce phénomène s’appelle la pollution lumineuse. Mais d’où vient-elle? Cette lumière provient principalement des lampadaires, qui, trop souvent, diffusent inutilement de la lumière au dessus de l’horizon, là où personne ne peut en profiter. Cette lumière indésirable interfère donc avec la lumière naturelle de la nuit (des étoiles, de la lune …) et peut ainsi nuire à la faune, à la flore, mais également au travail des astronomes. C’est donc pour essayer de mieux cerner ce problème, que le modèle ILLUMINA fut créé.

Le modèle ILLUMINA permet de modéliser la pollution lumineuse pour un endroit donné, en affichant sa provenance (si elle vient des villages avoisinant, ou des villes plus éloignées...). Toutefois, ILLUMINA a ceci de spécial : contrairement à la plupart des autres modèles déjà existant dont plusieurs reposent sur celui de Garstang 1986, ILLUMINA est un modèle hétérogène (un autre modèle hétérogène est celui de Kocifaj 2007). Il permet donc aux villes de prendre n’importe quelle forme (puisque les villes ne sont généralement pas des cercles parfait), d’avoir une luminosité variable, une réflectance du sol variable, et une topographie variable. La seule condition étant que nous possédions toutes ces données pour la zone étudiée.

ILLUMINA calcule pour chaque pixel de la carte (représentant l’ensemble des lampadaires pour cette région), à partir du point d’observation, et ce pour chaque couche de l’atmosphère, la quantité théorique de pollution lumineuse pouvant venir de chaque région. Ainsi le calcul indiquera l'importance d'un pixel donné s’il y a une montagne entre les deux points (bloquant ainsi la lumière), dépendamment du type de lampadaire se trouvant dans cette zone (par exemple, une lanterne de ferme, qui émet environ 13% de sa lumière vers le haut comparativement a 6% pour un lampadaire normal de la réflectance du sol De plus, la hauteur des lampadaires, la distance des obstacles viennent aussi jouer un rôle dans le calcul de cette importance.

ILLUMINA prend aussi compte des diverses réflexions que peut subir la lumière, de la simple diffusion et de la double diffusion.

Figure 1 :

Img:schéma.jpg Δ

Nous voyons ici comment les différentes diffusions traverse les différentes couches d’atmosphère, est réfléchi vers le sol, pour finalement être capter par un photomètre.

Provenance : Aubé, M. (2010) Heterogeneous modeling of artificial sky luminance. [http://www.cegepsherbrooke.qc.ca/~aubema/index.php/Prof/IllumEn]

Problématique

ILLUMINA prend en considération un grand nombre de facteurs, rallongeant considérablement le temps de calcul, comparativement aux autres modèles non hétérogènes. Un des facteurs importants de ce ralentissement est la double diffusion. En effet, en ne tenant compte que de la simple diffusion, le temps de calcul sur le Mammouth série II est d’environ 12 heures. Or, si nous calculons la double diffusion, nous dépassons les 80 heures (temps maximal accordé par calcul Canada). Le but de notre présent projet est donc d’optimiser le modèle de telle sorte que nous puissions l’utiliser en intégrant sans difficulté la seconde diffusion à l'intérieur d'un temps de calcul de 80 heures pour un domaine géographique couvert par 300x400 pixels.

Cadre méthodologique

Fonctionnement du code

Principe général

Le but du code est de modifier les cartes de luminosité et de réflectance pour diminuer le temps de calcul d’ILLUMINA. Notre stratégie est donc de faire des cartes demandant un moins grand nombre de calcul. Pour ce faire, nous regroupons certains pixels de manière à former un ensemble de zones pour que chacune d’entre elles ne représente qu’un seul calcul pour ILLUMINA. Pour les cartes de luminosité, nous faisons la somme des pixels de chaque zone et la plaçons sur un des pixels de la zone en mettant tout les autres à zéro. La figure 1 montre un exemple d’une carte fictive très petite comportant seulement 4 zones. Les zones sont délimitées par les bordures et chaque chiffre représente la luminosité d’un pixel. Le tableau de gauche est la carte original et celui de droite la carte modifiée qui servira à faire les calculs d’ILLUMINA.

Figure 2 

      Carte de luminosité originale                     Carte de luminosité modifiée

Ainsi, nous devons commencer par choisir un pixel pour chaque zone puis le mettre égal à la somme des pixels sur cette zone. On conserve donc l’importance de la luminosité de chaque zone mais ILLUMINA ne fera qu’un seul calcul pour chacune d’entre elles. En effet, les pixels égaux à zéro ne sont pas pris en compte dans le calcul puisqu’ils ne contribuent pas à la pollution lumineuse.

Pour la carte de réflectance, nous utilisons les mêmes zones, mais en faisant la moyenne pondérée par l’intensité lumineuse de chaque zone et en appliquant cette valeur sur chaque pixel.

La figure 3 illustre la modification d’une carte de réflectance en utilisant les données de la figure 2.

Figure 3 

     Carte de réflectance originale                    Carte de réflectance modifiée

La valeur des pixels de chaque zone est calculée à l’aide de l’équation suivante.

{#R = \frac{\sum_{i} R_i L_i}{\sum_{i} L_i}#}

R est la réflectance moyenne qui sera mis sur tous les pixels d’une zone, Ri est la réflectance d’un pixel et Li la luminosité du même pixel. On somme pour tous les pixels d’une zone.

Afin de diminuer le nombre de calcul tout en gardant une bonne précision des résultats, il est important de garder une haute résolution près du point d’observation. Les pixels restent donc inchangés à l’intérieur d’un certain rayon autour du point d’observation, puis commencent à se regrouper en zones de plus en plus grandes à mesure qu’ils s’en éloignent. La fonction choisie pour calculer les grandeurs des zones est proportionnelle à l’inverse de la distance au point d’observation (1/r) car il s’agit d’un phénomène isotropique à deux dimensions. La figure suivante montre de telles zones. Les teintes sont aléatoires pour permettre de distinguer les zones.

Figure 4 : Un exemple de zones ne dépendant que de la distance au point d'observation

De plus, les régions géographiques possédant une forte émission lumineuse ont plus d’effet sur les résultats calculés par ILLUMINA et il est donc judicieux de garder une meilleure résolution pour ceux-ci. C’est pourquoi la fonction servant à calculer la grandeur des zones est aussi proportionnelle à la luminosité des pixels. L’exemple suivant montre l’effet de la luminosité sur la grandeur des zones. La première carte est la luminosité et la deuxième montre des zones qui ont été calculés avec la luminosité de la carte précédente comme seul paramètre.

Figure 5 : Carte de luminosité du sud du Québec

Figure 6 : Un exemple de zones ne dépendant que de la luminosité de la figure 4

Une fois, que nous avons lancé les calculs avec les nouvelles cartes, ILLUMINA génère deux cartes de résultat. Toutefois, ces cartes ne sont pas représentatives car la luminosité des zones a été placée sur un seul pixel. Nous avons alors créé un code permettant de reconvertir les cartes en ayant une valeur pour tous les pixels. Le principe est d’utiliser les cartes de luminosité originale pour s’avoir, pour chaque pixel, sa proportion de la luminosité totale de sa zone. Ainsi, on utilise ces proportions pour redistribuer la valeur de chaque zone des cartes de résultats. On retrouve donc des cartes d’une précision très semblable à ce qu’on aurait obtenu en faisant tous les calculs.

La figure suivant illustre comment le code s’y prend.

Figure 7 : Carte de résultat, carte de proportion lumineuse et carte de résultat modifié

La première carte est la carte de résultats qu’ILLUMINA a transmit. La deuxième carte à été réalisée avec les données du premier exemple, c’est la proportion de l’intensité lumineuse de chaque pixel par rapport au total de sa zone. La troisième carte a été réalisée en multipliant les valeurs de la carte de résultat de chaque zone avec la valeur de ses pixels correspondant sur la deuxième carte. On obtient donc une carte avec une valeur pour chaque pixel et qui représente bien ce qu’ILLUMINA aurait généré si on avait lancé tous les calculs. Par exemple, un pixel qui n’a aucune luminosité donne une proportion nulle sur la deuxième carte et donc un résultat nul sur la carte finale ce qui correspond bien à la réalité puisqu’il ne contribue pas à la pollution lumineuse. Dans le même ordre d’idée, un pixel fortement lumineux par rapport à sa zone à une proportion plus élevée et un résultat aussi plus élevé ce qui est conforme à la réalité.

Algorithme de détermination des zones

Le code est divisé en deux fichiers, rand_pre.c et rand_post.c, qui sont respectivement pour le traitement avant et après ILLUMINA. Ils se trouvent en annexe.

Le cœur du code est dans la création des zones de regroupement. C’est pourquoi nous expliquons en détail cet algorithme. Le reste du code est assez standard et une explication détaillée de son fonctionnement ne révèlerait pas de faits importants par rapport à la méthodologie de ce projet. La méthode employée est basée sur les probabilités. Nous définissons une fonction qui attribue à chaque pixel une certaine probabilité. Tel que défini plus haut, la fonction diminue selon l’inverse de la distance et augmente proportionnellement à l’intensité lumineuse. La voici :

p = pr2 * ( r2 - r1 ) / ( d * ( 1 - pr2 ) + r2 * pr2 - r1 ) + imp * l / g

p = probabilité attribuée au pixel

d = distance au point d’obervation

r1 = rayon où la fonction de probabilité en 1/r touche 1

r2 = rayon où la fonction de probabilité en 1/r touche pr2

pr2 = probabilité à r2

imp = importance relative de l'intensité lumineuse sur la fonction de probabilité.

l = luminosité du pixel

g = luminosité maximale

Ensuite, un nombre aléatoire entre 0 et 1 est tiré pour chaque pixel. Si ce nombre est en dessous de la probabilité du pixel, le pixel est choisi pour représenter une zone. Nous obtenons une carte ressemblant à la figure suivante où les points blancs représentent les pixels choisis et seront donc le centre d’une zone.

Figure 8 : Carte des pixels représentant une zone (randmap)

La prochaine étape du code est alors de déterminer pour chaque pixel de la carte qu’elle est le pixel choisi le plus près afin de créer les zones. C’est cette opération qui peut être la plus longue de l’algorithme car pour chaque pixel il faut scanner un grand nombre de pixels autour de lui-même afin de trouver le pixel choisi le plus près. La méthode utilisée est de scanner en spirale autour des pixels. Toutefois, puisque dans un univers de pixels carrés les spirales sont aussi de forme carrée, le premier pixel trouvé dans le scan n’est pas nécessairement le pixel le plus près. La figure suivante en montre un exemple.

Figure 9 : Un exemple de spirale pour attribuer une zone à un pixel

Dans cet exemple, le programme devrait effectuer un tour supplémentaire pour trouver le pixel le plus près car ce n’est pas le premier sur la trajectoire mais celui situé au centre gauche. Donc le code doit continuer quelques tours quand il trouve un premier pixel. Toutefois, il doit savoir s’arrêter au bon endroit pour éviter de scanner toute la carte à chaque pixel ce qui prendrait un temps énorme. La figure suivante illustre comment s’assurer de trouver le bon pixel tout en en scannant un minimum.

Figure 10 : Scan complet pour un pixel

Lorsque le scan trouve un premier pixel, on calcul sa distance au point de départ puis on continu la spirale jusqu'à avoir scanné tous les pixels d’un carré de côtés égales au double de la distance calculé précédemment. Ainsi, nous avons calculé tous les pixels à l’intérieur et en périphérie d’un cercle de rayon égal à la distance calculé pour le premier pixel. Nous pouvons donc être sûrs de trouver le pixel le plus près car il est nécessairement à l’intérieur ou en périphérie de ce cercle. Pour chaque pixel choisi que la spiral rencontre, le code calcule sa distance au point de départ et à la toute fin il renvoie les coordonnées du point ayant la plus petite distance.

Nous obtenons ainsi une carte de zones ressemblant à la suivante. (les teintes sont aléatoires)

Figure 11 : Carte de zones complète

Une fois toutes les zones déterminées, nous pouvons commencer à traiter les cartes de luminosité et de réflectance de la façon expliquée plus haut.

Modifications des intrants

Nous avons aussi remarqué que le programme ILLUMINA ne prenait pas en considération la variation de pollution lumineuse entre les régions rurales et les régions urbaines. En effet, les modèles et le nombre de lampadaire par km2 de lampadaire en campagne ne sont pas les mêmes qu’en ville. Nous avons donc décidé de prendre un échantillon des lampadaires dans une région rurale près de la ville de Sherbrooke. Le village que nous avons choisi pour l’échantillonnage est la petite municipalité de Stoke. Nous sommes allés durant le jour pour dénombrer les lampadaires dans la région, en supposant que cet échantillon est représentatif du reste des régions rurales du Québec.

Nous classifions les lampadaires dans deux catégories : Les bons lampadaires (Cobra Head et Hélios des lampadaires qui dirigent une bonne partie de leur lumière vers le bas) et les mauvais lampadaires (Lampes de ferme, qui dirigent beaucoup plus de lumière vers le haut). Après le recensement, on constate effectivement qu’il y a avait une grande différence entre les régions rurales et les régions urbaines au niveau des ratios, des nombres et de la pollution lumineuse des bons et mauvais lampadaires. Nous nous apprêtions à essayer d’intégrer ces nouvelles valeurs dans Illumina mais nous avions oublié qu’il fallait calculer les paramètres photométriques des lampadaires de ferme, puisque ces données ne se trouvent pas dans ILLUMINA puisqu’elles ne possédaient que des valeurs des régions urbaines qui ne contiennent pas de lampadaire de ferme (image ci-dessous).

Pour mesurer les paramètres photométriques du lampadaire de ferme, il a fallu trouver un endroit assez sombre pour éviter toute réflexion de la lumière. Nous avons ensuite utilisé un luxmètre que nous avons placé à une certaine distance du lampadaire que nous avions allumé et nous avons pris en note les mesures de luminosité à tout les 5°. Il aurait été idéal de prendre les mesures à tous les degrés, mais il est possible de crée une bonne courbe de tendance avec seulement les mesures à tous les 5 degré. Nous avons ensuite créé un programme en C appelé cubic_interpolation.c pour interpoler les données et avoir une valeur de la luminosité à n’importe quel angle. Pour chaque segment entre deux points, le code calcule un polynome de degré 3 passant par ces deux points et les deux autres points à gauche et à droite du segment. Ensuite, il calcule la valeur de la luminosité pour chaque degré avec ces polynomes. Le code se trouve en annexe.

En prenant des mesures photométriques d’un lampadaire de ferme qui nous a été fourni par le Cegep de Sherbrooke, nous avons constaté que ces valeurs étaient nettement différentes des valeurs des lampadaires Cobra Head. Cela nous a poussé à nous questionner sur la qualité de notre échantillon et sur ce que nous avions définit comme un lampadaire de ferme. Nous avons fini par prendre la décision de retourner à Stoke pour refaire notre échantillon mais cette fois nous sommes allés durant la nuit. La décision de retourner à Stoke la nuit était dûe au fait que durant la soirée, les lampadaires sont allumés et au lieu de classifier les lampadaires selon les modèles et leur forme, nous pourront les classifier selon leur formes.Les nouveaux résultats que nous avons obtenus étaient très différents des premiers. Cela est en partie dû au fait que notre taille échantillonnale était beaucoup plus grande au 2e dénombrement.

Avec nos nouvelles données et les nouvelles cartes, il était nécessaire de décider des nouvelles zones pour les calculs d’ILLUMINA, sinon les nouvelles données photométriques ne serviraient pas. Il fallait se servir des anciennes cartes et les zones ayant les plus grandes concentrations de lumière feront partie de la nouvelle zone « ville » à laquelle on attribuerait les anciennes données photométriques d’ILLUMINA. Après avoir estimé l’endroit où se trouvaient les grandes zones lumineuses (entre autres Montréal, Québec et Trois-Rivières), il fallait estimer leur superficie par des cercles. Pour calculer ces rayons, il fallait aller sur Google Maps et effectuer une capture d’écran de la région et calculer le rayon à partir de power Point (En dessinant un cercle autour de la région et avec deux rayons perpendiculaires, estimer le centre de la région) En retournant sur Google Maps, on sélectionne le point centre indiqué sur le Power Point sur la carte Google, pour en connaitre les coordonnées géographiques. On effectue les mêmes manipulations pour un point situé sur la circonférence du cercle et ensuite nous nous servions des coordonnées du graphique pour calculer le rayon avec l’aide d'un site Internet trouvé grâce à une recherche sur Google. Le reste de la carte, une fois que l’on exclut les zones « ville », « Sherbrooke » et « Mégantic », sera considérer comme étant la zone « campagne » à laquelle on demandera à ILLUMINA d’utiliser nos nouvelles mesures photométriques.

Matériel, instrumentation et expérimentation

  1. Modèle de la pollution lumineuse illumina
  2. Ordinateur Mammouth série II
  3. Logiciel Code::Blocks
  4. Logiciel GIMP
  5. google maps
  6. Powerpoint
  7. Données d'entrée
    • Élévation du sol (srtm (Shuttle Radar Topography Mission))
    • Luminosité émise vers le haut, la nuit (dmsp-ols)
    • Réflectance du sol dans le visible (modis)
    • Fonction photométrique des lampadaires (lumière en fonction de l'angle)
    • Hauteur typique des lampadaires (9 m)
    • Hauteur typique des obstacles (bâtiments et arbres) (9,2 m)
    • Distance moyenne entre les obstacles

Test du code

Un des premiers tests effectués pour vérifier si le code fonctionne correctement a été de le tester sans passer par ILLUMINA. Autrement dit, nous avons ajusté la résolution des cartes d'entrées (luminosité et réflectance) puis avons mis les nouvelles cartes de luminosité directement dans le programme rand_post.c (reconversion des cartes après ILLUMINA). Le résultat a été de retrouver exactement les cartes de luminosité originales. C'est ce que nous nous attendions et c'est signe que la reconversion fontionne très bien. En effet, le but est de retrouver une carte le plus précisément possible identique à ce qu'on aurait eu en entrant les cartes en résolution complète. Or, si aucun changement n'est fait dans ILLUMINA, nous ne devons pas avoir aucun changement dans nos cartes aussi, mêmes si elles passent par une étape où leur résolution est changée.

L'autre important test a été de lancer des calculs en simple diffusion avec différents paramètres du code d'ajustement de la résolution des cartes entrantes. En effet, les calculs en simple diffusion sont relativement courts à effectuer et il est alors possible de comparer différents ajustement avec les résultats en résolution complète. Nous avons donc lancé le même calcul plusieurs fois en ne changeant que les paramètres de rand_pre.c. Un de ces tests était en résolution complète. Pour les autres, nous avons testé différentes versions de la fonction en 1/r (en faisant varier r2 et pr2) tout en gardant la contribution à la luminosité constante ( imp = 0.7 ).

Observations, interprétations et résultats

Pour analyser les résultats, nous avons pris une expérience de première diffusion, car nous ne possédons aucun résultat de double diffusion avec 100% des calculs. Pour ce faire nous avons utilisé l'expérience été avec une longueur d'onde de 568 nm. Nous avons donc inséré différents paramètres dans le programme rand_pre pour voir la précision et le temps de calcul. Nous avons pour l'instant seulement analysé l'impact du PR2 sur les résultats. Le PR2 représente la probabilité de d'avoir un calcul à faire lorsque nous nous trouvons au rayon 2. Nous avons donc gardé plusieurs paramètres constant.

Nous avons pris une expérience avec point d'observation le Mont-Mégantic. Nous avons gardé le rayon 1 constant à 30km. Le rayon 1 représente le rayon autour du point d'observation ou nous effectuons la totalité des calculs. Nous avons aussi gardé le rayon 2 constant à 200km. Celui ci est le rayon autour du point d'observation ou nous faisons le moins de calculs. Nous avons aussi gardé l'importance accordée à la luminosité des villes constante à 0,7. L'importance de la luminosité des villes est utile pour augmenter le nombre de calcul dans les grandes villes comme Montréal, même si elles sont éloignées du point d'observation.

Nous avons donc testé 4 PR2 différents. Les PR2 que nous avons utilisés sont les suivants: 1; 0,1; 0,01; 0,001.

Note: Pour voir plus facilement les 4 images de contribution lumineuse ci-dessous, nous avons appliqué une normalisation et 2 corrections gamma de 1,6 chacune. Cette correction se retrouve uniquement dans ce rapport et ne fait pas part des étapes importantes notre analyse.

Premier Test

Pour le premier test nous avons utilisé un PR2 = 1. Avec un PR2 de 1 nous effectuons la totalité des calculs. Le résultat que nous avons obtenu est 4,258E-07(W/str/m**2). Les calculs ont duré 12h 37 min et 37s et pour nous cette résolution est considéré de 100%.

Deuxième test

Pour le deuxième test nous avons utilisé un PR2 = 0,1. Un PR2 de 0,1 signifie que au rayon 2, 200km, nous effectuons 1 calcul sur 10 plutôt de la totalité des calculs. Le résultat que nous avons obtenu est 4,247E-07(W/str/m**2). Les calculs ont duré 2h 26 min et 21 sec. Soit 19,31% du temps du calcul pleine résolution avec un résultat précis à 99,74% du premier résultat.

Troisième test

Pour le troisième test nous avons utilisé un PR2 = 0,01. Un PR2 de 0,01 signifie que au rayon 2, 200km, nous effectuons 1 calcul sur 100 plutôt que la totalité des calculs. Le résultat que nous avons obtenu est 4,291E-07(W/str/m**2). Les calculs ont duré 50 min et 21 sec. Soit 6,65 % du temps de calcul pleine résolution avec un résultat précis à 99,23% du premier résultat.

Quatrième test

Pour le quatrième test nous avons utilisé un PR2 = 0,001. Un PR2 de 0,001 signifie que au rayon 2, 200km, nous effectuons 1 calcul sur 1000 plutôt que la totalité des calculs. Le résultat que nous avons obtenu est 4,385E-07(W/str/m**2). Les calculs ont duré 38 min et 34 sec. 5,09 % du temps de calcul pleine résolution avec un résultat précis à 97,10% du premier résultat. Aussi, grâce à la modification des cartes, nous avons accès à la deuxième diffusion. On ne peut pas dire la précision exacte, mais nous obtenons le résultat de 5,251E-07(W/str/m**2). Donc l'a contribution de la deuxième diffusion dans ce cas de figure est de l'ordre de 19,75% de plus de la simple diffusion.

Voici ici une petite animation des 4 images de contribution lumineuse de chacun des tests.

======= =======

>>>>>>> Nous avons observé que la modifications des zones des cartes et des paramètres photométriques étaient très justifiés. En effet, nous avons remarquer que le rapport de lampadaire de ferme sur nombre total de lampadaire était de 2/3. Nos nouvelles valeurs photométriques se retrouve dans ce graphique:

Nous avons aussi réunis les cordonnées et les rayons des cercles des différentes ville qui compose la nousvelle zone "ville" d'ILLUMINA dans le tableau ci-dessous:

Villecoordonné en x (en pixel)coordonné en y (en pixel)rayon (en pixel)
Drummondville1981056
Granby180523
Magog225372
Victoriaville2401254
Montréal1066720
Trois-Rivières1921595
Québec29220913
Sorel1491233
Shawinigan1791794
Grand-mère1831873
Thetford Mines2891284
St-Georges2641323
St-Agathe-des-Monts611232
St-Sauveur661074
St-Jerôme79944
Blainville957612

Discussion

Comme nous l’avons vu dans les résutats de radiance spectrale modélisée, moins il y a de calcul plus le résultat est petit. Cependant, quand nous regardons les images de contribution lumineuse, les tests avec moins de calculs ont des taches supplémentaires. Au début cela nous semblait étrange, mais nous avons vite pensé à une raison possible qui expliquerait pourquoi les cartes avaient ces taches.Pour expliquer ce phénomène nous devons bien comprendre ces images. Elles représentent le pourcentage de contribution de chaque pixel au total de radiance spectrale modélisée. Puisqu'il y a un grand nombre de pixels, le poids de chaque pixel est faible. Puisque le résultat est plus petit, le poids de chaque pixel dans le total augmente légèrement. Cette petite augmentation est suffisante pour faire une différence parce que leur poids original est si petit. Cet effet est ensuite amplifié par la normalisation et par la correction gamma que nous effectuons sur les cartes.

Il nous est possible de tracer un graphique pour mettre en relation le temps de calcul en fonction du pourcentage d'erreur. Dans ce graphique nous pouvons voir qu'il y a une zone qui ce situe entre le point PR2 = 0,1 et PR2 = 0,01 ou le pourcentage d'erreur est approximativement 0,5% et que le temps de calcul est d'environ 1 heure.

On peut alors mettre en évidence dans le tableau 1 l'erreur divisé par le temps de calcul.

Tableau 1 Erreur/temps de calcul en fonction du PR2
PR2 Erreur/temps de calcul (%/h)
1,000 0
0,100 0,11
0,010 0,83
0,001 4,60

Il est donc possible d'avoir plus de précision si nous avons beaucoup de temps, mais on peut aussi avoir un résultat précis en diminuant le temps de calcul.

Dans l'expérience que nous avons analysé, nous avons seulement eu un résultat de double diffusion lorsque nous avions un PR2 = 0,001. Les calculs de double diffusion avec les PR2 plus grand dépassait le temps de calcul de 80 heures qui nous était alloué au Mammouth. Nous avons calculé avec la première diffusion de lorsque le PR2 = 0,001 le pourcentage d'erreur est 2,90%. De plus, avec la double diffusion notre pourcentage d'erreur est probablement un peu plus haut, car il y a beaucoup plus de calcul. Cependant, nous avons aussi calculé que l'apport de la deuxième diffusion est de 19,75% au total de la radiance spectrale modélisée. Dans cette situation, le pourcentage d'erreur représente au moins 14,68% de l'apport de la deuxième diffusion. On peut alors se faire une idée représentative de la double diffusion, mais pas suffisamment précise pour le moment.

Donc, si nous voulons avoir un pourcentage d'erreur plus petit que 10% de la contribution de la double diffusion, nous devons obtenir un pourcentage d'erreur inférieur à 1,975%. Revenons à notre graphique de la précision en fonction du temps de calcul. On remarque que le 1,975% d'erreur se trouve approximativement entre le PR2 = 0,01 et PR2 = 0,001. On aurait alors la précision minimal avec un PR2 de 0,0055 environ.

Actuellement, nous ignorons si un calcul de double diffusion peut être effectué avec un PR2 = 0,0055.

Conclusion

En conclusion, nous pourrions dire que notre projet fut un franc succès. Notre but était d'optimiser le temps de calcul d'ILLUMINA sans pour autant trop empiétrer sur la précision des valeurs. Nous avons en effet réussis à créer un programme qui, en modifiant les données intrantes, pouvait réduire le temps de calcul de 95%, tout en ne perdant que 3% de précision. Ce gain de temps nous à permis de pouvoir calculer la double diffusion et de constater qu'elle augmentait la quantité de pollution lumineuse observé de près de 20%, ce qui est très significatif. Cependant, au niveau de la deuxième diffusion, il reste toujours du travail à accomplir, nous en parlerons un peu plus loin. Il est aussi à remarquer que le programme qui modifie les cartes est fonctionnelle et intégrer à Illumina. Le travail maintenant est de tester des combinaisons pour trouver la plus avantageuse. L'intégration de ce programme ne faisait pas parti de nos but au départ, mais il représente tout de même un avancement dans le projet.

Dans le futur, nous devrions tester un test de double diffusion avec les mêmes paramètres que ceux que nous avons analysé et un PR2 = 0,0055. Avec ce résultat, si nous l'obtenons, nous aurons une valeur de pollution avec les deux diffusions et une précision acceptable. Que nous obtenons le résultat ou pas, nous devrions par la suite essayer de nouveaux paramètres pour se faire une idée des paramètres optimaux. Par exemple, un rayon 1 de 30km n'est peut-être pas nécessaire. On pourrait commencer la diminution de calculs immédiatement. Aussi, on pourrait diminuer le PR2 et augmenter le rayon 2. De cette façon, nous aurions le même nombre de calcul à 200km du Mont-Mégantic, mais nous aurions encore moins de calcul dans les régions à plus de 200km. Une autre piste intéressante serait de jouer avec l'importance des villes. Pour notre expérience, nous avions une importance de 0,7. En augmentant ce chiffre, nous pourrions densifier les calculs près des grandes villes, la ou il y a beaucoup de lumière, mais moins de calculs dans les bois.

L'autre volet de notre projet aussi fut ateint. Nous avons trouver que la proportion de lanterne de ferme dans les zones rurales était significative (2/3), et qu'elle ont un comportement qui diffère des cobrahead en ce qui à attret à la pollution lumineuse.

Nous pourrions donc, dans l'avenir, faire une autre modélisation de la pollution lumineuse pour l'OMM, mais cette fois, avec notre programme d'optimisation permettant de calculer la deuxième diffusion et la mise à jour des intrants.

De plus, ce qui pourrait également être amélioré dans le futur, c'est la manière de définir les zones (surtout les zones rurales des zones villes). En effet, nous y sommes allés de manière arbitraire, en choisissant les zones qui nous semblaient les plus lumineuses. Toutefois, il serait intéressant de faire un programme qui, comme ILLUMINA, aurait un aspect hétérogène et saurait délimiter les zones rurales/villes en fonction de la luminosité. Pour ce faire, nous pourrions utiliser un plafond de luminosité pour déterminer les zones de ville. De plus, pour déterminer la valeur de ce plafond, nous pourrions utiliser la moyenne des valeurs des pixels de luminosité, où leur dérivé seconde s’inverse.

Médiagraphie

Aubé, M. (2010) Heterogeneous modeling of artificial sky luminance. [http://www.cegepsherbrooke.qc.ca/~aubema/index.php/Prof/IllumEn]

Aubé, M. (2007)Light pollution modeling and detection in a heterogenous environment [http://cegepsherbrooke.qc.ca/~aubema/index.php/Prof/Page?action=download&upname=Light-Pollution-Modeling-Aube2007.pdf]

Aubé, M., Franchomme­-Fossé, L.,Robert­-Staehler P. ,Houle V. (2005) Light Pollution Modelling and Detection in a Heterogeneous Environment: Toward a Night Time Aerosol Optical Depth Retrieval Method. [http://cegepsherbrooke.qc.ca/~aubema/index.php/Prof/Page?action=download&upname=article-san-diego-corrige.pdf]

D'arcy-Lepage, A., Fontaine, M., Pfister, C. (2011) Contrôle de la pollution lumineuse à l'observatoire du mont Mégantic (OMM) [http://www.cegepsherbrooke.qc.ca/~aubema/index.php/A10Cours01GroupeAEquipe01/RapportGraphycs1]

Kocifaj M., Light-pollution model for cloudy and cloudless night skies with ground-based light sources [http://www.opticsinfobase.org/abstract.cfm?URI=ao-46-15-3013]

Annexe

Données photométriques d'un lampadaire de ferme

Intensité Lumineuse (lux)Angle (degré)
0.00006460795890536230
0.00007156573909517051
0.00007852351928497882
0.0000854812994747873
0.00009243907966459524
0.00009939685985440355
0.00009995348226958826
0.00009967517106199587
0.00009911854864681118
0.00009884023743921889
0.000099396859854403510
0.00010373056294405511
0.00010941606332772712
0.00011641360226147713
0.00012468342100136414
0.00013418576080344515
0.00014201823335997216
0.00015171936688176117
0.00016396506001582418
0.00017943121140916919
0.00019879371970880720
0.00023020312742279821
0.0002649920283718422
0.00030196766023767823
0.0003399372607020624
0.00037770806744673325
0.00041042951371080326
0.00044148109272931827
0.00047058449329468828
0.00049746140419931929
0.00052183351423561830
0.00054167312746255731
0.0005588886635893432
0.00057363915759173333
0.00058608364444550534
0.00059638115912642135
0.00060596301641638536
0.00061339790153349537
0.00061852677950198238
0.0006211906153460839
0.00062123037409002240
0.00061276176163042741
0.00060278231690104442
0.00059256436940644243
0.00058338004985746544
0.0005765017871555445
0.00058338004985746546
0.00059256436940644247
0.00060278231690104448
0.00061276176163042749
0.00062123037409002250
0.00061403404143656351
0.00060600277516032752
0.00059908475371446153
0.0005952281555521154
0.00059638115912642155
0.00062103158037031356
0.00065045305088721657
0.00079301064994908558
0.00093596593210257659
0.00097064543887153160
0.00097987929065486361
0.00099311888789284662
0.0010087440244830663
0.0010342582247968364
0.0010685193634935165
0.0011247268723632166
0.0011928620044154467
0.0012731822627462968
0.0013654301571687869
0.0014633882629313270
0.0015732740802364671
0.0017015623398790572
0.0018417782907972473
0.0020047195125877974
0.002179588497491775
0.0023822363923201276
0.0025944263230870277
0.0028203582246129978
0.0030510612251104179
0.0032608252670644780
0.0034610471110756581
0.0036734277484477682
0.0038858085975227183
0.0041159401724680284
0.0043460721082114685
0.0045989974701891786
0.0048535135407846487
0.0050897809276829288
0.0053308194671271889
0.0055413190229079690
0.0057883967702724191
0.0060539002698387792
0.0063162231951048993
0.0065494728028260294
0.0067684094008523495
0.0070272518489651496
0.007266613119898797
0.0075694990492461498
0.0078536995023145999
0.0081158397694364100
0.00834259474347897101
0.00856811276449167102
0.00878528147187596103
0.00935607653035246104
0.00992806454672207105
0.0105883929061436106
0.0112523014965391107
0.0114550131954645108
0.0116565305611024109
0.0115716353437232110
0.0114744156018098111
0.0114891941381309112
0.0114964170118364113
0.0113631051367781114
0.0112206499833837115
0.0111743794863117116
0.0111189635784044117
0.0111512325125819118
0.0111759481851392119
0.0111571728531997120
0.0111324343260543121
0.0111498383742055122
0.0111628702229878123
0.0112331339589206124
0.0113006146324758125
0.0114397021280529126
0.0115811757149534127
0.0117023869577548128
0.0118347299130737129
0.0118360107046686130
0.0118889775748875131
0.0120129137986335132
0.012149573784715133
0.0122635565530802134
0.0123831054958988135
0.0124163109110588136
0.0124522992046498137
0.0123842347760456138
0.012320543466543139
0.0123059145118885140
0.0123338257131913141
0.0123125923059982142
0.012280624306437143
0.0121861973145552144
0.0120643371075483145
0.011927928098349146
0.0117644835311666147
0.0115695806712807148
0.0113651341915779149
0.0112134191635169150
0.0111205484783595151
0.0110055517515814152
0.0108905550248033153
0.0107297907718065154
0.0105610742033346155
0.0103286057942332156
0.0100865938150133157
0.00982320697907165158
0.0095550494914444159
0.00940434756792757160
0.00926636886909665161
0.00907842728571053162
0.00888889512085982163
0.00869502647941861164
0.00849797748441302165
0.0083359806097006166
0.00816682728047962167
0.00801450640938429168
0.0078486668476286169
0.00768943538287311170
0.00745784295444517171
0.00722122184114136172
0.00698460132228066173
0.00675578608916676174
0.00654685067530961175
0.0065251289248664176
0.00650340694077821177
0.00650107257514637178
0.00649873837345699179
0.00643453674259696180

rand_pre.c

#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#define lh 6                /*Number of lines in the headers of the lumlp files before the line of the width, height and gain*/
#define seed 2              /*Seed of the random number generator*/
#define use_default_data 0  /*If not equal to zero, it will use default parameters without asking for it*/

void make_areas(char* nlumin)
{
    FILE *f=NULL;
    int ox,oy,x,y,z,w,h,g,i,j,it,jt,itt=0,jtt=0,pi,pj,t,n,k,lim,s,u;
    float r1,r2,pr2,rm,imp,d,p;
    char ch,name[50];

    srand(seed);

/*Acquiring the values of the parameters for making the areas*/
    if (use_default_data==0)
    {
    printf("\nox = ");          /*x coordinate of the observation point starting with one and from the left*/
        scanf("%d",&ox);
    printf("\noy = ");          /*y coordinate of the observation point starting with one and from the bottom*/
        scanf("%d",&oy);
    printf("\nr1 = ");          /*Radius that determine where the probability function (witch is proportional to 1/r) is equal to 1*/
        scanf("%f",&r1);
    printf("\nr2 = ");          /*Radius that determine where the probability function (witch is proportional to 1/r) is equal to pr2*/
        scanf("%f",&r2);
    printf("\npr2 = ");         /*Value of the probability function at r2*/
        scanf("%f",&pr2);
    printf("\nrm = ");          /*Radius under witch the probability is reajusted to 1 (facultative, put 0 to cancel its effect)*/
        scanf("%f",&rm);
    printf("\nimp = ");         /*Relative importance of the luminosity of each pixel to the probability of being choosed to represent an area*/
        scanf("%f",&imp);
    }

/*Default parameters. Will be active only if use_default_data is not defined as 0*/
    else
    {
        ox=302;
        oy=58;
        r1=1;
        r2=0;
        pr2=0;
        rm=0;
        imp=0;
    }

    if(rm<r1)rm=r1;

/*Open the file containing the luminosity of all the zones.*/
    sprintf(name,"%s_recombined.pgm",nlumin);
    f=fopen(name,"r");
/*Skip the headers*/
    while ((ch=getc(f))!='\n');
/*Read the width, height and gain*/
    fscanf(f,"%d %d %d",&w,&h,&g);
/*Create the required tables*/
    char m[h][w];   /*Will contain the values of the randmap. That is a table of 0 and 1 dermining witch pixel will be calculated in ILLUMINA.*/
    unsigned short int mi[h][w],mj[h][w];   /*Will contain the information of what pixel in m[][] each pixel will be merged. In other words, it is the information of the areas.*/
    unsigned short int m1[h][w];    /*Will contain the luminosity of each pixel*/
/*Transform the cartesian coordinates of the observation point in table coordinates*/
    pi=h-oy;
    pj=ox-1;
/*Put the values of the luminosity file in m1[][]*/
    for(i=0;i<h;i++)
        for(j=0;j<w;j++)
        {
            fscanf(f,"%d",&x);
            m1[i][j]=x;
        }
    fclose(f);
/*Make the randmap*/
    for(i=0;i<h;i++)
        for(j=0;j<w;j++)
        {
            d=sqrt((pi-i)*(pi-i)+(pj-j)*(pj-j));    /*Distance to the observator*/
            if (r2==0) p=r1/d+imp*(float)m1[i][j]/g;    /*Simplified 1/r function that pass through 1 at r1 and as a vertical asymptote at r=0*/
            else p=pr2*(r2-r1)/(d*(1-pr2)+r2*pr2-r1)+imp*(float)m1[i][j]/g; /*Main function. Respect all the parameters*/
            if (((float)(rand()%32769)/32768)<p || p<0 || d<=rm ) m[i][j]=1;    /*Put 1 in m[][] if a random float number between 0 and 1 is under the probability function.*/
            else m[i][j]=0;
        }
/*Save m[][] in a file*/
    f=fopen("randmap.pgm","w");
    fprintf(f,"P2\n%d %d %d\n",w,h,1);
    for(i=0;i<h;i++)
        for(j=0;j<w;j++)
            fprintf(f,"%d\n",m[i][j]);
    fclose(f);

/*Calculate the areas by determining for each pixel what is the nearest pixel equal to 1 in m[][]. Each pixel equal to one in m[][] will therefore be the center of an area containing all the pixels near to it.*/
    printf("\nGENERATING THE AREAS :\n");

    for(i=0;i<h;i++)
    {
        for(j=0;j<w;j++)
        {
            for(it=i,jt=j,n=1,z=(w+h)*(w+h),lim=w+h;n<=lim;n++)
                for(s=0;s<=1;s++)
                    for(k=0;k<n;k++)
                    {
                        if( it>=0 && it<h && jt>=0 && jt<w )
                            if(m[it][jt]==1)
                            {
                                y=(i-it)*(i-it)+(j-jt)*(j-jt);
                                if (y<z)
                                {
                                    if(lim==w+h)lim=2*sqrt(y)+1;
                                    z=y;
                                    itt=it;
                                    jtt=jt;
                                }
                            }
                        if(s==0)it+=pow(-1,n);
                        if(s==1)jt+=pow(-1,n);
                    }
            mi[i][j]=itt;   /*Put the information of the areas in mi and mj.*/
            mj[i][j]=jtt;
        }
        if (i%(h/100)==0) printf("%.0f%%\r",(float)i/h*100);    /*Print the percentage of the calculation that is done.*/
    }
    printf("%.0f%%\n\n",(float)i/h*100);

/*Save mi and mj in files.*/
    sprintf(name,"%s_mi.rand",nlumin);
    printf("Create %s\n",name);
    f=fopen(name,"w");
    for(i=0;i<h;i++)
        for(j=0;j<w;j++)
            fprintf(f,"%d\n",mi[i][j]);
    fclose(f);
    sprintf(name,"%s_mj.rand",nlumin);
    printf("Create %s\n",name);
    f=fopen(name,"w");
    for(i=0;i<h;i++)
        for(j=0;j<w;j++)
            fprintf(f,"%d\n",mj[i][j]);
    fclose(f);
}

void make_lumlp_out(char* nlumin)
{
    FILE *f=NULL,*fn=NULL,*fmi=NULL,*fmj=NULL;

    int mi,mj,w,h,g,x,i,j,nzone;
    char name[50],ch;

    for(nzone=1;nzone<=9;nzone++)
    {
        sprintf(name,"%s_0%d.pgm",nlumin,nzone);
        f=fopen(name,"r");

        if (f==0) break;

        sprintf(name,"%s_mi.rand",nlumin);
        fmi=fopen(name,"r");
        sprintf(name,"%s_mj.rand",nlumin);
        fmj=fopen(name,"r");
        sprintf(name,"%s_0%d_new.pgm",nlumin,nzone);
        fn=fopen(name,"w");

        printf("Create %s\n",name);

        for(x=0;x<lh;x++){
            while ((ch=getc(f))!='\n') fprintf(fn,"%c",ch);
                fprintf(fn,"\n");}

        fscanf(f,"%d %d %d",&w,&h,&g);

        int m[h][w];

        for(i=0;i<h;i++)
            for(j=0;j<w;j++)
                m[i][j]=0;

        for(i=0;i<h;i++)
            for(j=0;j<w;j++)
            {
                fscanf(fmi,"%d",&mi);
                fscanf(fmj,"%d",&mj);
                fscanf(f,"%d",&x);
                m[mi][mj]+=x;
            }

        fclose(fmi);
        fclose(fmj);
        fclose(f);

        for(i=g=0;i<h;i++)
            for(j=0;j<w;j++)
                if(m[i][j]>g) g=m[i][j];

        fprintf(fn,"%d %d %d\n",w,h,g);

        for(i=0;i<h;i++)
            for(j=0;j<w;j++)
                fprintf(fn,"%d\n",m[i][j]);

        fclose(fn);
    }

    sprintf(name,"%s_recombined.pgm",nlumin);
    f=fopen(name,"r");
    sprintf(name,"%s_mi.rand",nlumin);
    fmi=fopen(name,"r");
    sprintf(name,"%s_mj.rand",nlumin);
    fmj=fopen(name,"r");
    sprintf(name,"%s_recombined_new.pgm",nlumin);
    fn=fopen(name,"w");

    printf("Create %s\n",name);

    while ((ch=getc(f))!='\n') fprintf(fn,"%c",ch);
        fprintf(fn,"\n");

    fscanf(f,"%d %d %d",&w,&h,&g);

    int m[h][w];

    for(i=0;i<h;i++)
        for(j=0;j<w;j++)
            m[i][j]=0;

    for(i=0;i<h;i++)
        for(j=0;j<w;j++)
        {
            fscanf(fmi,"%d",&mi);
            fscanf(fmj,"%d",&mj);
            fscanf(f,"%d",&x);
            m[mi][mj]+=x;
        }

    fclose(fmi);
    fclose(fmj);
    fclose(f);

    for(i=g=0;i<h;i++)
        for(j=0;j<w;j++)
            if(m[i][j]>g) g=m[i][j];

    fprintf(fn,"%d %d %d\n",w,h,g);

    for(i=0;i<h;i++)
        for(j=0;j<w;j++)
            fprintf(fn,"%d\n",m[i][j]);

    fclose(fn);
}

void reflec(char* nlumin,char* nrefin)
{
    FILE *ffull=NULL,*fnew=NULL,*fp=NULL,*fref=NULL,*frefnew=NULL,*fmi=NULL,*fmj=NULL;
    int x,i,j,mi,mj,w,h,g,r;
    float p,pt;
    char name[50],ch;

    sprintf(name,"%s_mi.rand",nlumin);
    fmi=fopen(name,"r");
    sprintf(name,"%s_mj.rand",nlumin);
    fmj=fopen(name,"r");
    sprintf(name,"%s_recombined.pgm",nlumin);
    ffull=fopen(name,"r");
    sprintf(name,"%s_recombined_new.pgm",nlumin);
    fnew=fopen(name,"r");
    sprintf(name,"%s_p.rand",nlumin);
    printf("Create %s\n",name);
    fp=fopen(name,"w");
    printf("Create %s\n",name);
    sprintf(name,"%s_new.pgm",nrefin);
    frefnew=fopen(name,"w");
    printf("Create %s\n",name);

    while ((ch=getc(ffull))!='\n');
    fscanf(ffull,"%d %d %d",&w,&h,&g);

    while ((ch=getc(fnew))!='\n');
    fscanf(fnew,"%d %d %d",&w,&h,&g);

    int m[h][w];

    for(i=0;i<h;i++)
        for(j=0;j<w;j++)
            fscanf(fnew,"%d",&m[i][j]);

    fclose(fnew);

    for(i=0;i<h;i++)
        for(j=0;j<w;j++)
        {
            fscanf(fmi,"%d",&mi);
            fscanf(fmj,"%d",&mj);
            fscanf(ffull,"%d",&x);
            if (m[mi][mj]==0) p=0;
            else p=(float)x/m[mi][mj];
            fprintf(fp,"%f\n",p);
        }

    fclose(fmi);
    fclose(fmj);
    fclose(ffull);
    fclose(fp);

    sprintf(name,"%s.pgm",nrefin);
    printf("\nRead %s :\n\n",name);
    fref=fopen(name,"r");

    for(x=0;x<lh+1;x++){
        while ((ch=getc(fref))!='\n') {printf("%c",ch); fprintf(frefnew,"%c",ch);}
        printf("\n"); fprintf(frefnew,"%c",ch);}

    fscanf(fref,"%d %d %d",&w,&h,&g);
    printf("%d %d %d\n",w,h,g);
    fprintf(frefnew,"%d %d %d\n",w,h,g);

    float ref[h][w];

    sprintf(name,"%s_p.rand",nlumin);
    fp=fopen(name,"r");

    for(i=0;i<h;i++)
        for(j=0;j<w;j++)
            {
                fscanf(fp,"%f",&p);
                fscanf(fref,"%d",&r);
                ref[i][j]=p*r;
            }

    fclose(fref);
    fclose(fp);

    sprintf(name,"%s_mi.rand",nlumin);
    fmi=fopen(name,"r");
    sprintf(name,"%s_mj.rand",nlumin);
    fmj=fopen(name,"r");

    for(i=0;i<h;i++)
        for(j=0;j<w;j++)
            {
                fscanf(fmi,"%d",&mi);
                fscanf(fmj,"%d",&mj);
                ref[mi][mj]+=ref[i][j];
            }

    fclose(fmi);
    fclose(fmj);

    sprintf(name,"%s_mi.rand",nlumin);
    fmi=fopen(name,"r");
    sprintf(name,"%s_mj.rand",nlumin);
    fmj=fopen(name,"r");

    for(i=0;i<h;i++)
        for(j=0;j<w;j++)
            {
                fscanf(fmi,"%d",&mi);
                fscanf(fmj,"%d",&mj);
                ref[i][j]=ref[mi][mj];
            }

    fclose(fmi);
    fclose(fmj);

    for(i=0,pt=0;i<h;i++)
        for(j=0;j<w;j++){
            p=ref[i][j];
            if (p>pt) pt=p;}
    g=(int)(pt+1);

    for(i=0;i<h;i++)
        for(j=0;j<w;j++)
            fprintf(frefnew,"%d\n",(int)(ref[i][j]));

    fclose(frefnew);
}

void recombine(char* nlumin)
{
    FILE *f=NULL;

    int x,i,j,nz,v,w,h,g;
    char nend[50],name[50],ch;

    sprintf(nend,"_0%d.pgm",1);

    sprintf(name,"%s%s",nlumin,nend);

    f=fopen(name,"r");

    if (f==0) {printf("\nERROR : %s not found",name); getchar();}

    for(x=0;x<lh;x++)while((ch=getc(f))!='\n');

    fscanf(f,"%d %d %d",&w,&h,&g);

    fclose(f);

    int m[h][w];

    for(i=0;i<h;i++)
        for(j=0;j<w;j++)
            m[i][j]=0;

    for(nz=1;nz<=9;nz++)
    {
        sprintf(nend,"_0%d.pgm",nz);

        sprintf(name,"%s%s",nlumin,nend);

        f=fopen(name,"r");

        if (f==0) break;

        printf("\nRead %s :\n\n",name);

        for(x=0;x<lh;x++){
            while ((ch=getc(f))!='\n') printf("%c",ch);
                printf("\n");}

        fscanf(f,"%d %d %d",&w,&h,&g);
        printf("%d %d %d\n",w,h,g);

        for(i=0;i<h;i++)
            for(j=0;j<w;j++)
            {
                fscanf(f,"%d",&v);
                m[i][j]+=v;
            }

        fclose(f);
    }

    sprintf(name,"%s_recombined.pgm",nlumin);
    printf("\nCreate %s\n",name);
    f=fopen(name,"w");

    fprintf(f,"P2\n%d %d %d\n",w,h,g);

    for(i=0;i<h;i++)
        for(j=0;j<w;j++)
            fprintf(f,"%d\n",m[i][j]);

    fclose(f);
}

int main()
{
    int i;
    char nlumin[50],nrefin[50];

/*Acquiring the name of the luminosity and reflectance map that will be modified by this program*/
    printf("lumlp input : ");
        for(i=0;(nlumin[i]=getchar())!='\n';i++); nlumin[i]='\0';
    printf("reflect input : ");
        for(i=0;(nrefin[i]=getchar())!='\n';i++); nrefin[i]='\0';

/*Recombine all the zones into one file*/
    recombine(nlumin);

/*Determine the areas*/
    make_areas(nlumin);

/*Make the new luminosity maps*/
    make_lumlp_out(nlumin);

/*Make the new reflectance map*/
    reflec(nlumin,nrefin);
    return 0;
}
 
[$[Get Code]]

rand_post.c

#include<stdio.h>
#define lh 6

main()
{
    FILE *data,*mi,*mj,*f,*fw;
    int x,y,w,h,g,i,j,im,jm;
    float p;
    char ch,nlumpl[50],in[50],name[50];

    printf("original lumlp file : ");
        for(i=0;(nlumpl[i]=getchar())!='\n';i++); nlumpl[i]=0;

    printf("\ninput file : ");
        for(i=0;(in[i]=getchar())!='\n';i++); in[i]=0;

    sprintf(name,"%s_p.rand",nlumpl);

    printf("\n%s\n",name);

    data=fopen(name,"r");

    sprintf(name,"%s_mi.rand",nlumpl);

    printf("%s\n",name);

    mi=fopen(name,"r");

    sprintf(name,"%s_mj.rand",nlumpl);

    printf("%s\n",name);

    mj=fopen(name,"r");

    sprintf(name,"%s.pgm",in);
    f=fopen(name,"r");
    printf("\n%s",in);
    printf("\n\n");

    sprintf(name,"%s_new.pgm",in);
    fw=fopen(name,"w");

    for(x=0;x<lh;x++){
        while ((ch=getc(f))!='\n') {printf("%c",ch);fprintf(fw,"%c",ch);}
        {printf("\n");fprintf(fw,"\n");}}

    fscanf(f,"%d %d %d\n",&w,&h,&g);
    printf("width = %d\nheight = %d\ngain = %d\n",w,h,g);

    unsigned int m1[h][w],m2[h][w];

    for(i=0;i<h;i++)
        for(j=0;j<w;j++)
            m1[i][j]=0;

    for(i=0,x=0;i<h;i++)
        for(j=0;j<w;j++)
        {
            fscanf(f,"%d",&x);
            fscanf(mi,"%d",&im);
            fscanf(mj,"%d",&jm);
            m1[im][jm]+=x;
        }

    fclose(f);
    fclose(mi);
    fclose(mj);

    sprintf(name,"%s_mi.rand",nlumpl);
    mi=fopen(name,"r");
    sprintf(name,"%s_mj.rand",nlumpl);
    mj=fopen(name,"r");

    for(i=0;i<h;i++)
        for(j=0;j<w;j++)
        {
            fscanf(mi,"%d",&im);
            fscanf(mj,"%d",&jm);
            fscanf(data,"%f",&p);
            m2[i][j]=(float)p*(float)m1[im][jm]+0.5;
        }

    fclose(data);
    fclose(mi);
    fclose(mj);

    for(i=g=0;i<h;i++)
        for(j=0;j<w;j++)
            if (m2[i][j]>g) g=m2[i][j];

    fprintf(fw,"%d %d %d\n",w,h,g);

    for(i=0;i<h;i++)
        for(j=0;j<w;j++)
            fprintf(fw,"%d\n",m2[i][j]);

    fclose(fw);
}
 
[$[Get Code]]

cubic_interpolation.c

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

void printm(double **m,int w,int h)
{
    int i,j;
    printf("\n");
    for(i=0;i<h;i++){
        for(j=0;j<w;j++)
            printf("%f\t",m[i][j]);
        printf("\n");}
    printf("\n");
}

void switch_row(double **m,int w,int h,int r1,int r2)
{
    int i,j;
    double rt[w];
    for(j=0;j<w;j++)
    {
        rt[j]=m[r1][j];
        m[r1][j]=m[r2][j];
        m[r2][j]=rt[j];
    }
}

void solve(double **m,int n)
{
    int i,j,col;
    double p,c;

    for(col=0;col<4;col++)
    {
        p=m[col][col];
        if(p==0)
        {
            switch_row(m,n+1,n,col,col+1);
            col--;
            continue;
        }
        for(j=0;j<n+1;j++)m[col][j]/=p;

        for(i=0;i<n;i++)
        {
            c=m[i][col];
            if(i==col)continue;
            else for(j=0;j<n+1;j++) m[i][j]=m[i][j]-m[col][j]*c;
        }
    }
}

main()
{
    FILE *r,*w;
    float x[37],y[37],m,b,l;
    int i,j,ii,jj,a;
    double det,c1,c2,c3,c4;

    double **mat;
    mat=malloc((4)*sizeof(double*));
    for(i=0;i<4;i++)mat[i]=malloc((5)*sizeof(double));

    r=fopen("lum_angle_original.dat","r");
    w=fopen("lum_angle_cubic_interpolation.dat","w");

    for(i=0;i<37;i++)
    {
        fscanf(r,"%f",&y[i]);
        fscanf(r,"%f",&x[i]);
    }

    for(a=i=0;i<37;i++)
    {
        if(i==0||i==35||i==36)
        {
            c2=(y[i+1]-y[i])/(x[i+1]-x[i]);
            c1=y[i]-c2*x[i];
            c3=0;
            c4=0;
            goto A;
        }

        else
        {
            for(ii=0;ii<4;ii++)
                {
                    for(jj=0;jj<4;jj++)
                    {
                        mat[ii][jj]=pow(x[i-1+ii],jj);
                    }
                    mat[ii][4]=y[i-1+ii];
                }
        }

        solve(mat,4);

        c1=mat[0][4];
        c2=mat[1][4];
        c3=mat[2][4];
        c4=mat[3][4];

        A:
        for(a=i*5;a<(i+1)*5;a++)
        {
            l=c1+c2*a+c3*a*a+c4*a*a*a;
            fprintf(w,"%f\t%d\n",l,a);
            if (a==180)break;
        }
    }

    fclose(r);
    fclose(w);
}

 
[$[Get Code]]
Edit - History - Print - Recent Changes - Search
Page last modified on September 04, 2012, at 11:44 pm UTC