% -------------------------------
%
% Download a NRRD reader
% For example:
% http://www.mathworks.com/matlabcentral/fileexchange/34653-nrrd-format-file-reader
%
% Requires: MATLAB 7.13 (R2011b)
%
% Download average_template_25.nrrd,
%          ara_nissl_25.nrrd,
%          ccf_2015/annotation_25.nrrd
%
% ---------------------------------

%
% Read image volume with NRRD reader
% Note that reader swaps the order of the first two axes
%
% AVGT = 3-D matrix of average_template
% NISSL = 3-D matrix of ara_nissl
% ANO = 3-D matrix of ccf_2015/anntotationannotation
%
[AVGT, metaAVGT] = nrrdread('average_template_25.nrrd');
[NISSL, metaNISSL] = nrrdread('ara_nissl_25.nrrd');
[ANO, metaANO] = nrrdread('annotation_25.nrrd');


% Display one coronal section
figure;imagesc(squeeze(AVGT(:,264,:)));colormap(gray(256)); axis equal;
figure;imagesc(squeeze(NISSL(:,264,:)));colormap(gray(256)); axis equal;
figure;imagesc(squeeze(ANO(:,264,:)));
caxis([1,2000]); colormap(lines(256)); axis equal;

% Display one sagittal section
figure;imagesc(squeeze(AVGT(:,:,220)));colormap(gray(256)); axis equal;
figure;imagesc(squeeze(NISSL(:,:,220)));colormap(gray(256)); axis equal;
figure;imagesc(squeeze(ANO(:,:,220)));
caxis([1,2000]); colormap(lines(256)); axis equal;
 

# -------------------------------
#
# Install pynrrd: https://github.com/mhe/pynrrd
#
# Download average_template_25.nrrd,
#          ara_nissl_25.nrrd,
#          ccf_2015/annotation_25.nrrd
#
# ---------------------------------
import nrrd
import numpy as np
import matplotlib.pyplot as plt
import Image

#
# Read image volume with NRRD reader
# Note that reader swaps the order of the first two axes
#
# AVGT = 3-D matrix of average_template
# NISSL = 3-D matrix of ara_nissl
# ANO = 3-D matrix of ccf_2015/anntotationannotation
#

AVGT, metaAVGT = nrrd.read('average_template_25.nrrd');
NISSL, metaNISSL = nrrd.read('ara_nissl_25.nrrd');
ANO, metaANO = nrrd.read('annotation_25.nrrd');

# Save one coronal section as PNG
slice = AVGT[264,:,:].astype(float)
slice /= np.max(slice)
im = Image.fromarray(np.uint8(plt.cm.gray(slice)*255))
im.save('output/avgt_coronal.png')

slice = NISSL[264,:,:].astype(float)
slice /= np.max(slice)
im = Image.fromarray(np.uint8(plt.cm.gray(slice)*255))
im.save('output/nissl_coronal.png')

slice = ANO[264,:,:].astype(float)
slice /= 2000
im = Image.fromarray(np.uint8(plt.cm.jet(slice)*255))
im.save('output/ano_coronal.png')

# Save one sagittal section as PNG
slice = AVGT[:,:,220].astype(float)
slice /= np.max(slice)
im = Image.fromarray(np.uint8(plt.cm.gray(slice)*255))
im.save('output/avgt_sagittal.png')

slice = NISSL[:,:,220].astype(float)
slice /= np.max(slice)
im = Image.fromarray(np.uint8(plt.cm.gray(slice)*255))
im.save('output/nissl_sagittal.png')

slice = ANO[:,:,220].astype(float)
slice /= 2000
im = Image.fromarray(np.uint8(plt.cm.jet(slice)*255))
im.save('output/ano_sagittal.png')