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

import cv2 as cv
In [2]:
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 [3]:
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 showResult_offset(model, r1, s1):
    gridsize = r1[0].shape[0]
    gridsizeF = float(r1[0].shape[0])
    
    amps = []
    phs = []
    for i in range(len(r1)):
        amp = np.abs(r1[i])
        max_amp = np.max(amp)
        ph = np.arctan2(np.imag(r1[i]), np.real(r1[i]))
        #ph = np.where(amp < 0.01*max_amp, 0.0, ph)
        amps.append(amp)
        phs.append(ph)
    
    fig, ( (ax11, ax12), (ax13, ax14),(ax21, ax22), (ax23, ax24), (ax31, ax32), (ax33, ax34) ) = plt.subplots(6,2, figsize=(10, 30))
    ax11.matshow( amps[0] )
    ax12.matshow( phs[0] )
    ax21.matshow( amps[1] )
    ax22.matshow( phs[1] )
    ax31.matshow( amps[2] )
    ax32.matshow( phs[2] )
    FN = [0] * len(r1)
    axplots = [ax13, ax23, ax33]
    axplotsP = [ax14, ax24, ax34]

    for i in range(len(r1)):
        amp = np.abs(r1[i])
        max_amp = np.max(amp)
        ph = np.arctan2(np.imag(r1[i]), np.real(r1[i]))
        ph = np.where(amp < 0.01*max_amp, 0.0, ph)
        section = (amp[:, gridsize//2])
        sectionP = (ph[:, gridsize//2])
        x_coords = np.array(np.linspace( 0, s1[i]*gridsizeF, gridsize)) - s1[i]*gridsizeF/2.0 
        #axplots[i].set_yticks(np.arange(0, 1., 0.05))
        axplots[i].plot(x_coords, section)
        #axplotsP[i].set_yticks(np.arange(-2*np.pi, 2*np.pi, 0.3))
        axplotsP[i].plot(x_coords, sectionP)
        axplots[i].grid()
        axplotsP[i].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 [5]:
from skimage.transform import resize
from skimage import data, img_as_float, exposure
from skimage.draw import circle, rectangle

def lensletArrayScreen( gridsize, sampling, num_lens_per_side, FN, mm ):
    wv = 600e-9
    k = 2.0*np.pi/wv
    print('sampling:', sampling)
    si = sampling 
    print('si:', si)
    lens_diam = gridsize // num_lens_per_side
    print('lens_diam:', lens_diam)
    radiusInPixels = lens_diam // 2
    M = radiusInPixels
    print('radiusInPixels:', radiusInPixels)
    complex_image = np.ones((gridsize, gridsize)).astype(np.complex)
    (x1,y1) = np.meshgrid(si*np.arange(-M,M), si*np.arange(-M,M))
    r1sq = (x1**2 + y1**2)
    focale = FN * mm * lens_diam * si 
    phi = k*(- r1sq/(2.0*focale) )

    print('np.max(phi):', np.max(phi))
    print('np.min(phi):', np.min(phi))

#    phi1 = np.exp( -1j*2*np.pi*( np.floor(phi/(2*np.pi)) ) )
#    phi2 = np.exp( -1j*np.mod(phi, 2*np.pi) )
    phi = np.where(r1sq < si*si*M*M, np.exp( 1j*phi), 0.0)
    for i in range(int(num_lens_per_side)):
        for j in range(int(num_lens_per_side)):
            rr, cc = circle(i*lens_diam+radiusInPixels, j*lens_diam+radiusInPixels, radiusInPixels)
            complex_image[rr, cc] =  1.0 # + np.exp(-1j*k*si*M/(2.0*FN))
            hd = 0 # radiusInPixels # int(lens_diam//2)
            dd = lens_diam+hd - phi.shape[0] 
            complex_image[i*lens_diam+hd:(i+1)*lens_diam+hd-dd, j*lens_diam+hd:(j+1)*lens_diam+hd-dd] = phi
    return complex_image

def lensletArrayScreenPhase( gridsize, sampling, num_lens_per_side, FN ):
    wv = 600e-9
    k = 2.0*np.pi/wv
    print('sampling:', sampling)
    si = sampling # / float(gridsize)
    print('si:', si)
    lens_diam = gridsize // num_lens_per_side
    print('lens_diam:', lens_diam)
    radiusInPixels = lens_diam // 2
    M = radiusInPixels
    print('radiusInPixels:', radiusInPixels)
    complex_image = np.ones((gridsize, gridsize)).astype(np.float)
    (x1,y1) = np.meshgrid(si*np.arange(-M,M), si*np.arange(-M,M))
    r1sq = (x1**2 + y1**2)
    focale = FN * lens_diam * si
    phi = -k*r1sq/(2.0*focale)
    # phi = np.exp(1j*phi)
    for i in range(int(num_lens_per_side)):
        for j in range(int(num_lens_per_side)):
            rr, cc = circle(i*lens_diam+radiusInPixels, j*lens_diam+radiusInPixels, radiusInPixels)
            complex_image[rr, cc] =  1.0 # + np.exp(-1j*k*si*M/(2.0*FN))
            hd = 0 # radiusInPixels # int(lens_diam//2)
            dd = 0
            complex_image[i*lens_diam+hd:(i+1)*lens_diam+hd-dd, j*lens_diam+hd:(j+1)*lens_diam+hd-dd] = phi
    return complex_image

def complexDisplay(image):
    f, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize = (5, 15))
    imageAmp = np.abs(image)
    ax1.imshow(imageAmp, cmap = 'hot')
    ax2.imshow(np.arctan2(image.imag, image.real), cmap = 'hot')
    N = image.shape[0]
    xx = 1.0* np.arange(-N/2,N/2)
    ax3.plot(xx, imageAmp[N//2,:])
    plt.show()    
    
def simple_telescope(wavelength, gridsize, d_objective, offset, phaseOn=True, period = 2.0, dist_amp = 1.0, lenslet=False):
    d_collimator = (5.0/17.7) * (d_objective/40.0)
    fl_objective = 17.7 * d_objective # objective focal length in meters
    fl_collimator = 17.7 * d_collimator # eyepiece focal length
    beam_ratio = 1.0 # 0.75 # initial beam width/grid width

    wfos = []
    samplings = []
    ff = 1.3
    
    wfo = proper.prop_begin((d_objective)*ff, wavelength, gridsize, beam_ratio)    
    proper.prop_circular_aperture(wfo, (d_objective/2.0)*ff)
    proper.prop_define_entrance(wfo)

    if phaseOn:
        m = 2.0 * np.pi * dist_amp
        x = (np.arange(gridsize, dtype = np.float64) - gridsize/2) * proper.prop_get_sampling(wfo)
        grating =  m * (np.cos((2*np.pi/period)*x)+1) * 0.5
        # create 2-D amplitude grating pattern
        grating = np.dot(grating.reshape(gridsize,1), np.ones([1,gridsize], dtype = np.float64))
        grating_pahse = np.exp(1j*grating)
        proper.prop_multiply(wfo, grating_pahse)

    #wfos.append(proper.prop_get_wavefront(wfo))
    #samplings.append(proper.prop_get_sampling(wfo))
    
    proper.prop_propagate(wfo, 10e3, "from atmosphere layer", TO_PLANE=False)
    
    wfos.append(proper.prop_get_wavefront(wfo))
    samplings.append(proper.prop_get_sampling(wfo))

    proper.prop_circular_aperture(wfo, d_objective/2)
    proper.prop_lens(wfo, fl_objective, "objective")
    proper.prop_propagate(wfo, fl_objective+fl_collimator, "collimator", TO_PLANE=False)

    wfos.append(proper.prop_get_wavefront(wfo))
    samplings.append(proper.prop_get_sampling(wfo))
    
    proper.prop_circular_aperture(wfo, d_collimator/2)
    proper.prop_lens(wfo, fl_collimator, "collimator")

    #exit_pupil_distance = fl_collimator / (1.0 - fl_collimator/(fl_objective+fl_collimator))
    exit_pupil_distance = fl_collimator - 10e3* (fl_collimator/fl_objective) ** 2
    print(exit_pupil_distance) # layer distance
    print('exit_pupil_distance:', exit_pupil_distance-offset)

    proper.prop_propagate(wfo, exit_pupil_distance-offset, "exit pupil", TO_PLANE=True)

    if lenslet:
        n_lens = 10
        FN_lenslet_array = 3.0 # 40.0 
        FN_one_lenslet = FN_lenslet_array * n_lens
        #lenslet_screen_phase = lensletArrayScreenPhase(gridsize, proper.prop_get_sampling(wfo), n_lens*2, FN_lenslet)
        #proper.prop_add_phase(wfo, lenslet_screen_phase)

        mm = 1
        
        lenslet_screen = lensletArrayScreen(gridsize, proper.prop_get_sampling(wfo), n_lens*1.3 + 1, FN_one_lenslet, mm)
        complexDisplay(lenslet_screen)
        proper.prop_multiply(wfo, lenslet_screen)
        
        proper.prop_propagate(wfo, mm * FN_lenslet_array * d_collimator, "exit pupil", TO_PLANE=True)    

    
    wfos.append(proper.prop_get_wavefront(wfo))
    samplings.append(proper.prop_get_sampling(wfo))

    (wfo, sampling) = proper.prop_end(wfo, NOABS = True)    

    return (wfos, samplings)
In [6]:
def runmodel_mono(model, system=simple_telescope):
    gridsize = model['grid_size'] 
    wavelengths = model['wavelengths'] 
    beam_ratio = model['beam_ratio']
    pahse_on = model['pahse_on']
    r1, s1 = system(wavelengths[0], gridsize, model['aperture_diameter'], model['offset'][0], pahse_on,  model['dist_period'], model['dist_amp'], model['lenslet'])
    return r1, s1
In [7]:
model0a = {}
model0a['grid_size'] = 1024*16
model0a['wavelengths'] = np.linspace(600*1e-9, 600*1e-9, 1)
model0a['beam_ratio'] = 0.5
model0a['offset'] = np.linspace(0, 0, 1)
model0a['aperture_diameter'] = 40.0
model0a['pahse_on'] = True
model0a['dist_period'] = 1.0
model0a['dist_amp'] = 0.05
model0a['lenslet'] = True

#rrt1, r1 = runmodel_w(model1)
r0a, s0a = runmodel_mono(model0a)
print(s0a)

Propagating to from atmosphere layer
Applying lens at objective
Propagating to collimator
Applying lens at collimator
4.501260812665581
exit_pupil_distance: 4.50126081267
Propagating to exit pupil
sampling: 2.2414040429115756e-05
si: 2.2414040429115756e-05
lens_diam: 1170.0
radiusInPixels: 585.0
np.max(phi): -0.0
np.min(phi): -2288.51300435

/home/frossi/anaconda3/envs/cuda/lib/python3.6/site-packages/ipykernel_launcher.py:34: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future

Propagating to exit pupil
[0.003173828125, 2.2414040429115756e-05, 2.2414040429115756e-05]
In [8]:
showResult_offset(model0a, r0a, s0a)
In [9]:
model0a = {}
model0a['grid_size'] = 1024*8
model0a['wavelengths'] = np.linspace(600*1e-9, 600*1e-9, 1)
model0a['beam_ratio'] = 0.5
model0a['offset'] = np.linspace(0, 0, 1)
model0a['aperture_diameter'] = 40.0
model0a['pahse_on'] = True
model0a['dist_period'] = 1.0
model0a['dist_amp'] = 0.05
model0a['lenslet'] = True

#rrt1, r1 = runmodel_w(model1)
r0a, s0a = runmodel_mono(model0a)
print(s0a)

Propagating to from atmosphere layer
Applying lens at objective
Propagating to collimator
Applying lens at collimator
4.501260812665581
exit_pupil_distance: 4.50126081267
Propagating to exit pupil
sampling: 4.482808085823151e-05
si: 4.482808085823151e-05
lens_diam: 585.0
radiusInPixels: 292.0
np.max(phi): -0.0
np.min(phi): -2280.69571543

/home/frossi/anaconda3/envs/cuda/lib/python3.6/site-packages/ipykernel_launcher.py:34: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future

Propagating to exit pupil
[0.00634765625, 4.482808085823151e-05, 4.482808085823151e-05]
In [10]:
showResult_offset(model0a, r0a, s0a)
In [11]:
model0a = {}
model0a['grid_size'] = 1024*2
model0a['wavelengths'] = np.linspace(600*1e-9, 600*1e-9, 1)
model0a['beam_ratio'] = 0.5
model0a['offset'] = np.linspace(0, 0, 1)
model0a['aperture_diameter'] = 40.0
model0a['pahse_on'] = True
model0a['dist_period'] = 1.0
model0a['dist_amp'] = 0.05
model0a['lenslet'] = True

#rrt1, r1 = runmodel_w(model1)
r0a, s0a = runmodel_mono(model0a)
print(s0a)

Propagating to from atmosphere layer
Applying lens at objective
Propagating to collimator
Applying lens at collimator
4.501260812665581
exit_pupil_distance: 4.50126081267
Propagating to exit pupil
sampling: 0.00017931232343292605
si: 0.00017931232343292605
lens_diam: 146.0
radiusInPixels: 73.0
np.max(phi): -0.0
np.min(phi): -2284.60101631

/home/frossi/anaconda3/envs/cuda/lib/python3.6/site-packages/ipykernel_launcher.py:34: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future

Propagating to exit pupil
[0.025390625, 0.00017931232343292605, 0.00017931232343292605]
In [12]:
showResult_offset(model0a, r0a, s0a)
In [13]:
model3a = {}
model3a['grid_size'] = 1024*4
model3a['wavelengths'] = np.linspace(600*1e-9, 600*1e-9, 1)
model3a['beam_ratio'] = 0.5
model3a['offset'] = np.linspace(2.5, 2.5, 1)
model3a['aperture_diameter'] = 40.0
model3a['pahse_on'] = True
model3a['dist_period'] = 1.0
model3a['dist_amp'] = 1.0
model3a['lenslet'] = False
r3a, s3a = runmodel_mono(model3a)
print(s3a)

Propagating to from atmosphere layer
Applying lens at objective
Propagating to collimator
Applying lens at collimator
4.501260812665581
exit_pupil_distance: 2.00126081267
Propagating to exit pupil
[0.0126953125, 8.965616171646303e-05, 8.965616171646303e-05]
In [14]:
showResult_offset(model3a, r3a, s3a)
In [15]:
model2a = {}
model2a['grid_size'] = 1024*4
model2a['wavelengths'] = np.linspace(600*1e-9, 600*1e-9, 1)
model2a['beam_ratio'] = 0.5
model2a['offset'] = np.linspace(0, 0, 1)
model2a['aperture_diameter'] = 40.0
model2a['pahse_on'] = True
model2a['dist_period'] = 1.0
model2a['dist_amp'] = 1.0
model2a['lenslet'] = False

#rrt1, r1 = runmodel_w(model1)
r2a, s2a = runmodel_mono(model2a)
print(s2a)

Propagating to from atmosphere layer
Applying lens at objective
Propagating to collimator
Applying lens at collimator
4.501260812665581
exit_pupil_distance: 4.50126081267
Propagating to exit pupil
[0.0126953125, 8.965616171646303e-05, 8.965616171646303e-05]
In [16]:
showResult_offset(model2a, r2a, s2a)
In [17]:
model1a = {}
model1a['grid_size'] = 1024*4
model1a['wavelengths'] = np.linspace(600*1e-9, 600*1e-9, 1)
model1a['beam_ratio'] = 0.5
model1a['offset'] = np.linspace(0, 0, 1)
model1a['aperture_diameter'] = 40.0
model1a['pahse_on'] = True
model1a['dist_period'] = 1.0
model1a['dist_amp'] = 0.2
model1a['lenslet'] = False

#rrt1, r1 = runmodel_w(model1)
r1a, s1a = runmodel_mono(model1a)
print(s1a)

Propagating to from atmosphere layer
Applying lens at objective
Propagating to collimator
Applying lens at collimator
4.501260812665581
exit_pupil_distance: 4.50126081267
Propagating to exit pupil
[0.0126953125, 8.965616171646303e-05, 8.965616171646303e-05]
In [18]:
showResult_offset(model1a, r1a, s1a)
In [ ]:
In [ ]:
In [ ]: