; 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.spi 
; 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  ; Width of the circular line, in microns

x99 = 3.72  ; Pixelsize

x80 = 1     ; 1 = save intermediate files

; ----- Input files -----

[mic] = './mic{***x85}'             ; Full-sized micrograph

[filenums] = '../filenums'

; ----- Output files -----

[outdir] = 'cmask'                  ; Output directory

[circle] = '[outdir]/circ{***x85}'  ; Image of the circle

[small]  = '[outdir]/small{***x85}' ; Small version of micrograph

[ring]   = '[outdir]/sup{***x85}'     ; Circle superimposed on small micrograph
                                    ; CHECK THIS FILE TO ENSURE CORRECT MASK

[outdoc] = 'coords'                 ; 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  input image  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