[2] | 1 | #!/bin/env python |
---|
| 2 | # coding: utf-8 |
---|
| 3 | |
---|
| 4 | # Calcul du dépointage à partir d'images du collimateur |
---|
| 5 | # Frederic.Meynadier@aerov.jussieu.fr Janvier 2008 |
---|
| 6 | |
---|
| 7 | |
---|
| 8 | import params |
---|
| 9 | |
---|
| 10 | import pyfits |
---|
| 11 | from optparse import OptionParser |
---|
| 12 | import sys |
---|
| 13 | |
---|
| 14 | from numpy import * |
---|
| 15 | from math import * |
---|
| 16 | |
---|
[5] | 17 | #from fitgaussian import * |
---|
[2] | 18 | |
---|
| 19 | VERSION='0.0' |
---|
| 20 | |
---|
| 21 | label = sys.argv[1] |
---|
| 22 | |
---|
| 23 | (ref,img,liste_points,sh_guess) = params.pset(label) |
---|
| 24 | |
---|
| 25 | |
---|
| 26 | sh_guess = array(sh_guess) |
---|
| 27 | |
---|
| 28 | size=128 |
---|
| 29 | subsize=size/2 |
---|
| 30 | |
---|
| 31 | for point in liste_points : |
---|
| 32 | limites = [size,2048-size] |
---|
| 33 | target = array(point) + sh_guess |
---|
| 34 | print target |
---|
| 35 | if (target[0] < limites[0] or target[0] > limites[1]\ |
---|
| 36 | or target[1] < limites[0] or target[1] > limites[1]): |
---|
| 37 | print "Point ",point," hors cadre dans l'image destination" |
---|
| 38 | sys.exit(1) |
---|
| 39 | |
---|
| 40 | cor_list = [] |
---|
| 41 | for coords in liste_points : |
---|
| 42 | |
---|
| 43 | tgt_ref=array(coords) |
---|
| 44 | tgt_img = tgt_ref + sh_guess |
---|
| 45 | |
---|
| 46 | |
---|
| 47 | |
---|
| 48 | |
---|
| 49 | ## parser = OptionParser() |
---|
| 50 | ## parser.add_option("-i","--input",dest="InputFilename", |
---|
| 51 | ## help="Input filename", metavar="RAWIMAGE") |
---|
| 52 | ## parser.add_option("-o","--output",dest="OutputFilename", |
---|
| 53 | ## help="Output Filename", metavar="FITSIMAGE") |
---|
| 54 | |
---|
| 55 | #(options, args) = parser.parse_args () |
---|
| 56 | |
---|
| 57 | |
---|
| 58 | # Ouverture de l'image de référence |
---|
| 59 | refhdulist = pyfits.open(ref) |
---|
| 60 | refdata = array(refhdulist[0].data[tgt_ref[1]-size/2:tgt_ref[1]+size/2,\ |
---|
| 61 | tgt_ref[0]-size/2:tgt_ref[0]+size/2]) |
---|
| 62 | refhdulist[0].data = refdata |
---|
| 63 | refhdulist.writeto('ref'+str(len(cor_list))+'.fits',clobber=True) |
---|
| 64 | |
---|
| 65 | |
---|
| 66 | # Ouverture de l'image à mesurer |
---|
| 67 | |
---|
| 68 | imghdulist = pyfits.open(img) |
---|
| 69 | imgdata = imghdulist[0].data[tgt_img[1]-subsize/2:tgt_img[1]+subsize/2,\ |
---|
| 70 | tgt_img[0]-subsize/2:tgt_img[0]+subsize/2] |
---|
| 71 | |
---|
| 72 | imghdulist[0].data = imgdata |
---|
| 73 | imghdulist.writeto('img'+str(len(cor_list))+'.fits',clobber=True) |
---|
| 74 | |
---|
| 75 | # Etablissement de la carte de corrélation |
---|
| 76 | |
---|
| 77 | corsize = size-subsize+1 |
---|
| 78 | |
---|
| 79 | cordata = zeros([corsize,corsize]) |
---|
| 80 | |
---|
| 81 | motif = sqrt(sum(sum(imgdata*imgdata))) |
---|
| 82 | |
---|
| 83 | for i in range(corsize): |
---|
| 84 | for j in range(corsize): |
---|
| 85 | refpart = refdata[i:i+subsize,j:j+subsize] |
---|
| 86 | cordata[i,j] = \ |
---|
| 87 | sum(sum( imgdata * refpart))\ |
---|
| 88 | / (motif * sqrt(sum(sum(refpart*refpart)))) |
---|
| 89 | |
---|
| 90 | |
---|
| 91 | cor_list.append(cordata) |
---|
| 92 | |
---|
| 93 | refhdulist[0].data = cordata |
---|
| 94 | prihdr = refhdulist[0].header |
---|
| 95 | |
---|
| 96 | prihdr.update('CRPIX1',(corsize-1)/2,'crpix1') |
---|
| 97 | prihdr.update('CRPIX2',(corsize-1)/2,'crpix2') |
---|
| 98 | prihdr.update('CRVAL1',sh_guess[0],'crval1') |
---|
| 99 | prihdr.update('CRVAL2',sh_guess[1],'crval2') |
---|
| 100 | prihdr.update('CDELT1',-1.0,'crdelt1') |
---|
| 101 | prihdr.update('CDELT2',-1.0,'crdelt2') |
---|
| 102 | prihdr.update('CROTA1',0.0,'crota1') |
---|
| 103 | prihdr.update('CROTA2',0.0,'crota2') |
---|
| 104 | prihdr.update('CTYPE1',' ','ctype1') |
---|
| 105 | prihdr.update('CTYPE2',' ','ctype2') |
---|
| 106 | |
---|
| 107 | |
---|
| 108 | refhdulist.writeto('cor'+str(len(cor_list)-1)+'.fits',clobber=True) |
---|
| 109 | |
---|
| 110 | |
---|
| 111 | |
---|
| 112 | cor_prod = cor_list[0] |
---|
| 113 | |
---|
| 114 | for i in range(1,len(cor_list)): |
---|
| 115 | cor_prod *= cor_list[i] |
---|
| 116 | |
---|
| 117 | |
---|
| 118 | refhdulist[0].data = cor_prod |
---|
| 119 | prihdr = refhdulist[0].header |
---|
| 120 | |
---|
| 121 | crpix1 = (corsize-1)/2 |
---|
| 122 | crpix2 = (corsize-1)/2 |
---|
| 123 | crval1 = sh_guess[0]+1 |
---|
| 124 | crval2 = sh_guess[1]+1 |
---|
| 125 | cdelt1 = -1.0 |
---|
| 126 | cdelt2 = -1.0 |
---|
| 127 | |
---|
| 128 | prihdr.update('CRPIX1',crpix1,'crpix1') |
---|
| 129 | prihdr.update('CRPIX2',crpix2,'crpix2') |
---|
| 130 | prihdr.update('CRVAL1',crval1,'crval1') |
---|
| 131 | prihdr.update('CRVAL2',crval2,'crval2') |
---|
| 132 | prihdr.update('CDELT1',cdelt1,'crdelt1') |
---|
| 133 | prihdr.update('CDELT2',cdelt2,'crdelt2') |
---|
| 134 | prihdr.update('CROTA1',0.0,'crota1') |
---|
| 135 | prihdr.update('CROTA2',0.0,'crota2') |
---|
| 136 | prihdr.update('CTYPE1',' ','ctype1') |
---|
| 137 | prihdr.update('CTYPE2',' ','ctype2') |
---|
| 138 | |
---|
| 139 | |
---|
| 140 | refhdulist.writeto('cor_prod.fits',clobber=True) |
---|
| 141 | |
---|
| 142 | refhdulist.close() |
---|
| 143 | imghdulist.close() |
---|
| 144 | |
---|
| 145 | # On cherche le maximum du tableau -> passage en 1D |
---|
| 146 | lin_max = argmax(cor_prod.ravel()) |
---|
| 147 | x_max = lin_max % corsize |
---|
| 148 | y_max = floor(lin_max / corsize) |
---|
| 149 | #print "Pixel max :",(int(x_max),int(y_max)) |
---|
| 150 | #print "Décalage (entier) correspondant : (%d,%d)" %\ |
---|
| 151 | # (int((x_max-crpix1)*cdelt1+crval1-1), |
---|
| 152 | # int((y_max-crpix2)*cdelt2+crval2-1)) |
---|
| 153 | |
---|
| 154 | # On cherche le centroïde (bloc 2*(npfit)+1**2) |
---|
| 155 | npfit = 1 |
---|
| 156 | |
---|
| 157 | |
---|
| 158 | bloc = cor_prod[y_max-npfit:y_max+npfit+1, |
---|
| 159 | x_max-npfit:x_max+npfit+1] |
---|
| 160 | bloc = bloc - min(bloc.ravel()) |
---|
| 161 | |
---|
| 162 | og = array([0.0,0.0]) |
---|
| 163 | |
---|
| 164 | for i in range(2*npfit+1): |
---|
| 165 | for j in range(2*npfit+1): |
---|
| 166 | vec = array([j+1,i+1]) |
---|
| 167 | og += vec*bloc[i,j] |
---|
| 168 | |
---|
| 169 | og /= sum(sum(bloc)) |
---|
| 170 | og -=1 |
---|
| 171 | |
---|
| 172 | |
---|
| 173 | print "ref :",ref |
---|
| 174 | print "img :",img |
---|
| 175 | |
---|
| 176 | #print "Centroïde en : %f %f" % (og[0]+x_max,og[1]+y_max) |
---|
| 177 | print "Décalage centroïde : %f %f" % (((og[0]+x_max)-crpix1)*cdelt1+crval1, |
---|
| 178 | ((og[1]+y_max)-crpix2)*cdelt2+crval2) |
---|
| 179 | |
---|
| 180 | |
---|
[5] | 181 | #params = fitgaussian(bloc) |
---|
| 182 | #fit = gaussian(*params) |
---|
| 183 | #(height, x, y, width_x, width_y) = params |
---|
| 184 | #print "Décalage gaussienne : %f %f" % (((x+x_max)-crpix1)*cdelt1+crval1, |
---|
| 185 | # ((y+y_max)-crpix2)*cdelt2+crval2) |
---|
[2] | 186 | |
---|