source: mire/process.py @ 16

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

mise au point pour fonctionnement batch

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