; CLASSIFY PARTICLES ASSIGNED TO EACH REFERENCE VIEW ; ----------------- Parameters ----------------- x81 = 1 ; first reference-view x83 = -1 ; last reference-view (<1==all) x13 = 1 ; correspondence analysis (1), PCA (2), IPCA (3) x75 = 75 ; particle-to-class ratio x40 = 40 ; minimum number of particles for 2 classes x91 = 0 ; flag to save eigenimages (1==save) reconstituted (images commented out) x97 = 9 ; number of eigenvalues to use x26 = 2 ; reduction factor applied to input particles x28 = 115 ; window size (temporary) for labeling ; (Class number and size unlikely to fit in label after size-reduction.) ; Suggestions: 115 for 3-digit class# + 3-digit class-size, 128 for 3+4 or 4+3 x39 = 0 ; additive constant for correspondence analysis (0==automatic) ; (CorAn only works on non-negative pixel values.) ; if unsure, check the range of a few particles or the noise-reference for positivity x76 = 0 ; flag to save variance (1==save) ; ------------------- Inputs ------------------- fr l [parameter_doc]../params ; parameter document fr l [ref_view_list]../Alignment/projlist ; list of projection views fr l [df_group_list]../Power_Spectra/order_defgrps ; (optional) defocus-group list fr g [ref_proj]../Alignment/projs/prj_{***x19}@ ; reference-projection prefix fr g [view_dir]select/prj{***x80} ; I/O subdirectory pattern, for each reference-view fr g [selection_doc]select/sel{***x80} ; selection doc filename prefix ; VIEW_SLICE PARTICLE DF_SLICE CC_ROT MIRROR GROUPNUM VIEW fr g [particles][view_dir]/stkfilt@****** ; particle template ; ------------------- Outputs ------------------- fr g [class_doc][view_dir]/docclass ; prefix for list of particles for each class fr g [class_stats_doc][view_dir]/listclasses ; contains CCC for each class-average fr g [class_avg][view_dir]/classavg{***x16} ; class-average file pattern ;[class_avg]select/classavg{***x80} ; class-average (one class per view) fr g [class_var][view_dir]/classvar*** ; variance for each class (optional, see X76) fr g [cas_prefix][view_dir]/verify ; CA S output prefix fr l [eigenvalue_doc][view_dir]/listeigenvalues ; list of eigenvalues fr l [eigen_img][view_dir]/eigenimg*** ; eigenimage template ;fr l ;[reconstituted_img][view_dir]/reconstituted*** ; reconstituted image template ;; lines for reconstituted images commented out fr l [summary_doc]summary-classify ; summary doc file ; -------------- END BATCH HEADER -------------- fr g [class_map][view_dir]/classmap ; temporary list of each particle with assigned class ; check if defocus-group list exists iq fi x14 [df_group_list] ; if exists, then if(x14.eq.1) then ; get #defocus-groups (x18) ud n,x18 [df_group_list] ; get last defocus-group# (x19) ud x18,x19 [df_group_list] else x19=1 endif vm echo "Using projections from defocus group #{***x19} for CCC" ; echo ; get image dimension (x65) ud 17,x65 [parameter_doc] x65=x65/x26 ; x65==reduced-image dimension, x26==reduction factor ; calculate mask radius x32=int(x65/2) ; x65==reduced image dimension ; make circular mask mo _6 ; output mask file x65,x65 ; dimensions c ; _C_ircle x32 ; radius if(x83.lt.1) then ; get number of views ud n,x83 [ref_view_list] endif ; initialize particle-counter [total-particles] = 0 ; label summary file SD / VIEW_NUM NUMPARTS MAX_CLASS_SIZE [summary_doc] vm echo "Beginning classification"; date md set mp 0 md vb off if(x91.eq.1) then vm echo ; echo "Will generate eigenimages" ; echo endif ; loop through reference views do lb1 x48=x81,x83 ; get reference-view# (x80) ud ic,x48,x80 [ref_view_list] ; get number of particles (x27) in current reference view ud n,x27 [selection_doc] ; initialize maximum class size x68 = 0 ; skip unpopulated views if(x27.eq.0) then x20 = 0 ; #classes x68 = 0 ; max.class-size goto lb19 endif ; clean up de a [class_doc]001 x16=1 de a [class_avg] de a [class_var]x16 de [eigenvalue_doc] ; shrink reference-projection dc s [ref_proj]{****x80} ; INPUT: normalized reference-projection _5 ; OUTPUT: reduced reference-projection x26,x26 ; reduction factor in x,y ; determine number of classes to use (x20) x20=int((x27/x75)+0.5) ; x27==#particles, x75==ratio ; trap for tiny classes if(x20.lt.2) then if(x27.ge.x40) then ; force two classes if greater than specified threshold x20=2 else x20=1 ; minimum number of classes will be 1 x16=1 ; dummy variable: class# ; sort particles by cross-correlation coefficient doc ren ; WAS doc key [selection_doc] ; INPUT: selection doc (unsorted) [class_doc]{***x16}_unsort ; OUTPUT: class doc, sorted by CCC vm echo "Reference view: #{***x80}, number of classes: {***x20}" x69 = x27 ; #particles in class goto lb3 ; skip past correspondence analysis + classification endif endif ; if CA selected: if(x13.eq.1) then ; CHECK FOR SMALL PIXEL INTENSITIES x79=9999 ; initialize minimum ; loop through particles do lb5 x50=1,x27 ; get particle number (x52) ud ic,x50,x52 [selection_doc] ; get max (x88) fs [particles]x52 ; WAS x50 ; WAS {******x52} fi x88 [particles]x52 ; WAS x50 ; WAS {******x52} (8) ; header location for minimum intensity ; update if necessary if(x88.lt.x79) then x79=x88 ; minumum intensity x49=x50 ; slice # ; WAS x52 ; particle# endif lb5 ; end particle-loop ud ice ; close document [selection_doc] ; if additive constant set to automatic, set it if (x79.lt.0.05) then if(x39.eq.0) x29 = 0.05-x79 else x29 = x39 endif vm echo "Reference view: #{***x80}, minimum intensity: {%f9.4%x79} (image #{******x49}), additive constant: {%f6.3%x29}" endif ; RUN MULTIVARIATE STATISICAL ANALYSIS vm echo "Reference view: #{***x80}, number of classes: {***x20}" if(x13.eq.1) then ; run correspondence analysis ca s [particles] [selection_doc] ; WAS (1-x27) ; WAS _6 ; mask x97 ; number of eigenvalues C ; _C_orrespondence analysis x29 [cas_prefix] ; output file prefix endif if(x13.eq.2) then ; run iterative principle component analysis ca s [particles] [selection_doc] ; WAS (1-x27) ; WAS _6 ; mask x97 ; number of eigenvalues P ; _Principle component analysis [cas_prefix] ; output file prefix endif if(x13.eq.3) then ; run principle component analysis ca s [particles] [selection_doc] ; WAS (1-x27) ; WAS _6 ; mask x97 ; number of eigenvalues I ; _I_terative PCA [cas_prefix] ; output file prefix endif ; run k-means classification cl km x21,x22,x23,x24,x25 [cas_prefix]_IMC x20 ; number of classes (1-9) ; factors to use (0) ; no factor weighting (0) ; no random seed [class_doc]***_noccc ; OUTPUT (temp): class-list doc [class_map] ; if eigenimage-flag is 1, then save if(x91.ge.1) then ; GENERATE EIGENIMAGES (OPTIONALLY) ; ; calculate dimensions for reconstituted images ; [double-idim] = x65*2 ; double image-dimension ; [idim-plus1] = x65+1 ; image-dimension + 1 ; loop through eigenvalues do x90 = 1,x97 ; if (I)PCA, extra question asked when generating eigenimages if(x13.ne.1) then ca sre [cas_prefix] N ; subtract average? x90 ; eigenvector [eigen_img]x90 else ; correspondence analysis ca sre [cas_prefix] x90 ; eigenvector [eigen_img]x90 endif ; ; generate positive reconstituted image ; ca sra ; [cas_prefix] ; x90 ; eigenvector ; (0.2) ; eigenvalue ; _91 ; OUTPUT ; ; ; pad image to twice the height ; pd ; _91 ; INPUT: positive reconstituted image ; [reconstituted_img]x90 ; x65,[double-idim] ; xdim, ydim ; B ; set background to _B_order ; (1,1) ; top-left coordinates ; ; ; generate negative reconstituted image ; ca sra ; [cas_prefix] ; x90 ; eigenvector ; (-0.2) ; eigenvalue ; _93 ; OUTPUT ; ; ; insert negative reconstituted image into larger image ; in ; _93 ; INPUT: negative reconstituted image ; [reconstituted_img]x90 ; INPUT (overwritten) ; (1,[idim-plus1]) ; x-,y-coords sd x90,x90,x90 [eigenvalue_doc] enddo ; end eigenvalue-loop sd e [eigenvalue_doc] endif ; GENERATE CLASS AVERAGES ; loop through classes do lb2 x16=1,x20 ; ADD PARTICLE#, CCC, ETC. TO CLASS DOC ; add info from selection doc doc and [selection_doc] [class_doc]{***x16}_noccc [class_doc]{***x16}_unsort (1) ; column# to intersect: view-slice# lb3 ; skip to here if only one class ; sort individual particles by CCC doc sort [class_doc]{***x16}_unsort ; INPUT (temp): w/CCC, unsorted [class_doc]{***x16}_noccc ; OUTPUT: class doc, sorted by CCROT (4) ; column# to sort: CCROT Y ; renumber? ; get #particles in class ud n,x69 [class_doc]{***x16}_noccc ; if size greater than maximum, then save if(x69.gt.x68) x68=x69 x34=0 ; cumulative CCC-sum ; loop through particles do lb6 x60=1,x69 ; read view-slice#, other parameters, from selection file ud ic x60,[slice-num],[part-num],[df-slice],[ccrot],[mirror],[defocus-group],[view-num] [class_doc]{***x16}_noccc ; WAS [selection_doc] ; calculate CCC cc c [ccc] _5 ; INPUT: reduced reference-projection [particles][slice-num] ; INPUT: particle _6 ; INPUT: reduced mask x34 = x34 + [ccc] ; write to doc sd x60,[slice-num],[part-num],[ccc],[ccrot],[mirror],[defocus-group],[df-slice] [class_doc]{***x16} lb6 ; clean up ud ice [class_doc]{***x16}_noccc sd e [class_doc]{***x16} SD / SLICE PARTICLE CCC CCROT MIRROR DF_GROUP DF_SLICE [class_doc]{***x16} sd e [class_doc]{***x16} ; calculate unlabeled average as r [particles] [class_doc]{***x16} ; INPUT: class-list doc A ; _A_ll images _4 ; OUTPUT: class average _7 ; OUTPUT: class variance ; LABEL AVERAGES ; expand to fit text label, if necessary ip _4 ; INPUT: class-average _2 ; OUTPUT: expanded class-average x28,x28 ; get class size ud n,x15 [class_doc]{***x16} [total-particles] = [total-particles] + x15 ; label with class number & size la b _2 ; INPUT: expanded class-average _3 ; OUTPUT: expanded, labeled class-average {***x16},n={***x15} x51=(x28+36)*x65/x28 ; y-dimension after labeling+shrinking ; label is 36 pixels tall ; shrink back down ip _3 ; INPUT: expanded, labeled class-average [class_avg] x65,x51 ; delete variance if not needed if(x76.eq.1) then cp _7 [class_var]x16 endif ; calculate CCC x33 = x34/x69 ; x34==cumulative CCC-sum, x69==#particles ; get standard deviation of the variance (x44) fs m,x41,x42,x43,x44 _7 ; class-variance _6 ; reduced mask ; trap for variance of a single image if(x15.eq.1)x44=999 ; x15==class-size ; otherwise, the variance of a single image would be NaN ; write CCC (x33) and variance-SD (x44) to doc file sd x16,x16,x33,x44 [class_stats_doc]unsort ; clean up intermediate doc files de [class_doc]{***x16}_noccc de [class_doc]{***x16}_unsort lb2 ; end class-loop vm echo "Reference view: #{***x80}, maximum class size: {***x68}" ; echo ; close document sd e [class_stats_doc]unsort ; sort by CCC doc sort [class_stats_doc]unsort [class_stats_doc] (2) ; column for CCC Y ; renumber? SD / CLASSNUM CCC VARIANCE_SD [class_stats_doc] sd e [class_stats_doc] de [class_stats_doc]unsort ; remove class-map vm rm -f [class_map].$DATEXT ; class-map doc contains each particle with its assigned class vm m rm -f [cas_prefix]_SEQ.$DATEXT ; rm -f [cas_prefix]_SET.$DATEXT ; rm -f [cas_prefix]_PIX.$DATEXT ; rm -f [cas_prefix]_MAS.$DATEXT ; . lb19 ; skip to here if view unpopulated sd x48,x80,x27,x68 [summary_doc] ; x80==view#, x27==#particles, x68=max.class-size lb1 ; end reference-view loop SD / TOTAL_PARTS [summary_doc] [dummy] = -x83 sd [dummy], [total-particles] [summary_doc] ; close docs ud ice [ref_view_list] sd e [summary_doc] vm echo "Done -- classified {******[total-particles]} particles"; date; echo en d ; Modified 2010-07-12 ; TO DO: write unsorted class doc in core ; TO DO: add #classes to summary doc ; 2009-07-13 (trs) -- added summary-file output ; 2009-07-03 (trs) -- prints maximum class size to screen ; 2009-06-22 (trs) -- extra question added by CA SRE when using (I)PCA ; 2009-06-09 (trs) -- can operate on non-consecutive views from list ; 2009-05-27 (trs) -- keeping _EIG file from CA S ; 2009-05-26 (trs) -- option to save eigenimages and reconstituted images, adapted from ca-pca.spi ; 2009-05-26 (trs) -- number of eigenvalues now user-defined ; 2009-05-22 (trs) -- uses selection doc in CA S instead of first N particles ; 2008-11-12 (trs) -- now a parameter to specify CA, PCA, or IPCA ; 2007-11-27 (trs) -- calculates averaged CCC instead of CCC of class-average ; 2007-11-27 (trs) -- can force poorly-populated views into two classes ; 2007-10-11 (trs) -- last reference-view now an input parameter ; 2007-10-01 (trs) -- defocus-group list now optional ; 2007-09-06 (trs) -- input particles are now in stacks ; 2007-09-06 (trs) -- uses iterative PCA instead of CorAn ; 2006-08-29 (trs) -- gets last defocus-group number automatically ; 2006-07-27 (trs) -- bug fix in reference-projection file-pattern ; 2006-04-05 (trs) -- uses last defocus group projections for CCC ; 2006-02-06 (trs & pp) -- updated for changes in projection-matching ; 2005-03-27 (trs) -- no longer needs how_many file ; 2005-01-27 (trs & gra) -- deals with variances of single-image classes ; 2005-01-24 (trs & js) -- bug fix in loop that checks for low intensities ; 2004-12-22 (trs) -- checks for images with low intensities ; 2004-12-22 (trs & jl) -- prints standard deviation of the class-variance to stats document ; 2004-12-22 (trs) -- first reference-view is a parameter, in case you need to resume ; 2004-05-11 (trs) -- handles poorly-populated classes differently ; 2004-05-04 (trs) -- adds option to save/delete variance ; 2004-04-22 (trs) -- sorts individual images by CCC (worst to best) ; 2004-04-16 (trs) -- uses CorAn instead of PCA ; 2004-04-14 (trs) -- deletes unused 'CA S' output (all except _IMC)