# Inspecting atlas volume data using SimpleITK and NumPy

import SimpleITK as sitk
import numpy as np

# First download the BMA 2017 zip file and unzip

# Load some of the data 
img_mri = sitk.ReadImage("./bma-1-all/bma-1-mri.nii.gz")
img_label = sitk.ReadImage("./bma-1-all/bma-1-region_seg.nii.gz")

# Print some information about the volumes:
spacing = img_mri.GetSpacing() # Will use this later.
print("Spacing (voxel size):")
print(spacing)
print("Origin:")
print(img_mri.GetOrigin())
print("Get the directions describing the volume's coordinate system:")
print(img_mri.GetDirection())
size = img_mri.GetSize()
print("Size of the atlas volume:")
print(size)

# Get a numpy array from the label map image
nda_label = sitk.GetArrayFromImage(img_label)

# Get a list of the IDs
label_ids = np.unique(nda_label)
print("Number of unique labels (0 is the background):")
print(len(np.unique(label_ids)))

# Get the volume in voxels for V1 of the right hemisphere (region with ID = 142)
v1_voxel_count = np.count_nonzero(nda_label == 142)
print("V1 (right hemisphere) voxel count:"+str(v1_voxel_count))

# Convert to mm^32
v1_size = v1_voxel_count*(spacing[0]*spacing[1]*spacing[2])
print("V1 (right hemisphere) size mm^3: "+str(round(v1_size,2)))

# Get a slice from the volume
A = 500
slice = sitk.GetArrayFromImage(img_mri)[:,A,:]
# Get the array data type
print(slice.dtype)

# See more commands at http://simpleitk.org/SimpleITK-Notebooks/01_Image_Basics.html