source: mire/process.py @ 14

Last change on this file since 14 was 14, checked in by meynadie, 12 years ago

Ajouts logiciels mire

  • Property svn:executable set to *
File size: 6.9 KB
Line 
1#!/bin/env python
2# coding: utf-8
3# Un programme pour exploiter les données des mires
4# -- Frédéric Meynadier
5
6import string
7from numpy import *
8
9from functions import *
10
11from optparse import OptionParser
12parser = OptionParser()
13parser.add_option("-c","--catalog",dest="catalogfile",
14                  help="Name of the catalog file")
15parser.add_option("-o","--outputfits",dest="outputfits",
16                  default="output.fits",
17                  help="Name of the output fits file")
18
19
20(options, args) = parser.parse_args ()
21
22# Donner les coordonnées approximatives des 4 repÚres
23
24rep = zeros((4,2))
25i = 0
26for line in open('reperes.txt'):
27    if ((line[0] != "#") and (line != '\n')):
28        spl = string.split(line)
29        rep[i][0] = float(spl[0])
30        rep[i][1] = float(spl[1])
31        i+=1
32
33# Lecture du fichier catalogue
34fields = ["FLUX_MAX",
35          "X_IMAGE",
36          "Y_IMAGE",
37          "THETA_IMAGE",
38          "ELONGATION",
39          "ELLIPTICITY",
40          "FWHM_IMAGE",
41          "ERRX2_IMAGE",
42          "ERRY2_IMAGE"]
43
44catalog = []
45for line in open(options.catalogfile):
46    if ((line[0] != "#") and (line != '\n')):
47        tmp = {}
48        spl = string.split(line)
49        for i in range(len(fields)):
50            tmp[fields[i]] = float(spl[i])
51        catalog.append(tmp)
52
53
54# Recherche de la position des repÚres selon catalogue - attention aux doublons
55real_rep = zeros((4,2))
56for i in range(4):
57    detections = 0
58    for source in catalog:
59        pos_source = array([source['X_IMAGE'],source['Y_IMAGE']])
60        if ((modulus(rep[i,:] - pos_source)) < 5):
61            detections +=1
62            real_rep[i,:] = pos_source
63    if (detections !=1):
64        print "Attention ! ambiguïté pour repÚre ",i,\
65              ", nb de détections :",detections
66       
67
68# Détermination du centre : barycentre des 4 repÚres
69center = array([sum(real_rep[:,0])/4,sum(real_rep[:,1])/4])
70
71# Vecteurs unitaires de référence
72vec_i = ((real_rep[0,:] - real_rep[1,:])/20 \
73         + (real_rep[3,:] - real_rep[2,:])/20) / 2 
74vec_j = ((real_rep[1,:] - real_rep[2,:])/20 \
75         + (real_rep[0,:] - real_rep[3,:])/20) / 2 
76
77# Détection des noeuds de grille les plus proches de chaque source
78vec_i_2 = dot(vec_i,vec_i)
79vec_j_2 = dot(vec_j,vec_j)
80for source in catalog:
81    pos_source = array([source['X_IMAGE'],source['Y_IMAGE']])-center
82    # dans le référentiel grille (en projetant sur les vecteurs unitaires) :
83    pos_source_grid = array([dot(pos_source,vec_i)/vec_i_2,
84                             dot(pos_source,vec_j)/vec_j_2])
85    # on arrondit au noeud le plus proche
86    knot = pos_source_grid.round()
87    # et on calcule la distance entre les deux
88    r = modulus(pos_source_grid - knot)
89    # enregistrement des résultats
90    source['I_GRID'] = knot[0]
91    source['J_GRID'] = knot[1]
92    source['KNOT_DIST'] = r
93    #if (modulus(knot) <1):
94    #    print source['I_GRID'], source['J_GRID'], source['KNOT_DIST'], pos_source_grid, pos_source,source['X_IMAGE'], source['Y_IMAGE']
95
96# Chasse aux doublons
97# tableau 81*81 contenant :
98# champ 0 : nb de sources rattachée au noeud
99# champ 1 : distance actuelle entre le noeud et la source
100# champ 2 : index courant de la source dans catalog
101# champ 3 : distance moyenne aux voisins (calculé plus tard)
102lookup_sources = zeros((81,81,4))
103for si in range(len(catalog)):
104    source = catalog[si]
105    i = source['I_GRID']+40
106    j = source['J_GRID']+40
107    if ((lookup_sources[i,j,0] < 1)
108        or (lookup_sources[i,j,1] > source['KNOT_DIST'])):
109        lookup_sources[i,j,2] = si
110        lookup_sources[i,j,1] = source['KNOT_DIST']
111    lookup_sources[i,j,0] += 1
112
113# On élimine le noeud (1,1), trop prÚs du "P"
114lookup_sources[41,41,0] = 0
115
116# Calcul de la distance moyenne entre noeuds voisins
117
118for i in range(1,80):
119    for j in range(1,80):
120        if (lookup_sources[i,j,0] > 0):
121            knot_idx = int(lookup_sources[i,j,2])
122            catalog[knot_idx]['MEAN_DIST_TO_NGBRS'] = 0
123            nb_voisin = 0
124            distance = 0
125            pos_source = array([catalog[knot_idx]['X_IMAGE'],
126                                catalog[knot_idx]['Y_IMAGE']])
127            for shifts in [[0,1],[1,0],[-1,0],[0,-1]]:
128                i_voisin = i + shifts[0]
129                j_voisin = j + shifts[1]
130                if (lookup_sources[i_voisin,j_voisin,0] > 0):
131                    neigh_idx = int(lookup_sources[i_voisin,j_voisin,2])
132                    nb_voisin +=1
133                    pos_voisin = array([catalog[neigh_idx]['X_IMAGE'],
134                                        catalog[neigh_idx]['Y_IMAGE']])
135                    distance += modulus(pos_source-pos_voisin)
136            if (nb_voisin !=0):
137                catalog[knot_idx]['MEAN_DIST_TO_NGBRS'] = distance / nb_voisin
138                lookup_sources[i,j,3] = distance / nb_voisin
139       
140cleaned_catalog = []
141for i in range(81):
142    for j in range(81):
143        if (lookup_sources[i,j,0] > 0):
144            cleaned_catalog.append(catalog[int(lookup_sources[i,j,2])])
145
146# Sortie du catalogue nettoyé des doublons
147
148of = open('clean.cat','w')
149
150fields = ["FLUX_MAX",
151          "X_IMAGE",
152          "Y_IMAGE",
153          "THETA_IMAGE",
154          "ELONGATION",
155          "ELLIPTICITY",
156          "FWHM_IMAGE",
157          "ERRX2_IMAGE",
158          "ERRY2_IMAGE",
159          "I_GRID",
160          "J_GRID",
161          "KNOT_DIST",
162          "MEAN_DIST_TO_NGBRS"]
163
164fmts = ["%10.3f ",
165        "%8.3f ",
166        "%8.3f ",
167        "%8.3f ",
168        "%8.3f ",
169        "%8.3f ",
170        "%8.3f ",
171        "%8.5f ",
172        "%8.5f ",
173        "%5d ",
174        "%5d ",
175        "%5.3f ",
176        "%8.3f"]
177
178for source in cleaned_catalog:
179    outs = ""
180    for i in range(len(fields)):
181        outs += fmts[i] % source[fields[i]]
182    outs += "\n"
183    of.write(outs)
184
185of.close()
186
187# Génération d'un fichier fits contenant les résultats
188
189
190import pyfits
191
192#a = zeros((2048,2048))
193#
194# patches width in pixel:
195## width = int(modulus(vec_i))
196
197## for i in range(len(cleaned_catalog)):
198##     print i
199##     source = cleaned_catalog[i]
200##     pos_source = array([int(source['X_IMAGE']),int(source['Y_IMAGE'])])
201##     pos_mean_dist = source['MEAN_DIST_TO_NGBRS']
202##     for rel_x in range(-width,width):
203##         for rel_y in range(-width,width):
204##             rel_vec = array([rel_x,rel_y])
205##             # Forme du patch : carré
206##             if ( (dot(rel_vec,vec_i) < vec_i_2/2) and
207##                  (dot(rel_vec,vec_j) < vec_j_2/2)):
208##                 (x,y) = pos_source+rel_vec
209##                 if ( x>0 and x<2048 and y>0 and y<2048):
210##                     a[x,y] = pos_mean_dist
211
212
213
214a = lookup_sources[:,:,3]
215           
216hdu = pyfits.PrimaryHDU(a)
217hdulist = pyfits.HDUList([hdu])
218
219hdu.writeto(options.outputfits,clobber=True)
220
221# Facteur d'échelle moyen
222
223tmp = []
224for i in range(81):
225    for j in range(81):
226        if (lookup_sources[i,j,0] > 0):
227            tmp.append(lookup_sources[i,j,3])
228
229print "Moyenne des mesures : ",sum(tmp)/len(tmp)
230
231# Module des vecteurs unitaires
232
233print "Module vec_i :",modulus(vec_i)
234print "Module vec_j :",modulus(vec_j)
Note: See TracBrowser for help on using the repository browser.