#!/usr/bin/env python

# ****************************************************
# **    CMOS and Hybrid pCT analysis routine        **
# **    written by Michela Esposito                 **
# **    for the PRaVDA collaboration                **
# **                                                **
# **    Edit by Laurent Kelleter                    **
# ****************************************************

import numpy as np
import os

#read binary image format. Returns:images,info,no_images (numpy ndarray containing all readout image, header info and number of images)
def BinaryReader(path,filename):
	print('Reading %s from %s'%(filename,path))
	f=open(os.path.join(path,filename),'rb')
	f.seek(0)
	f.seek(35)
	icd_len=np.fromfile(f,dtype=np.uint32,count=1)[0]
	f.seek(39+icd_len+34)
	stream_path_len=np.fromfile(f,dtype=np.uint32,count=1)[0]
	dt_header=np.dtype([
		('version',np.uint16),
		('size',np.int64),
		('image',[
			('left',np.uint32),
			('top',np.uint32),
			('width',np.uint32),
			('height',np.uint32),
			('force',np.bool_),
			('type',np.uint32),
		]),
		('icd',[
			('lead',(np.uint8,4)),
			('size',np.uint32),
			('path',(np.uint8,icd_len)),
		]),
		('stream',[
			('systems',np.int32),
			('boards',np.int32),
			('halves',np.int32),
			('images',np.float64),
			('rate',np.int32),
			('samples',np.uint32),
			('mask',np.uint16),
			('path',[
				('lead',(np.uint8,4)),
				('size',np.uint32),
				('path',(np.uint8,stream_path_len)),
			]),
			('mode',np.uint16),
			('max_time',np.int32),
		])

	])

	f.seek(0)
	header=np.fromfile(f,dtype=dt_header,count=1)

	dt_images=np.dtype(
		(np.uint16,(
			header['stream']['samples'][0],
			header['image']['height'][0],
			int(header['image']['width'][0]/2)
		))
	)
	f.seek(39+icd_len+38+stream_path_len+6)
	images=np.fromfile(f,dtype=dt_images,count=1)[0]
	return images,header['image'],header['stream']['samples']

#discard noisy images at the start of a streaming files
#remove defective rows is optional because it distorts the image
def DiscardImageRows(images):
  images=np.delete(images,[0,1,2],axis=0)
  ##delete a couple of defect rows
  #nY=len(images[0])
  #images=images[:,5:,:]
  #for i in range(5):
  #  images = np.delete(images, (nY/2-5), axis=1)
  return images

#rotate 2nd half (bottom half) of the image to give correct physical orientation
def Rotate2Half(image):
	image_rot=np.zeros(image.shape)
	half1=image[:int(image.shape[0]/2),:]
	half2=image[(int(image.shape[0]/2)):image.shape[0]:,:]
	half2=np.fliplr(np.flipud(half2))
	image_rot[:(int(image.shape[0]/2)),:]=half1
	image_rot[(int(image.shape[0]/2)):image.shape[0],:]=half2
	return image_rot

#rotate second half of list of images
def Rotate2Halfs(images):
  images_rot=[]
  for image in images:
    images_rot.append(Rotate2Half(image))
  return images_rot

#Do the standard treatment to the images: Remove noisy images, sum them up, scale and rotate
def StandardTreatment(images):
  no_images=len(images)
  #rotate bottom half of the image to give correct physical orientation - this is due to both sensor halves being readout row by row in opposite directions
  images_rot=Rotate2Halfs(images)
  #discard first 3 noisy images and noisy rows
  images_rot=DiscardImageRows(images_rot)
  image_rot_sum=np.sum(images_rot,axis=0)
  #Scale to the number of images = no_images minus number of deleted images
  image_rot_sum=1./(no_images-3)*image_rot_sum
  return image_rot_sum