first commit folder mask_package

; $Id: angle.pro, v 1.0 Aug 1999 e.d. $

;

;+

; NAME:

;	ANGLE

;

; PURPOSE:

;	Compute the position angles of a set of points on a plane with

;	respect to the horizontal axis of a reference frame passing through

;	a fixed origin. The angles are measured counter-clockwise in radians

;	and belong to the interval [0, 2*pi[.

;	The computations are performed in floating-point arithmethic.

;

; CATEGORY:

;	Mathematics.

;

; CALLING SEQUENCE:

;	Result = ANGLE(X0, Y0, X, Y)

;

; INPUTS:

;	X0, Y0:	Couple of scalars, representing coordinates of the origin

;

;	X, Y:	Coordinates of the points for which the position angle

;		must be computed

;

; OUTPUTS:

;	Result:	Array of position angles, with the same size as the input

;		arrays X and Y.

;

; MODIFICATION HISTORY:

; 	Written by:	Emiliano Diolaiti, August 1999.

;-

FUNCTION angle, x0, y0, x, y

	on_error, 2

	a = float(x) - x  &  dx = x - float(x0[0])  &  dy = y - float(y0[0])

	w = where(dx eq 0 and dy ne 0, n)

	if  n ne 0  then  a[w] = !pi/2

	w = where(dx ne 0, n)

	if  n ne 0  then  a[w] = atan(dy[w] / dx[w])

	w = ((dx lt 0) or ((dx eq 0) and (dy lt 0))) and 1B

	a = a + w * !pi

	w = where(a lt 0, n)

	if  n ne 0  then  a[w] = a[w] + 2*!pi

	w = where(a eq 2*!pi, n)

	if  n ne 0  then  a[w] = 0

 	return, a

end

;ROUTINE NAME:

;asong_center_pupil

;

;PURPOSE:

;Computes the coordinate center of the pupil/pupils

;This procedure should be used if a pinhole mask is applied in the pupil plane

;

;CALLING SEQUENCE:

;centres = asong_center_pupil(ima, npupils, box=box, click=click, x=x, y=y, gauss=gauss)

;

;INPUT:

;ima : frame 

;npupils: number of pupils of wich compute the centres. 

;;        When the pyramid is mounted npupils should be set to 4

;

;KEYWORDS

;box = box:   dimention of the subarray in wich compute the centroid. 

;             When a pupil mask is mounted, this value can be set

;             around 9, otherwise it should be set > then the pupil diameter

;click=click: if this keyword is not set, the subarrays are built 

;             around the default values or around the x and y values

;             when defined. Otherwise an interactive window is opened. 

;             The values around wich the subarray of box=box dimension

;             are cut to compute the centroid are defined by clicking 

;             with the mouse.

;x=x, y=y:    If x and y are defined, the subarray of dimension box=box 

;             is cut Around these values

;gauss=gauss: compute the centroid by means of goussian fit of the 

;             frame in the subarray having dimesion=box

; 

;median_subtract= median_subtract: set this keyword to subtract the 

;                                  median of the image 

;

;OUTPUT:

;centres : npupil X 2 array containing the pupil/pupils centres

;

FUNCTION asong_center_pupil, ima, npupils, $

                             median_subtract=median_subtract, box=box, click=click, x=x, y=y, gauss=gauss, threshold = threshold

on_error,2

s=size(ima,/dim)

if keyword_set(median_subtract) then begin

   ima=ima-median(ima)

endif   

  if not keyword_set(box) then box=9

  if keyword_set(click) then begin

      window,xs=s[0],ys=s[1]

      tvscl,congrid(ima,s[0],s[1])

      click_on_max, ima, /mark, nselect=npupils, boxsize=box, x, y

  endif else begin

      if not keyword_set(x) then begin

        x = s[0]/2 

      endif

      if not keyword_set(y) then begin

        y = s[1]/2

      endif  

  endelse       

  centres=fltarr(npupils,2)

  box_cube=fltarr(box,box,npupils)

  for i=0,npupils-1 do begin

     box_cube[*,*,i]=sub_array(ima,box,reference=[x[i],y[i]])

     if keyword_set(gauss) then begin 

        g = gauss2dfit(box_cube[*,*,i],a)

        centres[i,0]=a[4]+x[i]-box/2 &  centres[i,1]=a[5]+y[i]-box/2

     endif else begin

        a = centroid(box_cube[*,*,i])

        centres[i,0]=a[0]+x[i]-box/2 &  centres[i,1]=a[1]+y[i]-box/2

     endelse   

  endfor

return,centres

end

 No newline at end of file

;+ Author: Laura Schreiber - August 2022

;

; ROUTINE NAME:

;asong_check_mask

;

;PURPOSE:

;To verify the mask alignement as size with the 4 measured pupils image.

;Given a mask and an image, it visualize the inverse of the multiplication of the image for the mask and

;for its inverted value.

;

;CALLING SEQUENCE:

;asong_check_mask, image, qmask, xw_size = xw_size, yw_size = yw_size

;

;INPUT:

;image: 4 pupils image taken with ASONG

;qmask: 4 pupils mask to be visualised

;

;KEYWORDS:

;wscale    = wscale    :  Variable type: Scalar. Graphycal window scaling factor w.r.t. image size. DEFAULT VALUE = 0.5    

;;log      = log       : set this keyword to visualise the image logarithm and bust the constrast

;

pro asong_check_mask, image, qmask     , $                         ;input

                      wscale  = wscale , $

                      log     = log    , $

                      _EXTRA  = extra    ;keyword

if not keyword_set(log) then log = 0

if not keyword_set(wscale)  then wscale  =  0.5

xw_size = wscale * (size(image,/dim))[0]

yw_size = wscale * (size(image,/dim))[1]

if log eq 1 then log = 0.1                     

if n_elements(size(qmask,/dim)) eq 3 then $ 

  mask = total(qmask,3)  else $

  mask = qmask

window,1, title = 'mask check direct ', xs = xw_size, ys = yw_size

if not keyword_set(log) then $

   tvscl, congrid(image*mask,xw_size,yw_size) else $

   tvscl, congrid(image^0.1*mask,xw_size,yw_size) 

window,2,title = 'mask check inverse ', xs = xw_size, ys = yw_size

if not keyword_set(log) then $

   tvscl, congrid(image*(1-mask),xw_size,yw_size) else $

   tvscl, congrid(image^log*(1-mask),xw_size,yw_size)

return

end

 No newline at end of file

;+ Author: Laura Schreiber - August 2022

;

; ROUTINE NAME:

; asong_define4quad

;

;PURPOSE:

;For the construction of the qmak.

;

;

;CALLING SEQUENCE:

;asong_define4quad, image

;

;INPUT:

;image: 4 pupils image taken with ASONG

;

;KEYWORDS:

;wscale    = wscale    :  Variable type: Scalar. Graphycal window scaling factor w.r.t. image size. DEFAULT VALUE = 0.5

;silent    = silent    :  set this keyword to not have printed outputs in the IDL console

;log       = log       : set this keyword to visualise the image logarithm and bust the constrast

;

;OUTPUT:

;quad: cube containing the 4 quadrants to select each of the 4 pupils in the correct order

;for the qmask definition.

FUNCTION asong_define4quad, image,             $                                              ;input

                            wscale = wscale  , $

                            silent = silent  , $

                            log = log        , $

                            _EXTRA = extra     ;keyword

quad = fltarr((size(image,/dim))[0],(size(image,/dim))[1],4)

if not keyword_set(wscale)  then wscale  =  0.5

xw_size = wscale * (size(image,/dim))[0]

yw_size = wscale * (size(image,/dim))[1]

window,0, title = 'quad definition',xs = xw_size, ys = yw_size

if not keyword_set(log) then begin

   tvscl, congrid(image,xw_size,yw_size)

endif else begin

   tvscl, congrid(image^0.1,xw_size,yw_size)

endelse

if  not keyword_set(silent)  then begin

  print, 'Select the center of the mask by clicking with the left button of your mouse...'

  print, 'Push right button to exit'

endif

!MOUSE.button = 1 ; &  n_click = 0

;B = ''

;while strupcase(B) ne 'Y' do begin

while  !MOUSE.button eq 1  do begin

;  if  n_click eq 1 then  !MOUSE.button = 4

cursor, x_click, y_click, /NORMAL, /DOWN

  x_click = round(x_click * (size(image,/dim))[0])  &  y_click = round(y_click * (size(image,/dim))[1])

  crosses, tvrd(), x_click/float((size(image,/dim))[0] - 1)*!D.x_size, y_click/float((size(image,/dim))[1] - 1)*!D.y_size, $

    /EXISTING, _EXTRA = extra, /DEVICE, symsize = 200 * wscale, thick = 3

;  read,B,prompt='PRESS Y if THE CROSS IS WELL CENTERED IN THE 4 PUPILS CENTER : '

 ; if  !MOUSE.button eq 1  then  cursor, x_click, y_click, /NORMAL, /DOWN

endwhile

quad[0:x_click,0:y_click,0] = 1

quad[x_click+1:(size(image,/dim))[0]-1,0:y_click,1] = 1

quad[0:x_click,y_click+1:(size(image,/dim))[1]-1,2] = 1

quad[x_click+1:(size(image,/dim))[0]-1,y_click+1:(size(image,/dim))[1]-1,3] = 1

;for i = 0,3 do begin

;  window,i+1, xs = xw_size, ys = yw_size

;  if not keyword_set(log) then begin

;    tvscl, congrid(image * quad[*,*,i],xw_size,yw_size)

;  endif else begin

;    tvscl, congrid(image^0.1 * quad[*,*,i],xw_size,yw_size)

;  endelse  

;endfor

return,quad

end

 No newline at end of file

;+ Author: Laura Schreiber - August 2022

;

; ROUTINE NAME:

;asong_mask

;

;PURPOSE:

;Given a the centres of the four pupils, builds the binary mask for the signal computation and the pupil mask

;

;CALLING SEQUENCE:

;asong_mask, centres, px = px, dim = dim , dimpup = dimpup, rad = rad, qmask, pup, x, y, opup, ox, oy, odimpup

;

;INPUT:

;centres: 4 X 2-elements vector of coordinates of metapupils centroids

;

;KEYWORDS:

;px = px :   oversampling factor for the numerical derivative 

;            computation in lz_derivative. Default value is 1

;dim = dim : The dimension of the frame containing the pupils images.

;            Default value is 3072X2048 

;dimpup = dimpup: The dimension of the array that will contain the 

;                 pupil for the wavefront calculation. It can be 

;                 slightly larger than the calculated pupil. Default

;                 value is 900

;rad = rad :    Measured pupil radius. Default value is 450 pixels     

;

;OUTPUT:

;qmask:  3D array with four planes, each representing a binary mask

;        associated to one quadrant of the sensor

;pup:    pupil mask, logical AND of four masks

;x, y:   coordinates

;opup:   pupil mask defined on a finer grid (defined by the px keyword)

;ox, oy: coordinates in the finer reference

;odimpup:dimension of opup

;-

;MODIFICATION HISTORY:

; 29/09/2022: added +1 to the rad value in the qmask definition to avoid edge issues

PRO asong_mask, centres,                $                                                            ;input

                px        = px        , $

                dimension = dimension , $

                dimpup    = dimpup    , $

                rad       = rad       , $   

                _EXTRA    = EXTRA     , $    ;keywords  

                qmask, pup, x, y, opup, ox, oy, odimpup, c                            ;output

on_error,2

if (size(centres))[0] eq 1 then c = transpose(rebin(centres,2,4)) else $

c=centres

if not keyword_set(rad) then rad=450

if not keyword_set(dimension) then dimension=[3072,2048]

if not keyword_set(dimpup) then dimpup = 2 * rad

if not keyword_set(px) then px=1

qmask=fltarr(dimension[0],dimension[1],4)

y=rebin(transpose(findgen(dimension[1])),dimension[0],dimension[1])

x=rebin(findgen(dimension[0]),dimension[0],dimension[1])

for i=0,3 do begin

   qmask[*,*,i]=sqrt((x-c[i,0])^2+(y-c[i,1])^2) le rad + 1

endfor

y=rebin(transpose(findgen(dimpup)),dimpup,dimpup)-dimpup/2

x=rebin(findgen(dimpup),dimpup,dimpup)-dimpup/2

pup=sqrt(x^2+y^2) le rad

odimpup = px * dimpup

oy=rebin(transpose(findgen(odimpup)),odimpup,odimpup)-odimpup/2

ox=rebin(findgen(odimpup),odimpup,odimpup)-odimpup/2

d = odimpup + 2

oy = rebin(transpose(findgen(d)),d,d)-d/2

ox = rebin(findgen(d),d,d)-d/2

opup = sqrt(ox^2+oy^2) le rad * px

return

end

 No newline at end of file