# Imports x,y,z positions in a  text file as points in a seperate meshes with material settings controlled by image value
# Meshes each contain multiple data points of similar values
# Designed to work on text files of x,y,z,v format produced by FitsToText
# Produces cubes within each mesh directly as this should be faster and use much less memory

import Blender as B

import Blender
from Blender import *
import math
from math import *
import bpy
from sys import stdout 
import numpy
from numpy import *
import pyfits
from pyfits import *

srange = 50.0							# For colour scaling, maximum pixel value above which everything is white
n = 2000.0								# No. meshes to create
S = 3.0									# For colour scaling, minimum pixel value below which nothing is displayed
count = 0.0

# Count the number of cubes in each mesh created
cc = 0

deltaS = (srange-S)/n
Sf = S + deltaS

scn = Blender.Scene.GetCurrent()

me=B.Mesh.New('myMesh')				# Create initial mesh

ma = Blender.Material.New()		# Create initial material
ma.setMode('Halo')
a = float(S) / srange
ma.setAlpha(a)
ma.setAdd(a)
ma.rgbCol =1.0, 1.0, 1.0
ma.setHaloSize(1.0)

print 'Converting fits file to internal array...'
	
FitsFile = pyfits.open('AGESM33_SN.fits')
image = FitsFile[0].data

z = image.shape[0]	# No. channels					
y = image.shape[1]	# No. x pixels
x = image.shape[2]	# No. y pixels

temp = array(transpose(nonzero(image)),float)	# Make an array of x,y,z of all non-zero values in the fits file
vals = array(image[nonzero(image)])					# Make a list of all those non-zero values

#print vals
l = len(vals)
print l, 'non-zero values found'

temp.resize(l,4)							# Enlarge the array to hold the pixel values

for i in range(0,l):						# Add in the pixel values in the new column
	temp[i,3] = vals[i]
		

# Remove NaNs :
print 'Removing NaNs...'
temp=temp[~numpy.isnan(temp).any(1)]

l = len(temp)

print l,'good values remaining'

print 'Sorting array by value'
temp = temp[temp[:,3].argsort(),:]	# Sorts the array by values in the 4th column. Got it from the net. Don't know how it works.

maxv = temp[l-1,3]

print 'Maximum value :', maxv

#print temp[0:10,:]

templ = temp.shape[0]
tempm = temp[templ-1,3]									# Maximum value since array has already been sorted

print 'Displaying data'


# Read in array of x,y,z,v format
for i in range(0,templ):
	V = temp[i,3]											# x,y,z axes are more or less arbitrary, depends how you want Blender to display
	x = temp[i,2]											# stuff
	y = temp[i,1]
	z = temp[i,0]
	
	#print (float(i)/float(nlines))*100.0			

	Vc = float(V)											# x,y,z axes are more or less arbitrary, depends how you want Blender to display
	xc = float(x)											# stuff
	yc = float(z)
	zc = float(y)
	

	if (Vc >= S) and (Vc < Sf):				# If S/N correct, create a cube at these coordinates with this material
	
		me.verts.extend(xc+0.5,yc-0.5,zc-0.5)
		me.verts.extend(xc-0.5,yc-0.5,zc-0.5)
		me.verts.extend(xc-0.5,yc+0.5,zc-0.5)
		me.verts.extend(xc+0.5,yc+0.5,zc-0.5)

		me.verts.extend(xc-0.5,yc-0.5,zc+0.5)
		me.verts.extend(xc+0.5,yc-0.5,zc+0.5)
		me.verts.extend(xc+0.5,yc+0.5,zc+0.5)
		me.verts.extend(xc-0.5,yc+0.5,zc+0.5)

		facelist = [[0+cc,1+cc,2+cc,3+cc],[0+cc,1+cc,4+cc,5+cc],[2+cc,3+cc,6+cc,7+cc],[1+cc,2+cc,7+cc,4+cc],[3+cc,0+cc,5+cc,6+cc],[5+cc,6+cc,7+cc,4+cc]]
		me.faces.extend(facelist)
		cc = cc + 8
															
	if Vc >= Sf:									# If S/N too high, join the cubes and set up the next material
		count = count + 1.0
		cc = 0
		print (count/n)*100.0
		
		me.materials = [ma]
		for f in me.faces:
			f.mat = 0

		# Create the object and remove duplicate verticies to save memory
		ob=scn.objects.new(me,'Gas')
		me = ob.getData(mesh=1)
		me.remDoubles(0.0)
		ob.drawMode |= Blender.Object.DrawModes.TRANSP	# Enable transparency in realtime view
		
		ma = Blender.Material.New()					# Setup a new material and mesh for the next set of data
		ma.setMode('ZTransp')

		
		if float(Vc) < srange :							# Checks for S/N cutoff
			#a =(float(count)+1)/n
			a = Vc/maxv
		else:
			a = 1.0
			
		ma.setAlpha(a)
		ma.setAdd(a)
		ma.rgbCol =1.0, 1.0, 1.0
		ma.ref = 1.0
		ma.spec = 0.0
		ma.mode |= Material.Modes.SHADELESS
		
		me=B.Mesh.New('myMesh')
		
		S = S + deltaS
		Sf = Sf + deltaS
	
		if S >= srange :					# If above S/N cutoff, set the check to a huge value so only create 1 more object
												# and material			
		#	sys.exit()
			S = srange						
			Sf = maxv-1.0		
				

B.Redraw()

print 'Done !'

# Have empty and camera position moved based on average x,y,z values
# Add axes !!!