; Copyright (c) 1998-2014 A.P. Hitchcock  All rights reserved
;+
;NAME:
;		READ_MRC_FILE.PRO
;
;LAST CHANGED: ----------------------------------- 22-Jul-14 (aph)
;
; PURPOSE:
;	This procedure reads a whole stack (3d (x,y,z) array) from a *.MRC binary file
;   details of the mrc format can be found at:
;	http://bio3d.colorado.edu/imod/doc/mrc_format.txt
;	http://ami.scripps.edu/prtl_data/mrc_specification.htm
;
; CATEGORY:
;	stack processing
;
; CALLING SEQUENCE:
;	READ_MRC_FILE, FILE
;
; INPUTS:
; 	FILE	name of the input file
;
; KEYWORDS: none
;
; OUTPUTS: none
;
; COMMON BLOCKS:
; AXIS_COM	standard set of common blocks
; STACK_PROCESS_COM
; BSIF_COM
; volume_data, image_stack
;
; SIDE EFFECTS:
;
;
; RESTRICTIONS: none
;
; MODIFICATION HISTORY:
; (12-Jun-06 gaj) first version based on WRITE_MRC_FILE and STACK_RB
; (13-Jun-06 gaj) corected some variable names and bugs
; (30-jul-09 aph) changed analcom to stack_process_com (update to stack_process)
; (22-Jul-14 aph) corrected header
;-

PRO read_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(/READ, FILTER='*.mrc', /LPATH, DEFPATH=defpath)
endif

if strlen(file) GT 0 then begin


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

	; Define 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, $
				ISPG:0, NSYMBT:0, NEXT:0L, CREATID:0, $
				EXTRA:bytarr(30), $
				NINT:0, NREAL:0, $
				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) $
				}


	openu, iunit, file, /get_lun

; Reads the MRC header
	readu, iunit, header
;	print, header


; Get extended header
	if header.NEXT NE 0 then begin
		ext_header = bytarr(header.NEXT)
		readu, iunit, ext_header
	endif


;Want to do this, but don't know how:
;	for i=0 to header.NLABL
;		axis_log, header.LABEL?
;	endfor
;Have to do this insted:
	axis_log, 'MRC labels:'
	axis_log, string(header.LABEL1)
	axis_log, string(header.LABEL2)
	axis_log, string(header.LABEL3)
	axis_log, string(header.LABEL4)
	axis_log, string(header.LABEL5)
	axis_log, string(header.LABEL6)
	axis_log, string(header.LABEL7)
	axis_log, string(header.LABEL8)
	axis_log, string(header.LABEL9)
	axis_log, string(header.LABEL10)


	; define stack array size and read the stack
	if header.MODE EQ 1 then begin
		image_stack = intarr(header.NX, header.NY, header.NZ)
    	readu, iunit, image_stack
    	axis_log, 'Read Binary MRC file from ' + fileshort
    endif


	; setup up aXis stack parameters from header

	n_cols = header.NX
	n_rows = header.NY
	n_sect = header.NZ
	text = string(format='(" Image size: ", i3," x ",i3," pixels")',n_cols,n_rows)
	axis_log, text

	x_start = header.NXSTART
	x_stop = header.XLEN / 10000
	y_start = header.NYSTART
	y_stop = header.YLEN /10000
	text = string(format='("X-range (um):  ",f7.3,"  to  ",f7.3)',x_start, x_stop)
	axis_log, text
	text = string(format='("Y-range (um):  ",f7.3,"  to  ",f7.3)',y_start, y_stop)
	axis_log, text


	n_energies = header.NZ
	print, ' # of images = ', fix(n_energies)

	ev = fltarr(n_energies)
	dwell = fltarr(n_energies)
	filename_ev_msec_list = strarr(n_energies)

	; Setting up the Z dimension to go in the energy scale
	for i = 0, n_energies-1 do begin
		ev(i) =  header.ZLEN/header.NZ*i/10000 ; '/10000 is to convert from Angstrom to um
		filename_ev_msec_list(i) = string(ev(i)) + ' um'
	endfor

	text = string(format='("Z-range (um):  ",f7.2,"  to  ",f7.2)', min(ev), max(ev))
	axis_log, text

   	close, iunit & free_lun, iunit


ENDIF
END
