source: cherche_centre/04_cherche_centre.py @ 34

#!/bin/env python
# encoding: utf-8

import pyfits
from optparse import OptionParser

import numpy
from scipy import *

parser = OptionParser()

parser.add_option("-t","--threshold",dest="threshold", default=0.5,
                  help="threshold for edge detection (relative to max, default 0.5)")

parser.add_option("-w","--width",dest="hw_fen_int", default=5,
                  help="window for measuring central intensity, size will be 2*hw_fen_int+1")

(options, args) = parser.parse_args ()

if (len(args) != 1):
    print "Usage : cherche_centre.py [-t threshold] [-w width] image.fits"
    exit(1)


# Ouverture de l'image fits en argument
hdulist = pyfits.open(args[0])
scidata = hdulist[0].data


# Estimation fond (4 coins)
tf = 256 # taille fond, en pix
fond_bas = .5 * (mean(scidata[:tf,:tf]) + mean(scidata[:tf,2048-tf:]))
fond_haut= .5 * (mean(scidata[2048-tf:,:tf]) + mean(scidata[2048-tf:,2048-tf:]))


if (abs(fond_haut - fond_bas) > .5 * (fond_haut + fond_bas) ):
    print "Ereur : offsets déséquilibrés. Corriger le fond."
    exit(3)

def filtre_median(profil,hw):
    tmp = numpy.zeros(2048)
    for i in range(hw,2048-hw-1):
        tmp[i] = median(profil[i-hw:i+hw+1])
    return tmp

# Moyenne suivant lignes et colonnes
moyenne_h = filtre_median(sum(scidata, axis=0) / 2048.,5) 
moyenne_v = filtre_median(sum(scidata, axis=1) / 2048.,5)



of = open("profils"+args[0][:-5]+".csv",'w')
for i in range(2048):
    of.write("%4i %f %f\n" % (i, moyenne_h[i], moyenne_v[i]))
of.close()

def get_rayon_centre(profil):
    # Détermination du seuil, ici 10%
    seuil = max(profil) * float(options.threshold)
    # Méthode bête
    bord1 = 0
    for i in range(2047):
        if ((profil[i] < seuil) and (profil[i+1] >= seuil)):
            bord1 = i

    bord2 = 2048
    for i in range(2047,1,-1):
        if((profil[i] < seuil) and (profil[i-1] >= seuil)):
            bord2 = i

    flag=""
    if ((bord1 == 0) or (bord2 == 2048)) :
        flag = "*"
        print "Attention ! un seul bord détecté !"
        rayon = 0
        centre = 0
    else:
        rayon = (bord2 - bord1) / 2.
        centre = (bord1 + bord2) / 2.
    return (rayon,centre,flag)


(rayon_h,centre_h,flag) = get_rayon_centre(moyenne_h)
(rayon_v,centre_v,flag) = get_rayon_centre(moyenne_v)


hw_fen_int = int(options.hw_fen_int)

int_centre = sum(sum(scidata[int(centre_v)-hw_fen_int:int(centre_v)+hw_fen_int+1,
                             int(centre_h)-hw_fen_int:int(centre_h)+hw_fen_int+1]))/((2*hw_fen_int + 1)**2)
                             


print "%13s %4i %6.1f %4i %6.1f %6.0f %1s" % (args[0],
                                              rayon_h,
                                              centre_h+1, # fits commence à 1
                                              rayon_v,
                                              centre_v+1, # fits commence à 1
                                              int_centre,
                                              flag)