md
set mem
800
md
tr off
md
vb off
md
set mp
2
; multireference alignment and classification based on AP MD
; separates shift from rotation
; uses the following procedures:
; alqmd.mar
; centr.mar
; combat.mar



; number of images
x19=4467
; image size
x20=75
; file code
x22=6


; radius for alignment and mask
x21=35
; first radius for alignment
x26=5
; One can do the classification either on factors or on the raw data
; It is controlled by x27
; x27=1   classification on raw data
; x27=2   classification on factors, their number is given by x23
x27=1
; number of factors in CA S
x23=12
; number of groups
x24=100
; number of iterations of the multireference alignment
x25=6

;  Center all the images
@centr[x19]
../win/ser*****
center/prj****
mares/shiftr001


x91=0
do lb18 x92=1,x19
sd x92,x91,x91,x91,x91
mares/trns001
lb18
sd e
mares/trns001

; create mask for PCA
mo
_9
x20,x20
C
x21

; if classification on raw data find number of points under
; the mask
if(x27.eq.1) then
fs x81,x82,x83
_9
x23=x83*x20*x20
x23=INT(x23)
endif
; copy all the centered images to the stack
do lb21 x31=1,x19
cp
center/prj{****x31}
mares/data001@{****x31}
lb21
cp
mares/data001@
mares/data002@


;
; Iterate mutlireference alignment
do lb8 x81=1,x25
x82=x81+1
; Create SEQ file for PCA
ca si
mares/data002@****
1-x19
_9
SEQ{***x22}
Y
(0.0)
; Do the PCA
if(x27.eq.2) then
ca s
SEQ{***x22}
x23,x22
Y
0
endif
; create subdirectory for the classification results
; and the averages
vm
echo "if(-d cluster{**x81}) rm -r cluster{**x81}"|csh
vm
mkdir cluster{**x81}
;
; classification
cl km
x27
x22
x24
1-x23
(0)
(0)
cluster{**x81}/sel***
cluster{**x81}/cluster

; Get the averages
do lb1 x71=1,x24
as dc
mares/data002@****
cluster{**x81}/sel{***x71}
a
cluster{**x81}/avg{***x71}
cluster{**x81}/var{***x71}
lb1


; Optional
; Put averages in a decent orientation
; centr would use rotational average to center the images
;@centr[x24]
;cluster{**x81}/avg***
;cluster{**x81}/cvg***
;cluster{**x81}/centar
;
; Another possibility is to use CG PH to center the averages
do lb51 x51=1,100
cg ph,x10,x10,x41,x42
cluster{**x81}/avg{***x51}
sh f
cluster{**x81}/avg{***x51}
cluster{**x81}/cvg{***x51}
-x41,-x42
lb51



; Align averages roationaly
ap ra
cluster{**x81}/cvg***
1-x24
x26,x21
(0)
F
cluster{**x81}/haldoc

do lb41 x90=1,x24
ud ic,x90,x41,x42,x43,x44
cluster{**x81}/haldoc
rt sq
cluster{**x81}/cvg{***x41}
cluster{**x81}/rvg{***x41}
x42
(0,0)
lb41
ud ice
cluster{**x81}/haldoc

; end of optional alignment of the averages


; do the multireference alignment
;
;  Rotation first
ap md
cluster{**x81}/rvg***
1-x24
x26,x21
(1)
N
mares/data002@****
1-x19
cluster{**x81}/apmd001

; use crosscorelation to center 2D images.
; crosscorelate each projection only with the reference image
; found by APMD
@alqmd
cluster{**x81}/apmd001
mares/data001@****
mares/data002@****
cluster{**x81}/rvg***
mares/trns{***x81}
mares/trns{***x82}



lb8

en