; Copyright (c) 1998-2014 A.P. Hitchcock  All rights reserved
;+
;NAME:
;		WRITE_MRC_FILE.PRO
;
;LAST CHANGED: ----------------------------------- 22-Jul-14  [aph,30-Jlu-09,from 06-Jul-06 (gaj)
;
; PURPOSE:
;	This procedure writes  a whole stack (3d (E,x,y) array) to a *.MRC binary file
; 	It assumes Common volume_data, image_stack exists and contains the data
;
;	For definition of the MRC format see:
;	http://bio3d.colorado.edu/imod/doc/mrc_format.txt
;   http://www.msg.ucsf.edu/IVE/IVE4_HTML/IM_ref2.html#Introduction
;
;
;
; CATEGORY:
;	stack processing
;
; CALLING SEQUENCE:
;	WRITE_MRC_FILE, FILE
;
; INPUTS:
; 	FILE	name of the output file
;
; KEYWORDS: none
;
; OUTPUTS:
;	*.mrc file
;
; COMMON BLOCKS:
; AXIS_COM	standard set of common blocks
; STACK_PROCESS_COM
; BSIF_COM
; volume_data, image_stack
;
; MODIFICATION HISTORY:
; (12-apr-06 gaj) first version based on WRITE_WB
; (01-may-06 aph) clean-up
; (06-jul-06 gaj) added extended header to store tilt angles. Fixed bug with scaling of min, max, mean
; (30-jul-09 aph) changed analcom to stack_process_com (update to stack_process)
; (22-Jul-14 aph) corrected header
;-

PRO write_MRC_file, file
@axis_com		; include COMMON from file
@stack_process_com
COMMON volume_data, image_stack
@bsif_com

on_error,2

if strlen(file) EQ 0 then begin
	file = pickfile2(/WRITE, FILTER='*.mrc', /LPATH, DEFPATH=defpath)
endif

if strlen(file) GT 0 then begin

; ****************  For Macintosh computers ***********
;	Note that Intel processors and VAX machines are little endian,
;; while PowerPC processors (e.g., Macs, plus SGI/MIPS, Sun/SPARC),
;; are big endian.  For a number 0x01020304, a big endian machine
;; stores it starting from the "big" end as successive bytes 0x01,
;; 0x02, 0x03, and 0x04.  A little endian machine stores it starting
;; from the "little" end as successive bytes 0x04, 0x03, 0x02, and
;; 0x01.  In any case, we will assume that the binary data was written
;; from an Intel machine in little endian order, so we need to swap
;; the bytes around if we are running this program on a big endian
;; machine.

if !version.OS_family NE 'Windows' AND !version.OS_family NE 'UNIX' then $
	byteorder,image_stack,/swap_if_big_endian

	;image_stack = transpose(image_stack, [1, 0, 2])

	t = ax_name(file)
	fileshort = t(1) + '.' + t(2)

	openw, iunit, file, /get_lun

; --------------- Define and initialize the MRC file header ---------------

	header = {	NX:0L, NY:0L, NZ:0L, $					;# of col, row, sect
				MODE:1L, $								; type of data 1 <-> 16bit signed
				NXSTART:0L, NYSTART:0L, NZSTART:0L, $
				MX:1L, MY:1L, MZ:1L, $
				XLEN:1.0, YLEN:1.0, ZLEN:1.0, $			; cell dimension
				ALPHA:90.0, BETA:90.0, GAMMA:90.0, $	; cell angle
				MAPC:1L, MAPR:2L, MAPS:3L, $
				AMIN:0.0, AMAX:1.0, AMEAN: 0.5, $		; stack statistic
				ISPG:0, NSYMBT:0, NEXT:0L, CREATID:0, $
				EXTRA:bytarr(30), $
				NINT:0, NREAL:0, $						; external header
				EXTRA2:bytarr(28), $
				IDTYPE:0, LENS:0, ND1:0, ND2:0,VD1:1,VD2:2, $
				TILTANGLES:bytarr(24), $
				XORIGIN:0.0, YORIGIN:0.0, ZORIGIN:0.0, $
				CMAP:'MAP ', $
				STAMP:'    ', $
				RMS:0.0, $
				NLABL:0L, $
				LABEL1:bytarr(80), $
				LABEL2:bytarr(80), $
				LABEL3:bytarr(80), $
				LABEL4:bytarr(80), $
				LABEL5:bytarr(80), $
				LABEL6:bytarr(80), $
				LABEL7:bytarr(80), $
				LABEL8:bytarr(80), $
				LABEL9:bytarr(80), $
				LABEL10:bytarr(80) $
				}

; --------------- Setting the parameters in the header ---------------
	im_size = size(image_stack)
	if im_size(0) EQ 2 then n_img = 1 else n_img = im_size(3)

	header.NX = im_size(1)
	header.NY = im_size(2)
	header.NZ = im_size(3)

	header.MODE = 1 	; type of data is 16bit signed

	header.MX = im_size(1)
	header.MY = im_size(2)
	;header.MZ = im_size(3)

	nn = (x_stop-x_start)/im_size(1)
	nn = get_num(val=nn,Prompt='Pixel size [um] (' + strtrim(string(nn),2) + ')', group = Axis_ID)
	print, 'pixel size', nn*1.0e4
	header.XLEN = nn*1.0e4*im_size(1)
	header.YLEN = nn*1.0e4*im_size(2)


; --------------- INTEGER format to save space - MUST re-set scale !! ------------
	tmax = max(image_stack, min = tmin)
	test = max([abs(tmin), abs(tmax)])
	scale = 1.
	if test GT 3e4 OR test LT 1e3 then  scale = 10.^(3-fix(alog10(test)))

	header.AMIN = tmin * scale
	header.AMAX = tmax * scale
	header.AMEAN = mean(image_stack) * scale


; --------------- Creating the MRC extended header ---------------

	; Setting up extended header for tilt angular data
	header.NINT = 2				; 2 bytes for storing tilt angles.
	header.NREAL = 1			; 1 = tilt angle

	ext_header = indgen(n_img)	; Angles stored as (tilt angle * 100) in 2 bytes signed integer
	ta_offset = -90				; tilt angle offset
	ta_offset = get_num(val=ta_offset,Prompt='Tilt angle offset [Deg] ', group = Axis_ID)

	for i = 0, n_img-1 DO ext_header(i) = fix( (ev(i) + ta_offset) * 100)
	axis_log, 'Tilt angles: ' + strtrim(string(ext_header(0)/100),2) + ' to ' + strtrim(string(ext_header(n_img-1)/100),2)

	header.NEXT = header.NINT * n_img ; Calculating and storing the size of extended header



; --------------- Storing information in labels ---------------
	header.LABEL1 = byte('aXis2000: Created MRC stack' + '     ' + systime())
	header.NLABL = header.NLABL + 1
	header.LABEL2 = byte('aXis2000: Scale = ' + strtrim(string(scale),2))
	header.NLABL = header.NLABL + 1

; --------------- Writing data to file
	writeu, iunit, header					; writing header
	writeu, iunit, ext_header				; writing extended header. (Note: header.NEXT must be right)
	writeu, iunit, fix(scale*image_stack)	; writing image stack

	close, iunit

	axis_log, ' Binary MRC file wrote to  ' + fileshort
	if scale NE 1 then 	axis_log, ' scaled by  ' + strtrim(string(scale),2)

ENDIF
END
