; cmask
;
; Create a circular mask for the round Quantifoil hole in a micrograph.
; This assumes the hole is surrounded by a slightly brighter ring (inverted 
; contrast).
;
; This batch file creates small intermediate files that should be checked
; to make sure the circle is located correctly. The main output is a 
; coordinate doc file that is read by the batch file cmfull.bat, which
; applies the mask to the full-sized micrograph.

;
; NB : while Quantafoil guarantees consistent hole-to-hole spacing, the
; diameter of the ring and its width vary from one batch to another. Therefore,
; you may need to adjust ring parameters x97 and x98


; In addition to normalization, it uses an inner negative ring.
;
; Inputs: 
;     filenums : file numbers
;     micrographs in spider format
;
; Outputs
;     coords : doc file of alignment parameters (to be used in cmfull.bat)
;              columns: X_center, Y_center, radius
;     The following files are used to check correctness. 
;     set x80 = 1 to save these files (default), otherwise they are not saved.
;        circle : image of circle
;        small  : tiny image of micrograph
;        ring   : circle superimposed on little micrograph.
;
; Check the ring images, if ok, go to cmfull.bat. 
; If not ok, adjust RING PARAMETERS

MD
SET MP
0

x40 = 10 ;  reduction factor

; ===== RING PARAMETERS =====

x97 = 2.65  ; diameter of circle, in microns
x98 = 0.16  ; the width of the circular line, in microns

x99 = 3.72  ; pixelsize

x80 = 1     ; 1 = save intermediate files

; ----- input files -----

FR G
[mic]./mic{***x85}     ; full-sized micrograph

FR G
[filenums]../filenums

; ----- output files -----

FR G
[outdir]cmask       ; output directory

FR G
[circle][outdir]/circ{***x85}   ; image of the circle

FR G
[small][outdir]/small{***x85}    ; small version of micrograph

FR G
[ring][outdir]/sup{***x85}     ; circle superimposed on small micrograph
                               ; CHECK THIS FILE TO ENSURE CORRECT MASK

FR G
[outdoc]coords    ; output docfile of coordinates: mic x y radius

; ------------ END BATCH HEADER ------------------------------------

DE                 ; deletes the output document at start
[outdoc]

VM                      ; Make sure output dir. is present
mkdir -p [outdir]


x99 = x99 * x40  ; new pixelsize

UD N,x90
[filenums]

SD /       MIC          X COORD      Y COORD      RADIUS
[outdoc]

DO LB1 x84 = 1,x90

UD x84,x85      ; x85 is now the file number
[filenums]

VM
echo working on [mic]

; shrink the micrograph
FI x21,x22
[mic]
(12,2)

x31 = INT(x21/x40)
x32 = INT(x22/x40)

IP
[mic]
_1
(x31,x32)

IF (x80.EQ.1) THEN
   CP
   _1
   [small]
ENDIF

FI x11,x12    ; get input image size
_1
(12,2)

; convert to pixels
x20 = 10000 * x97/x99   ; diameter of circle
x21 = 10000 * x98/x99   ; width of line
x23 = INT(x20 + x21) + 1         ; width, height of circle image

x24 = x23 - x11   ; pad x
x25 = x23 - x12   ; pad y

IF (x24.LT.0) x24 = 0
IF (x25.LT.0) x25 = 0

PD
_1    ; pad the input image to be as big as the circle
_2    ; padded image
(x11+x24, x12+x25)
Y
(1,1)

CP 
_2
_1

; make the circle masks

x19 = x20/2.0   ; radius
x19
x23

MO
_2        ; outer circle
(x23,x23)
C
(x19)

x19 = (x20-x21)/2.0   ; radius of inner circle

MO
_3        ; inner circle
(x23,x23)
C
(x19)

SU
_2
_3
_5   ; subtract small from large circle to make ring 
*

CP
_5
_44   ; save regular ring for later, without inner negative ring


AR IF    ; change small circle, to create inner negative ring
_3
_4
IF (P1.eq.1) THEN P2 = 2




CP
_4
_3

SU
_2
_3
_5   ; subtract small from large circle to make ring (-1 in center)
*

x69 = (x20 - (2*x21)) / 2.0
x69

MO  
_61   ; inner inner circle
(x23,x23)
C
(x69)

AD
_5
_61 
_4   ; subtract even smaller  from ring to give negative inner ring
*

IF (x80.EQ.1) THEN
   CP
   _44
   [circle]    ; only positive circle, without negative inner ring
ENDIF

; ********** Cross-correlation ***************
CC N
_1     ; mic
_4     ; ring
_5

;CP
;_5
;cca{***x85}

; "normalize" the CC (sums) by the # pixels under the ring at each point

BL
_33        
(x23,x23)     ; blank image size of padded img
N
(0.0)

PT
_33
R
(1,1)
(x11,x12)    ;  draw rectangle the size of unpadded img
N

CC 
_44    ; ring
_33   ; rectangle image
_34   ; count of pixels when circle overlaps rectangle

RT 90
_34    ; input
_35    ; output




(180)  ; rotation angle for RT 90


MU D
_5         ; cc result
_35        ; divided by pixel counts
_36
*

CP
_36
_5

;CP
;_5
;ccb{***x85}

; --- end normalization
; get X,Y shifts from peak search
PK x41,x42
_5           ; cc result
(1,0)
(1)

SH
_44    ; regular ring, without inner negative ring  
_8
x41,x42

AR IF
_8
_9
IF (p1.GT.0) THEN p2 = 1.01
ELSE IF(P1.EQ.0) THEN p2=1.0

IF (x80.EQ.1) THEN
   MU
   _9
   _1
   [ring]
   *
ENDIF


; use the inner circle as a mask 
; (create a new one. Can't use SH since it does a circular shift)
BL
_3        
(x11,x12)
N
(0.0)

x25 = x23/2.0 + 1
;x26 = x12/2.0 + 1

x71 = x25+x41    ; x,y coords of circle center
x72 = x25+x42

PT
_3
C
(x25+x41, x25+x42)
(x19)       ; radius of the inner circle
N

SD x84,x85,x71,x72,x19
[outdoc]


LB1

EN D