# This script uses a 3d delaunay triangulation to draw the surface of
# a volume defined by an atomselection.  Here is an example where we
# draw the volume enclosed by the C-alpha carbons, removing any cells
# with edges > 10.
# 
# source 3d_delaunay_surface.tcl
# set sel [atomselect top "protein and name CA"]
# clear_surface
# graphics $3dgraphicsmol color green
# draw_surface $sel 10
#
# Edit this variable to point to your copy of Hull! 
set hullpath /home/pgrayson/bin/hull

set found 0
foreach m [molinfo list] {
    if { [regexp "3d surface" [molinfo $m get name]] } {
	set found 1
    }
}

if { ! $found } {
    set top [molinfo top]

    mol load graphics "3d surface"

    material add surface
    material change ambient surface 0.000000
    material change specular surface 1.000000
    material change diffuse surface 0.580000
    material change shininess surface 0.540000
    material change opacity surface 0.620000

    molinfo $top set top 1
    display resetview
}

foreach m [molinfo list] {
    if { [regexp "surface" [molinfo $m get name]] } {
        set 3dgraphicsmol $m
    }
}


proc clear_surface { } {
    global 3dgraphicsmol

    graphics $3dgraphicsmol delete all
}

# Draw a tetrahedron
proc draw_tetra { r0 r1 r2 r3 } {
    draw_tri $r0 $r1 $r2
    draw_tri $r0 $r2 $r3
    draw_tri $r0 $r3 $r1
    draw_tri $r3 $r2 $r1
}

# Draw a triangle
proc draw_tri { r0 r1 r2 n0 n1 n2 } {
    global 3dgraphicsmol
    graphics $3dgraphicsmol trinorm $r0 $r1 $r2 $n0 $n1 $n2
}

# Find the max edge of a tetrahedron
proc max { args } {
    set maxval [lindex $args 0]
    foreach arg $args {
        if { $arg > $maxval } {
            set maxval $arg
        }
    }
    return $maxval
}

proc edge_size2 { r0 r1 } {
    return [veclength2 [vecsub $r0 $r1]]
}

proc tri_size2 { r0 r1 r2 } {
    return [max [edge_size2 $r0 $r1] \
                [edge_size2 $r1 $r2] \
                [edge_size2 $r2 $r0]]
}

proc tetra_size2 { r0 r1 r2 r3 } {
    return [max [edge_size2 $r0 $r1] \
                [edge_size2 $r1 $r2] \
                [edge_size2 $r2 $r3] \
                [edge_size2 $r3 $r0] \
                [edge_size2 $r0 $r2] \
                [edge_size2 $r1 $r3]]
}

# keep track of external faces
# the key is the sorted list of vertices, but we preserve the ordering
# as the array member value.
proc xor_external { p0 p1 p2 } {
    global external_faces
    set val "$p0 $p1 $p2"
    set key [lsort $val]
    if { [array get external_faces "$key"] == "" } {
        set external_faces($key) $val
    } else {
        array unset external_faces "$key"
    }
}

# this one calls xor_external, but only after correctly setting the
# direction of the triangle, assuming it is part of the convex hull
proc xor_external_outer { p0 p1 p2 } {
    global coordarray
    foreach testp {0 1 2 3} {
        if { $p0 != $testp && $p1 != $testp && $p2 != $testp } {
            set p3 $testp
        }
    }

    if { [tetra_ordered $coordarray($p0) $coordarray($p1) \
            $coordarray($p2) $coordarray($p3)] } {
        set temp $p0
        set p0 $p1
        set p1 $temp
    }
    xor_external $p0 $p1 $p2
}

# draw the external faces
proc draw_external { } {
    global point_normals
    global coordarray
    global external_faces
    foreach facekey [array names external_faces] {
        set face $external_faces($facekey)
        set p0 [lindex $face 0]
        set p1 [lindex $face 1]
        set p2 [lindex $face 2]
        draw_tri $coordarray($p0) $coordarray($p1) $coordarray($p2) \
                $point_normals($p0) $point_normals($p1) $point_normals($p2) 
    }
}

# shade the external faces
proc shade_external { } {
    global point_normals
    global coordarray
    global external_faces
    array set point_normals ""
    array unset point_normals "*"

    # clear the array
    foreach point [array names coordarray] {
        set point_normals($point) "0 0 0"
    }
    # fill the array
    foreach facekey [array names external_faces] {
        set face $external_faces($facekey)
        set p0 [lindex $face 0]
        set p1 [lindex $face 1]
        set p2 [lindex $face 2]
        set vec1 [vecsub $coordarray($p0) $coordarray($p1)]
        set vec2 [vecsub $coordarray($p0) $coordarray($p2)]
        set normal [veccross $vec1 $vec2]
        
        set point_normals($p0) [vecadd $point_normals($p0) $normal]
        set point_normals($p1) [vecadd $point_normals($p1) $normal]
        set point_normals($p2) [vecadd $point_normals($p2) $normal]
    }
    # normalize the array
    foreach point [array names coordarray] {
        if { $point_normals($point) != "0 0 0" } {
            set point_normals($point) [vecnorm $point_normals($point)]
        }
    }
}

# check the order of the tetra
proc tetra_ordered { r0 r1 r2 r3 } {
    set vec1 [vecsub $r2 $r0]
    set vec2 [vecsub $r3 $r0]
    set vec3 [vecsub $r1 $r0]
    set normal [veccross $vec1 $vec2]
    set dot [vecdot $normal $vec3]
    return [expr $dot >= 0]
}

# Give me a bunch of points, and I will plot their surface
proc draw_surface { sel cutoff } {
    global coordarray
    global external_faces
    global tcl_platform
    global 3dgraphicsmol
    global hullpath

    graphics $3dgraphicsmol material surface

#    [ atomselect top "index [$lipidsel list] and name P N" ]
    set coords [ $sel get {x y z} ]
    set headsx [ $sel get x ]
    set headsy [ $sel get y ]
    set headsz [ $sel get z ]
    set num [llength [$sel list]]

    # write the points to a file, and set up the coordinate array
    array set coordarray ""
    array unset coordarray "*"
    set i 0
    foreach x $headsx y $headsy z $headsz {
        append points "$x $y $z\n"
        set coordarray($i) "$x $y $z"
        incr i
    }

    # compute the delaunay tetras
    set delaunay [exec $hullpath -m 20 -X/dev/null -d << $points]

    set delaunay [split $delaunay "\n"]
    
    # compute the squared cutoff
    set cutoff2 [expr $cutoff*$cutoff]

    # read out the triangles
    array set external_faces ""
    array unset external_faces "*"
    foreach line $delaunay {
        if { ! [regexp "^\%" $line] } {
            if { [llength $line] == 4 } {
                set p0 [lindex $line 0]
                set p1 [lindex $line 1]
                set p2 [lindex $line 2]
                set p3 [lindex $line 3]

                if { $p0 == -1 } {
                    xor_external_outer $p1 $p2 $p3
                } elseif { $p1 == -1 } {
                    xor_external_outer $p0 $p3 $p2
                } elseif { $p2 == -1 } {
                    xor_external_outer $p0 $p1 $p3
                } elseif { $p3 == -1 } {
                    xor_external_outer $p0 $p2 $p1
                } elseif { [tetra_size2 $coordarray($p0) $coordarray($p1) \
                        $coordarray($p2) $coordarray($p3)] > $cutoff2 } {
                    if { ! [tetra_ordered $coordarray($p0) $coordarray($p1) \
                            $coordarray($p2) $coordarray($p3)] } {
                        set temp $p0
                        set p0 $p1
                        set p1 $temp
                    }
                    xor_external $p0 $p1 $p2
                    xor_external $p0 $p2 $p3
                    xor_external $p0 $p3 $p1
                    xor_external $p1 $p3 $p2
                }
            } else { puts "Bad line: $line" }
        }
    }
    shade_external
    draw_external
}