In [ ]:
import proper
import numpy as np
import scipy as sc
import matplotlib.pyplot as plt
%matplotlib inline

import cv2 as cv
In [ ]:
from skimage import data, color
from skimage.transform import hough_circle
from skimage.feature import peak_local_max, canny
from skimage.draw import circle_perimeter, circle
from skimage.util import img_as_ubyte
from scipy import ndimage

from skimage.transform import resize, downscale_local_mean

# old algorithm levels
#threshold1 = 0.16
#threshold2 = 0.205

#threshold1 = 0.18
#threshold2 = 0.24

#threshold1 = 0.03
#threshold2 = 0.0454

#threshold1 = 0.11
#threshold2 = 0.118

#threshold1 = 0.3
#threshold2 = 0.3

# ok for unoccluded
threshold1 = 0.18
threshold2 = 0.27

# ok for occluded
#threshold1 = 0.27
#threshold2 = 0.38

import random as rng
rng.seed(12345)

def thresh_callback(src_gray, val):
    threshold = val
    _, contours, _ = cv.findContours(src_gray, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE)
    contours_poly = [None]*len(contours)
    centers = [None]*len(contours)
    radius = [None]*len(contours)
    for i, c in enumerate(contours):
        contours_poly[i] = cv.approxPolyDP(c, 3, True)
        centers[i], radius[i] = cv.minEnclosingCircle(contours_poly[i])
    drawing = np.zeros((src_gray.shape[0], src_gray.shape[1], 3), dtype=np.uint8)
    for i in range(len(contours)):
        colore = (rng.randint(0,256), rng.randint(0,256), rng.randint(0,256))
        cv.drawContours(drawing, contours_poly, i, colore)
        cv.circle(drawing, (int(centers[i][0]), int(centers[i][1])), int(radius[i]), colore, 2)
    print("Max r", np.max(radius))
    fig, ax = plt.subplots(1,1, figsize=(15,10))
    ax.imshow(drawing + color.gray2rgb(src_gray) )
    plt.show()
    return np.max(radius)
    
def plateauLevel(input_image):
    imm = np.copy(input_image)
    imax = (np.max(input_image))
    imm *= 1.0/imax
    imm_cap = np.where(imm < threshold1, 0.0, 1.0)
    count = np.sum(imm_cap)
    guessR = np.sqrt(float(count)/np.pi)
    guessCenter = ndimage.measurements.center_of_mass(imm_cap)
    plateauRadius = 0.8 * guessR # 0.8 * guessR
    print(guessR)
    print(guessCenter)    
    # fig, ax = plt.subplots(1,1, figsize=(15,10))
    # ax.imshow(input_image)
    # plt.show()
    immInnerDisc = input_image*0.0
    rr, cc = circle(guessCenter[0], guessCenter[1], plateauRadius)
    immInnerDisc[rr, cc] = 1.0
    plateauImage = np.where(immInnerDisc >=0.97, input_image, 0.0)
    plateauImageMask = np.where(immInnerDisc >=0.97, 1.0, 0.0)
    #fig, ax = plt.subplots(1,1, figsize=(15,10))
    #ax.imshow(plateauImage)
    #plt.show()
    plateauLevelOut = np.sum(plateauImage)/np.sum(plateauImageMask)
    print('plateauLevel:', plateauLevelOut/imax)
    return plateauLevelOut

def computeDiameter(input_image):    
    imm = np.copy(input_image)    
    threshold = plateauLevel(imm) 
    imm *= 1.0/threshold
    imm_cap = np.where(imm < threshold2, 0.0, 1.0)
    count = np.sum(imm_cap)
    guess = 2*int(np.sqrt(float(count)/np.pi))
    # print(guess)
    #plt.imshow(imm_cap, cmap='hot')
    soglie = [0.5] 
    cx = [0]
    cy = [0]
    ii = 0
    image = (input_image/np.max(input_image))
#    imm_cap = resize(imm_cap, (4096, 4096), order=1)
#    imm_cap = resize(imm_cap, (2*2472, 2*3296), order=1)
    imm_cap = resize(imm_cap, (2*imm_cap.shape[0], 2*imm_cap.shape[1]) , order=1)

    sx = ndimage.sobel(imm_cap, axis=0, mode='constant',cval=0)
    sy = ndimage.sobel(imm_cap, axis=1, mode='constant',cval=0)
    sob = np.hypot(sx, sy)
    sob = sob/np.max(sob)
    idx = sob[:,:] > soglie[ii]
    sob *= 0.0; 
    sob[idx] = 1
    edges = img_as_ubyte(sob)
    image_C = color.gray2rgb(edges)
    
#    thresh = 50 # 25 # initial threshold
#    rr_alt = thresh_callback(edges, thresh)
#    return rr_alt # radius-0.5
    
    hough_radii = np.arange(guess-1, guess+4, 1)
    hough_res = hough_circle(edges, hough_radii)
    centers = []
    accums = []
    radii = []
    image_C = downscale_local_mean(edges, (4,4))
    image_C = color.gray2rgb(image_C)

    for radius, h in zip(hough_radii, hough_res):
        num_peaks = 1
        peaks = peak_local_max(h, num_peaks=num_peaks)
        centers.extend(peaks)
        accums.extend(h[peaks[:, 0], peaks[:, 1]])
        radii.extend([radius] * num_peaks)
        # print(radii)
    for idx in np.argsort(accums)[::-1][:1]:
        center_x, center_y = centers[idx]
        center_x = center_x//4
        center_y = center_y//4
        radius = radii[idx]
        ccx, ccy = circle_perimeter(center_y, center_x, radii[idx]//5)
        image_C[ccy, ccx] = (255, 0, 0)
    center_x, center_y = centers[0]
    center_x = center_x//4
    center_y = center_y//4
    #radius = radii[0]
    print('radius', radius)
    sss = 900//4
    fig, ax = plt.subplots(1,1, figsize=(15,10))
#    ax.imshow(image_C[center_x-sss:center_x+sss, center_y-sss:center_y+sss])
    ax.imshow(image_C) # [center_x-sss:center_x+sss, center_y-sss:center_y+sss])
    plt.show()
    cx[ii] = center_x
    cy[ii] = center_y
    
    # print('centri', float(center_x)/2.0+0.5, float(center_y)/2.0+0.5)
    return radius-0.5
In [ ]:
from scipy import stats

def fitAndShow(xx, dd, method=1):
    if method==1:
        a,b = np.polyfit(dd, xx, 1)
    else:
        a, b, r_value, p_value, std_err = stats.linregress(dd, xx)
        # print('p_value**2:', p_value**2, 'std_err:', std_err)
    print('a:', format(a, '.4f'))
    # print('1/a:', format(1/a, '.4f'))
    print('b:', format(b, '.4f'))
    
    N = float(dd.shape[0])
    deltaPrime = N * np.sum(np.multiply(dd,dd)) - np.sum(dd)**2
    err = xx - a*dd - b
    sigmaSquare = (1.0/(N-2)) * np.sum(np.multiply(xx,xx))
    
    aSigma = (sigmaSquare/deltaPrime) * np.sum(np.multiply(dd,dd))
    print( 'aSigma:', format(np.sqrt(aSigma), '.4f') )
    print( '(aSigma/a)*100:', format(100*np.sqrt(aSigma)/a, '.3f') )
    
    
    plt.figure(figsize=(15,10))
    plt.plot(dd, xx, 'yo', dd, a*dd+b, '--k')
    plt.show()
In [4]:
def simple_system(wavelength, gridsize, br, offset, aperture_diameter, f_number):
    focal_ratio = f_number
    diam = aperture_diameter
    focal_length = diam * focal_ratio
    beam_ratio = br
    wfo = proper.prop_begin(diam, wavelength, gridsize, beam_ratio)
    proper.prop_circular_aperture(wfo, diam/2.0)
    proper.prop_define_entrance(wfo)
    proper.prop_lens(wfo, focal_length)
    proper.prop_propagate(wfo, focal_length + offset, TO_PLANE=False) 
    gridsize = proper.prop_get_gridsize( wfo )
    print( proper.prop_get_fratio( wfo ))
    (wfo, sampling) = proper.prop_end(wfo)  
    return (wfo, sampling, gridsize)
In [5]:
def occluded_system(wavelength, gridsize, br, offset, aperture_diameter, f_number):
    focal_ratio = f_number
    diam = aperture_diameter
    focal_length = diam * focal_ratio
    beam_ratio = br
    wfo = proper.prop_begin(diam, wavelength, gridsize, beam_ratio)
    proper.prop_circular_aperture(wfo, diam/2)                  # aperture (primary mirror)

    # Normalizzato a 1 il diametro esterno della pupilla si ha:
    # spessore dello spider: 9.4e-3
    # diametro oscuramento centrale: 0.094
    proper.prop_circular_obscuration(wfo, 0.094 * diam/2)       # secondary mirror obscuration                                                                       
                                                                # diam/2== fully obscured
    # secondary vane (vertical)
    proper.prop_rectangular_obscuration(wfo, 9.4e-3 * diam, diam, 0,  0, NORM=False, ROTATION=45.0)
    # secondary vane (horizontal)
    proper.prop_rectangular_obscuration(wfo, diam, 9.4e-3 * diam, 0,  0, NORM=False, ROTATION=45.0)

    proper.prop_define_entrance(wfo)
        
    proper.prop_lens(wfo, focal_length)
    proper.prop_propagate(wfo, focal_length + offset, TO_PLANE=False) 
    gridsize = proper.prop_get_gridsize( wfo )
    print( proper.prop_get_fratio( wfo ))
    (wfo, sampling) = proper.prop_end(wfo)  

    return (wfo, sampling, gridsize)
In [6]:
def funcFilter(x, y):
    return 0.997916 - 0.0491551*y - 0.414188*(y**2) + 0.0322198*x + 0.000797465*x*y - 0.0105691*x*(y**2)- 0.411941*(x**2) + 0.0124160*(x**2)*y + 0.0507806*(x**2)*(y**2)


def apoFilter():
    xaxis = np.linspace(-2, 2, 2048)
    yaxis = np.linspace(-2, 2, 2048)
    x, y = np.meshgrid(xaxis, yaxis)
    return funcFilter(x, y) 
    
    
def gaussian_occluded_system(wavelength, gridsize, br, offset, aperture_diameter, f_number):
    focal_ratio = f_number
    diam = aperture_diameter
    focal_length = diam * focal_ratio
    beam_ratio = br
    wfo = proper.prop_begin(diam, wavelength, gridsize, beam_ratio)
    proper.prop_circular_aperture(wfo, diam/2)                  # aperture (primary mirror)

    occrad_m = aperture_diameter/2.0
    r = proper.prop_radius(wfo)
    print('occrad_m:', occrad_m)
    h = np.sqrt(-0.5 * occrad_m**2 / np.log(1 - np.sqrt(0.5)))        
    print('h:', h)
    gauss_spot = np.exp(-0.5 * (r/h)**2)         
    custom_spot = gauss_spot.copy()
    custom_spot *= 0.0    
    custom_spot = apoFilter()
    custom_spot = np.where(custom_spot < 0.0, 0.0, custom_spot)
    plt.imshow(gauss_spot, cmap='hot')
    plt.show()
    plt.imshow(custom_spot, cmap='hot')
    plt.colorbar()
    plt.show()
    proper.prop_multiply(wfo, custom_spot)             
    proper.prop_define_entrance(wfo)
    proper.prop_lens(wfo, focal_length)
    proper.prop_propagate(wfo, focal_length + offset, TO_PLANE=False) 
    gridsize = proper.prop_get_gridsize( wfo )
    print( proper.prop_get_fratio( wfo ))
    (wfo, sampling) = proper.prop_end(wfo)  

    return (wfo, sampling, gridsize)
In [7]:
def runmodel_off(model, Fnumber, system=simple_system):
    gridsize = model['grid_size'] 
    wavelengths = model['wavelengths'] 
    beam_ratio = model['beam_ratio'] 
    offsets = model['offset'] 
    aperture_diameter = model['aperture_diameter'] 
    r1 = [0]*len(offsets)
    s1 = [0]*len(offsets)
    gs1 = [0]*len(offsets)
    ii=0
    tot= len(offsets)*len(wavelengths)
    for i, offset in enumerate(offsets): 
        r1[i] = np.zeros([gridsize, gridsize])
        for j, w in enumerate(wavelengths): 
            r1ij, s1j, gs = system(w, gridsize, beam_ratio, offset, aperture_diameter, Fnumber)
            ii+=1
            print("Done:%f2" % (100* ii/tot))
            print(gs)
        gs1[i] =gs            
        s1[i] =s1j            
        r1[i] += r1ij
    return r1, s1, gs1
In [8]:
def runmodel_w(model, Fnumber):
    gridsize = model['grid_size'] 
    wavelengths = model['wavelengths'] 
    beam_ratio = model['beam_ratio'] 
    offsets = model['offset'] 
    aperture_diameter = model['aperture_diameter'] 
    r1 = [0]*len(wavelengths)
    s1 = [0]*len(wavelengths)
    gs1 = [0]*len(wavelengths)
    ii=0
    tot = len(offsets)*len(wavelengths)
    offset = offsets[0]
    for j, w in enumerate(wavelengths): 
        r1[j], s1[j], gs1[j] = simple_system(w, gridsize, beam_ratio, offset, aperture_diameter, Fnumber)
        print("Done:%f2" % (100* j/tot))
    return r1, s1, gs1
In [9]:
def showResult_offset(model, r1, s1, showImages=False):
    gridsize = r1[0].shape[0]
    gridsizeF = float(r1[0].shape[0])
    if showImages:
        fig, (ax1, ax2, ax3, ax4) = plt.subplots(4,1, figsize=(15,60))
    else:
        fig, (ax4) = plt.subplots(1,1, figsize=(15,10))
    if showImages:
        ax1.matshow(r1[0])
        ax2.matshow(r1[1])
        ax3.matshow(r1[2])
#        ax5.matshow(r1[2]-r1[3])
#        ax6.matshow(r1[1]-r1[4])
#        ax7.matshow(r1[0]-r1[5])
#        print("Max diff:", np.mean( np.abs(r1[0]-r1[5]) ) )
#        print("Max image:", np.mean( np.abs(r1[0]) ) )

    FN = [0] * len(r1)
    for i in range(len(r1)):
        section = (r1[i][:, gridsize//2])
        norm_section = (1.0/np.max(section))*section
        x_coords = np.array(np.linspace( 0, s1[i]*gridsizeF, gridsize)) - s1[i]*gridsizeF/2.0 
        section_cap = np.where(norm_section < 0.155, 0.0, 1.0)
#        ax4.plot(np.linspace( -s1[i]*gridsize//2, s1[i]*(gridsize - gridsize//2), gridsize), section)
        ax4.set_yticks(np.arange(0, 1., 0.05))
        plt.ylim((0., 1))
        ax4.plot(x_coords, norm_section)
#        ax4.plot(x_coords, section_cap)
        section_cap_list = section_cap.tolist()
        section_cap_diff = np.array([section_cap_list[i] - section_cap_list[i+1] for i in range(len(section_cap_list)-1)])
        # ax4.plot(x_coords[:-1], section_cap_diff)
        #rx = x_coords[np.where(norm_section>=0.999)]        
        rx = x_coords[np.where(section_cap_diff>=0.999)]
        FN[i] = (-model['offset'][i])/(2.0*rx)
    print("FN", FN)      
    #print("gs1", gs1)
    plt.grid()
    plt.show()
    
def showResult_w(model, r1, s1):
    gridsize = r1[0].shape[0]
    gridsizeF = float(r1[0].shape[0])
    fig, (ax4) = plt.subplots(1,1, figsize=(15,10))
    FN = [0] * len(r1)
    for i in range(len(r1)):
        section = (r1[i][:, gridsize//2])
        norm_section = (1.0/np.max(section))*section
        x_coords = np.array(np.linspace( 0, s1[i]*gridsizeF, gridsize)) - s1[i]*gridsizeF/2.0 
        section_cap = np.where(norm_section < 0.157, 0.0, 1.0)
        ax4.plot(x_coords[gridsize//5:2*gridsize//5], norm_section[gridsize//5:2*gridsize//5])
        section_cap_list = section_cap.tolist()
        section_cap_diff = np.array([section_cap_list[i] - section_cap_list[i+1] for i in range(len(section_cap_list)-1)])
        rx = x_coords[np.where(section_cap_diff>=0.999)]
    print(FN)      
    print(gs1)
    plt.show()
In [10]:
model1a = {}
model1a['grid_size'] = 4096
model1a['wavelengths'] = np.linspace(500*1e-9, 1600*1e-9, 1)
model1a['beam_ratio'] = 0.5
model1a['offset'] = np.linspace(-0.09, -0.08, 3) # -0.12
model1a['aperture_diameter'] = 0.011

#rrt1, r1 = runmodel_w(model1)
r1a, s1a, gs1a = runmodel_off(model1a, 13.5)
print(s1a)

Applying lens
Propagating
13.499991759075865
Done:33.3333332
4096
Applying lens
Propagating
13.499991759075865
Done:66.6666672
4096
Applying lens
Propagating
13.499991759075865
Done:100.0000002
4096
[3.255207041713938e-06, 3.07436202391171e-06, 2.893517006109483e-06]
In [11]:
showResult_offset(model1a, r1a, s1a)
dda = [0] * len(s1a) 
for i in range(len(s1a)):
    dda[i] = computeDiameter(r1a[i])
    print((-model1a['offset'][i])/(dda[i]*s1a[i]))

FN [array([13.46383876]), array([13.4769686]), array([13.46384045])]

1023.9089701547953
(2048.0, 2048.0)
plateauLevel: 0.7228932003572651

D:\Fabio\Anaconda\lib\site-packages\skimage\transform\_warps.py:105: UserWarning: The default mode, 'constant', will be changed to 'reflect' in skimage 0.15.
  warn("The default mode, 'constant', will be changed to 'reflect' in "
D:\Fabio\Anaconda\lib\site-packages\skimage\transform\_warps.py:110: UserWarning: Anti-aliasing will be enabled by default in skimage 0.15 to avoid aliasing artifacts when down-sampling images.
  warn("Anti-aliasing will be enabled by default in skimage 0.15 to "
D:\Fabio\Anaconda\lib\site-packages\skimage\util\dtype.py:141: UserWarning: Possible precision loss when converting from float64 to uint8
  .format(dtypeobj_in, dtypeobj_out))

radius 2046

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

13.51650499649438
1023.3461278450552
(2048.0, 2048.0)
plateauLevel: 0.6969864371628179
radius 2045

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

13.523116952033979
1022.7680371611195
(2048.0, 2048.0)
plateauLevel: 0.6752141540029214
radius 2045

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

13.52311785341205
In [14]:
from astropy.io import fits

fitsdata = fits.open('C:\\Users\\Fabio\\JupyterNotebooks\\f13-11mm-out-1mm.fits')
fitsdata.info()
images = fitsdata[0].data
images_d = images.copy()
plt.figure(figsize=(10,10))
cube = images_d[1:14].copy()
#imm[images_d[0]<300] = 0
plt.imshow(cube[12], cmap='hot')

Filename: C:\Users\Fabio\JupyterNotebooks\f13-11mm-out-1mm.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU      10   (2048, 2048, 20)   float32
Out[14]:
<matplotlib.image.AxesImage at 0x1ed906a04e0>
In [15]:
dda = [0] * 13
for i in range(13):
    dda[i] = computeDiameter(cube[i])
    
print(dda)

399.98614350211795
(1116.4192710198558, 1161.8020950220844)
plateauLevel: 0.584169152137475
radius 804

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

405.0154298740586
(1117.0330151608941, 1162.9240732799187)
plateauLevel: 0.5806486970842638
radius 813

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

408.70658243059745
(1117.2876947263112, 1165.362999380687)
plateauLevel: 0.5260888745739952

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

radius 824

413.78988735482994
(1117.6410849398599, 1166.7521741555279)
plateauLevel: 0.5310345167650405

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

radius 833

418.17677127275516
(1118.3586632834342, 1168.7262257543102)
plateauLevel: 0.5046310953835111

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

radius 843

423.0124041534112
(1118.8548994494402, 1170.656400814722)
plateauLevel: 0.4992180380701166
radius 854

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

428.51712387207044
(1119.8344181902332, 1172.0991486979117)
plateauLevel: 0.5150709069379469
radius 864

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

433.47534219766305
(1119.7740230523727, 1174.1577769571138)
plateauLevel: 0.5312349920168218
radius 873

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

439.65756145355425
(1120.610622033837, 1175.8857024763447)
plateauLevel: 0.5730552299743928

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

radius 884

444.6558471704088
(1121.3700446911546, 1177.2141858997475)
plateauLevel: 0.5733741671345949
radius 893

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

448.3387114301919
(1122.678470080002, 1178.6677097123602)
plateauLevel: 0.5191332973595639
radius 904

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

454.32852859313465
(1123.1319814332198, 1180.7125757552392)
plateauLevel: 0.5557062439246394
radius 914

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

460.0999780553461
(1123.8254883091497, 1182.6701947222014)
plateauLevel: 0.586806301396305
radius 924

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

[803.5, 812.5, 823.5, 832.5, 842.5, 853.5, 863.5, 872.5, 883.5, 892.5, 903.5, 913.5, 923.5]
In [16]:
from scipy import stats

def fitAndShow(xx, dd, method=1):
    if method==1:
        a,b = np.polyfit(dd, xx, 1)
    else:
        a, b, r_value, p_value, std_err = stats.linregress(dd, xx)
        # print('p_value**2:', p_value**2, 'std_err:', std_err)
    print('a:', format(a, '.4f'))
    # print('1/a:', format(1/a, '.4f'))
    print('b:', format(b, '.4f'))
    
    N = float(dd.shape[0])
    deltaPrime = N * np.sum(np.multiply(dd,dd)) - np.sum(dd)**2
    err = xx - a*dd - b
    sigmaSquare = (1.0/(N-2)) * np.sum(np.multiply(xx,xx))
    
    aSigma = (sigmaSquare/deltaPrime) * np.sum(np.multiply(dd,dd))
    print( 'aSigma:', format(np.sqrt(aSigma), '.4f') )
    print( '(aSigma/a)*100:', format(100*np.sqrt(aSigma)/a, '.3f') )
    
    
    plt.figure(figsize=(15,10))
    plt.plot(dd, xx, 'yo', dd, a*dd+b, '--k')
    plt.show()
In [17]:
print(dda)
ll =len(dda)
xx = np.linspace(1, ll, ll) * 0.001
dd = ( np.asarray(dda)) * 7.4e-6
#dd -= np.mean(dd)
fitAndShow(xx, dd, 2)
#res = stats.theilslopes(xx, dd, 0.90)
#print(res)

[803.5, 812.5, 823.5, 832.5, 842.5, 853.5, 863.5, 872.5, 883.5, 892.5, 903.5, 913.5, 923.5]
a: 13.4670
b: -0.0790
aSigma: 0.0551
(aSigma/a)*100: 0.409
In [18]:
# fitting con campioni pari e dispari

ll = len(dda)//2
xx = np.linspace(1, ll, ll ) * 0.002
dda_d = [dda[2*i+1] for i in range(ll) ]
dd = np.asarray(dda_d) * 7.4e-6
fitAndShow(xx, dd, 2)
#res = stats.theilslopes(xx, dd, 0.90)
#print(res)

ll = len(dda)//2+1
xx = np.linspace(1, ll, ll ) * 0.002
dda_p = [dda[2*i] for i in range(ll) ]
dd = np.asarray(dda_p) * 7.4e-6
fitAndShow(xx, dd, 2)
#res = stats.theilslopes(xx, dd, 0.90)
#print(res)

a: 13.4346
b: -0.0788
aSigma: 0.0979
(aSigma/a)*100: 0.729

a: 13.4884
b: -0.0782
aSigma: 0.0863
(aSigma/a)*100: 0.640
In [19]:
model1a = {}
model1a['grid_size'] = 2048
model1a['wavelengths'] = np.linspace(750*1e-9, 750*1e-9, 1)
model1a['beam_ratio'] = 0.5
#model1a['offset'] = np.linspace(-0.09, -0.08, 3) # -0.12
model1a['offset'] = np.linspace(0.045, 0.09, 6) # -0.12
model1a['aperture_diameter'] = 0.011

#rrt1, r1 = runmodel_w(model1)
r1a, s1a, gs1a = runmodel_off(model1a, 13.5, gaussian_occluded_system)
print(s1a)

occrad_m: 0.0055
h: 0.0035096018400939556

Applying lens
Propagating
13.499981457934847
Done:16.6666672
2048
occrad_m: 0.0055
h: 0.0035096018400939556

Applying lens
Propagating
13.499981457934847
Done:33.3333332
2048
occrad_m: 0.0055
h: 0.0035096018400939556

Applying lens
Propagating
13.499981457934847
Done:50.0000002
2048
occrad_m: 0.0055
h: 0.0035096018400939556

Applying lens
Propagating
13.499981457934847
Done:66.6666672
2048
occrad_m: 0.0055
h: 0.0035096018400939556

Applying lens
Propagating
13.499981457934847
Done:83.3333332
2048
occrad_m: 0.0055
h: 0.0035096018400939556

Applying lens
Propagating
13.499981457934847
Done:100.0000002
2048
[3.25522755859126e-06, 3.906270119455969e-06, 4.557312680320679e-06, 5.208355241185387e-06, 5.859397802050097e-06, 6.5104403629148075e-06]
In [20]:
showResult_offset(model1a, r1a, s1a)
dda = [0] * len(s1a) 
for i in range(len(s1a)):
    dda[i] = computeDiameter(r1a[i])
    print((-model1a['offset'][i])/(dda[i]*s1a[i]))

FN [array([ 16.01062387,  16.16043087, -14.25920371]), array([-13.77585128]), array([ 15.61263987, -14.20060256]), array([-14.05614637]), array([-13.85877255]), array([-13.77587209])]

468.84731478491994
(1036.0042399265542, 1016.5573765744058)
plateauLevel: 0.3428415262230872

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

radius 987

-14.013095140343399
494.94252727784817
(1027.2059915019686, 1022.0757052456503)
plateauLevel: 0.6669937644313396

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

radius 990

-13.9706203120149
468.09956688886103
(1037.4789744573104, 1015.8236547705691)
plateauLevel: 0.33404867416576217
radius 991

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

-13.956523228555376
484.44045862822674
(1029.152749916246, 1021.0287137670102)
plateauLevel: 0.4044141954596211

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

radius 992

-13.942452700200102
492.58888184699157
(1027.006103992192, 1022.2153005163403)
plateauLevel: 0.5318360780655822

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

radius 994

-13.914389736712394
495.75541279972776
(1026.836069522872, 1022.337869761177)
plateauLevel: 0.6065451628160435

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

radius 995

-13.900401894787358
In [21]:
from astropy.io import fits

frames = [ 'extra11454.fit',
'extra10623(saturated_point).fit',
'extra9641(saturated_ghost).fit',
'extra8668(saturated_ghost).fit',
'extra7613(saturated_ghost).fit',
'extra6593.fit',
'extra5582.fit',
'extra4686.fit',
'extra3570.fit',
'extra2573.fit' ]

distances=[ 0,
8.31,
18.13,
27.86,
38.41,
48.61,
58.72,
67.68,
78.84,
88.81
]

cube = [0] * 10
        
fig, ax1 = plt.subplots(1,1, figsize=(15,10))

for ii, fname in enumerate(frames):
    fitsdata = fits.open('C:\\Users\\Fabio\\JupyterNotebooks\\' + fname)
    fitsdata.info()
    image = fitsdata[0].data
    images_d = image.copy()
    plt.figure(figsize=(10,10))
    backgroud_level = np.average(images_d[1800:,1800:])
##imm[images_d[0]<300] = 0
    images_d = images_d - np.ones(images_d.shape) * backgroud_level
    print(backgroud_level)
    ax1.plot(images_d[1076,:])
    #plt.imshow(images_d, cmap='hot')
    cube[ii] = images_d.astype(np.float32)
plt.show()
        

Filename: C:\Users\Fabio\JupyterNotebooks\extra11454.fit
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU      44   (3296, 2472)   int16 (rescales to uint16)   
1172.6968224789916
Filename: C:\Users\Fabio\JupyterNotebooks\extra10623(saturated_point).fit
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU      44   (3296, 2472)   int16 (rescales to uint16)   
1198.6709817449707
Filename: C:\Users\Fabio\JupyterNotebooks\extra9641(saturated_ghost).fit
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU      44   (3296, 2472)   int16 (rescales to uint16)   
1213.9205749060989
Filename: C:\Users\Fabio\JupyterNotebooks\extra8668(saturated_ghost).fit
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU      44   (3296, 2472)   int16 (rescales to uint16)   
1221.5916034027248
Filename: C:\Users\Fabio\JupyterNotebooks\extra7613(saturated_ghost).fit
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU      44   (3296, 2472)   int16 (rescales to uint16)   
1236.3570971002036
Filename: C:\Users\Fabio\JupyterNotebooks\extra6593.fit
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU      44   (3296, 2472)   int16 (rescales to uint16)   
1217.731504249427
Filename: C:\Users\Fabio\JupyterNotebooks\extra5582.fit
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU      44   (3296, 2472)   int16 (rescales to uint16)   
1231.424395610517
Filename: C:\Users\Fabio\JupyterNotebooks\extra4686.fit
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU      44   (3296, 2472)   int16 (rescales to uint16)   
1268.3909045152152
Filename: C:\Users\Fabio\JupyterNotebooks\extra3570.fit
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU      44   (3296, 2472)   int16 (rescales to uint16)   
1431.8868271740514
Filename: C:\Users\Fabio\JupyterNotebooks\extra2573.fit
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU      44   (3296, 2472)   int16 (rescales to uint16)   
1703.8831467245989

<Figure size 720x720 with 0 Axes>
<Figure size 720x720 with 0 Axes>
<Figure size 720x720 with 0 Axes>
<Figure size 720x720 with 0 Axes>
<Figure size 720x720 with 0 Axes>
<Figure size 720x720 with 0 Axes>
<Figure size 720x720 with 0 Axes>
<Figure size 720x720 with 0 Axes>
<Figure size 720x720 with 0 Axes>
<Figure size 720x720 with 0 Axes>
In [22]:
dda = [0] * 10
for i in range(10):
    dda[i] = computeDiameter(cube[i])
print(dda)

324.5277139351465
(1073.0901328932773, 1556.8185524697235)
plateauLevel: 0.5042974206753754

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

radius 657

392.9866988948991
(1074.352217616858, 1556.064054181618)
plateauLevel: 0.5568737932486633
radius 791

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

456.8260164109862
(1076.2046722257896, 1555.010169015846)
plateauLevel: 0.5166882429838491
radius 923

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

520.9401068018652
(1078.2330237953647, 1552.801449984224)
plateauLevel: 0.4783367130326404
radius 1053

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

584.5845408119758
(1079.3570894323332, 1549.165146399281)
plateauLevel: 0.3859844118210377
radius 1193

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

658.6866875338346
(1079.7472783203978, 1550.180409629379)
plateauLevel: 0.4476400787641929
radius 1333

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

724.6485361383913
(1081.8750766048836, 1548.7298088924101)
plateauLevel: 0.437282256795725
radius 1469

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

783.8014830569829
(1083.0567511959714, 1545.8534456360837)
plateauLevel: 0.42438807047633553
radius 1588

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

858.1065763229736
(1084.23993711154, 1542.8568103083817)
plateauLevel: 0.428198002542033
radius 1738

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

922.1369374231391
(1085.916348326783, 1534.0552801497035)
plateauLevel: 0.40697233679498585

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

radius 1869

[656.5, 790.5, 922.5, 1052.5, 1192.5, 1332.5, 1468.5, 1587.5, 1737.5, 1868.5]
In [23]:
ll =len(dda)
xx = np.asarray(distances) * 0.001 # np.linspace(1, ll, ll) * 0.001
dd = ( np.asarray(dda)) * 5.5e-6

print(xx, dd)


#dd -= np.mean(dd)
fitAndShow(xx, dd, 2)
#res = stats.theilslopes(xx, dd, 0.90)
#print(res)

[0.      0.00831 0.01813 0.02786 0.03841 0.04861 0.05872 0.06768 0.07884
 0.08881] [0.00361075 0.00434775 0.00507375 0.00578875 0.00655875 0.00732875
 0.00807675 0.00873125 0.00955625 0.01027675]
a: 13.4307
b: -0.0496
aSigma: 0.0627
(aSigma/a)*100: 0.467
In [24]:
model1 = {}
model1['grid_size'] = 2049
model1['wavelengths'] = np.linspace(500*1e-9, 600*1e-9, 1)
model1['beam_ratio'] = 0.5
model1['offset'] = np.linspace(-0.05, -0.035, 3)
model1['aperture_diameter'] = 1.116

#rrt1, r1 = runmodel_w(model1)
r1, s1, gs1 = runmodel_off(model1, 13.63, occluded_system)
print(s1)

Applying lens
Propagating
13.629999999176015
Done:33.3333332
2049
Applying lens
Propagating
13.629999999176015
Done:66.6666672
2049
Applying lens
Propagating
13.629999999176015
Done:100.0000002
2049
[3.5806525226202346e-06, 3.0435546343492215e-06, 2.5064567460782083e-06]
In [ ]:
showResult_offset(model1, r1, s1, True)

FN [array([-134.20307938, -142.41959445, -155.07911395, -166.15619352,
       -634.41455708, 1395.71202558,  872.32001599,  178.93743918,
        162.292096  ,  148.48000272,   13.73732309]), array([-134.20307982, -145.3866698 , -536.81231927, 1163.09335843,
        697.85601506,  148.4800032 ,   13.76441844]), array([ -60.16000158, -129.23259598,  872.32002287,   64.02348792,
         13.81893106])]
In [ ]:
dd = [0] * len(s1) 
for i in range(len(s1)):
    dd[i] = computeDiameter(r1[i])
    print((-model1['offset'][i])/(dd[i]*s1[i]))

502.20416027656756
(1022.5182068763584, 1022.5182068763584)
plateauLevel: 0.46378754435773323
In [ ]:
model3 = {}
model3['grid_size'] = 2049
model3['wavelengths'] = np.linspace(500*1e-9, 600*1e-9, 1)
model3['beam_ratio'] = 0.5
model3['offset'] = np.linspace(0.05, 0.035, 3)
model3['aperture_diameter'] = 1.116

#rrt1, r1 = runmodel_w(model1)
r3, s3, gs3 = runmodel_off(model3, 13.63, occluded_system)
print(s3)
In [ ]:
showResult_offset(model3, r3, s3, True)
In [ ]:
dd3 = [0] * len(s3) 
for i in range(len(s3)):
    dd3[i] = computeDiameter(r3[i])
    print((-model3['offset'][i])/(dd3[i]*s3[i]))
In [ ]:
dda = [0] * 13
for i in range(13):
    dda[i] = computeDiameter(cube[i])
    
print(dda)
In [ ]:
print(dda)
ll =len(dda)
xx = np.linspace(1, ll, ll) * 0.001
dd = ( np.asarray(dda)) * 7.4e-6
#dd -= np.mean(dd)
fitAndShow(xx, dd, 2)
#res = stats.theilslopes(xx, dd, 0.90)
#print(res)
In [ ]:
# fitting con campioni pari e dispari

ll = len(dda)//2
xx = np.linspace(1, ll, ll ) * 0.002
dda_d = [dda[2*i+1] for i in range(ll) ]
dd = np.asarray(dda_d) * 7.4e-6
fitAndShow(xx, dd, 2)
#res = stats.theilslopes(xx, dd, 0.90)
#print(res)

ll = len(dda)//2+1
xx = np.linspace(1, ll, ll ) * 0.002
dda_p = [dda[2*i] for i in range(ll) ]
dd = np.asarray(dda_p) * 7.4e-6
fitAndShow(xx, dd, 2)
#res = stats.theilslopes(xx, dd, 0.90)
#print(res)
In [ ]:
#model2 = {}
#model2['grid_size'] = 2049
#model2['wavelengths'] = np.linspace(500*1e-9, 800*1e-9, 11)
#model2['beam_ratio'] = 0.5
#model2['offset'] = np.linspace(-0.9, -0.09, 1) # -0.12
#model2['aperture_diameter'] = 0.011
#
#r2, s2, gs2 = runmodel_w(model2, 13.5)
In [ ]:
# showResult_w(model2, r2, s2)