piperone.py

gistfile1.py

# Cosa dovrebbe fare questo script.
#
# Questo script dovrebbe prendere in input qualcosa di simile 
# al file che contiene tutti i datapoint, con i riferimenti 
# al numero della foto, e per ognuno di questi dovrebbe prendere
# l'appropriato tsv, quindi calcolare il raggio della circonferenza
# che ci interessa, e sputare il tutto fuori su stdout, in modo 
# che gnuplot possa plottare.

import scipy as sp
from scipy.optimize import leastsq
from scipy import sqrt

datapoints = sp.loadtxt('./phdata.tsv', skiprows=1, delimiter="\t")

for line in datapoints:
    circodati=sp.loadtxt('./digitized/'+str(int(line[0]))+'.tsv',\
        skiprows=1, delimiter="\t")

    data=circodati.transpose()
    guess=[data[0].mean(), data[1].mean()]
    
    circ = lambda p, x, y:  sqrt( (p[0]-x)**2 + (p[1]-y)**2 )
    resid = lambda p, x, y: circ(p, x, y) - circ(p, x, y).mean() 
        
    re, s = leastsq(resid, guess, args=(data[0], data[1]))
    #print(re, s)
    
    if (s>4):
        print('# questo cerchio non esiste' + str(line[0]) )
    else: 
        radius = (circodati-(re[0], re[1])).transpose()
        radius = sqrt(radius[0]**2 + radius[1]**2).mean()
        sep="\t"
        print(str(int(line[0]))+sep+str(line[1])+sep+str(line[2])+sep+\
            str(line[3])+sep+str(radius))