#!/usr/bin/env python
# -*- encoding: utf-8 -*-

# Copyright (c) 2012, tamanegi (tamanegi@users.sourceforge.jp)
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# 
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
# 
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
# IN THE SOFTWARE.

import sys
import copy
import os

filedirname = os.path.dirname( os.path.abspath( __file__ ) )
commonpath = os.path.join( filedirname, "../common" )
if not commonpath in sys.path:
  sys.path.append( commonpath )
from vector3d import *

# Voxel data

# a Voxel
class Voxel:
  # class-static variables
  # something lattice-vector-like object
  lattice = [ [0,0,0], [1,0,0], [0,0,1], [1,0,1], [0,1,0], [1,1,0], [0,1,1], [1,1,1] ]
  # list of tetrahedrons
  #tetrahedrons = [ [0,1,3,4], [0,2,3,4], [2,3,4,6], [1,3,4,5], [3,4,5,7], [3,4,6,7] ]
  tetrahedrons = [ [3,4,0,1], [3,4,0,2], [3,4,2,6], [3,4,1,5], [3,4,5,7], [3,4,6,7] ]
  # list of midpoints
  midpoints = [ [0,1], [0,2], [0,3], [1,2], [1,3], [2,3] ]
  def __init__( self ):
    # first 4 vertexes locates in the same square, last 4 vertexes do so
    #  z   2___3  6___7    y
    #  ^   |   |  |   |   /
    #  |   |___|  |___|  /
    #  |   0   1  4   5  -> x
    # 0-4, 1-5, 2-6, and 3-7 are also connected
    self.mydata = [ 0, 0, 0, 0, 0, 0, 0, 0 ]
    # origin of the voxel; position of grid point 0
    self.origin = Vector3D( 0.0, 0.0, 0.0 )
    # grid size for each direction
    # note that a voxel is always cubic, the bond angles are always 90 degrees.
    self.grid = Vector3D( 1.0, 1.0, 1.0 )

  def __getitem__( self, index ):
    return self.mydata[index]

  def __setitem__( self, index, value ):
    self.mydata[index] = value

  # pos must be Vector3D type
  def setOrigin( self, pos ):
    self.origin = Vector3D( pos.x, pos.y, pos.z )

  def setGridLength( self, vec ):
    self.grid = Vector3D( vec.x, vec.y, vec.z )

  def setValueFromList( self, values ):
    # check length of the list
    if len( values ) != 8:
      print >> sys.stderr, "Error(Voxel): voxel data length must be 8."
      sys.exit(1)
    self.mydata = copy.copy( values )

  # generate XML data by marching tetrahedron method
  def generateXml( self, isosurface, toWhat = sys.stdout, positive = True ):
    self.__generateXml( isosurface, toWhat, positive )

  def __generateXml( self, isosurface, toWhat, positive ):
    for hedron in Voxel.tetrahedrons:
      self.__gen( toWhat, isosurface, hedron, positive )

  def __gen( self, toWhat, isosurface, hedron, positive ):
    dat = [ 0, 0, 0, 0 ]
    num = 0 # number of found points
    counter = 0
    for vertex in hedron:
      if positive:
        if self.mydata[vertex] > isosurface:
          dat[counter] = 1
          num += 1
      else:
        if self.mydata[vertex] < -isosurface:
          dat[counter] = 1
          num += 1
      counter += 1

    if num == 0 or num == 4:
      return

    pos = []
    myval = isosurface
    if not positive:
      myval = - isosurface
    #mplist = []
    for mp in Voxel.midpoints:
      if dat[mp[0]] != dat[mp[1]]:
        #mplist.append( mp )
        #print Voxel.lattice[hedron[mp[0]]], Voxel.lattice[hedron[mp[1]]], self.mydata[hedron[mp[0]]], self.mydata[hedron[mp[1]]]
        pos.append( self.__getMidpointPos( myval, Voxel.lattice[hedron[mp[0]]], Voxel.lattice[hedron[mp[1]]], self.mydata[hedron[mp[0]]], self.mydata[hedron[mp[1]]] ) )

    ## debug message
    #print hedron, dat
    #for p in pos:
    #  print "      <ATOM pos=\"" + str(p) + "\" color=\"red\" />"

    if len(pos) == 3:
       # draw 0-1-2 trangles
      if self.__switchVector( dat[0], Voxel.lattice[hedron[0]], pos[0], pos[1], pos[2] ):
        self.__printTri( pos[0], pos[2], pos[1] )
      else:
        self.__printTri( pos[0], pos[1], pos[2] )
    elif len(pos) == 4:
      if self.__switchVector( dat[0], Voxel.lattice[hedron[0]], pos[0], pos[1], pos[3] ):
        self.__printTri( pos[0], pos[3], pos[1] )
      else:
        self.__printTri( pos[0], pos[1], pos[3] )
      if self.__switchVector( dat[0], Voxel.lattice[hedron[0]], pos[0], pos[2], pos[3] ):
        self.__printTri( pos[0], pos[3], pos[2] )
      else:
        self.__printTri( pos[0], pos[2], pos[3] )
    else:
      print >> sys.stderr, "Error: unexpected error in creating polygons."
      sys.exit(1)

  def __printTri( self, p0, p1, p2, toWhat = sys.stdout ):
    print >> toWhat, "      <TRI p0=\"" + str(p0) + "\" p1=\"" + str(p1) + "\" p2=\"" + str(p2) + "\" />"

  def __switchVector( self, dat0, gp, p0, p1, p2 ):
    mypos = Vector3D( self.grid.x * float( gp[0] ), self.grid.y * float( gp[1] ), self.grid.z * float( gp[2] ) ) + self.origin
    vec = mypos - p0
    cvec = ( p1 - p0 ).cross( p2 - p0 )
    val = vec.dot( cvec )
    if val == 0.0:
      print >> sys.stderr, "Warning: unexpected value."
    if dat0 == 0:
      # vertex[0] outside; face normal should stick toward vertex[0]
      if val > 0.0:
        return False
      else:
        return True
    else:
      # vertex[0] inside
      if val > 0.0:
        return True
      else:
        return False

  def __diagonal( self, m0, m1 ):
    if m0[0] == m1[0] or m0[0] == m1[1]:
      return False
    if m0[1] == m1[0] or m0[1] == m1[1]:
      return False
    return True

  def __getMidpointPos( self, myval, vertex0, vertex1, val0, val1 ):
    v = ( myval - val1 ) / ( val0 - val1 )
    # avoid vertex point itself
    #w = 1.0 - v
    w = min( max( 1.0 - v, 0.01 ), 0.99 )
    v = 1.0 - w
    x = self.grid.x * ( v * float( vertex0[0] ) + w * float( vertex1[0] ) )
    y = self.grid.y * ( v * float( vertex0[1] ) + w * float( vertex1[1] ) )
    z = self.grid.z * ( v * float( vertex0[2] ) + w * float( vertex1[2] ) )
    return self.origin + Vector3D( x, y, z )
