# ExportMultiSectColladaKMZ.sml

# write multiple geologic cross-section images to COLLADA files packaged in a
# KMZ file for viewing in Google Earth. The contained KML file animates extruding the
# sections from their true below-ground positions (where they are not visible in Google
# Earth). The elevation range for the animation is automatically set to raise 
# all sections above the terrain while keeping their correct relative elevations.

# Cross-sections must have 3D Piecewise Affine (manifold) georeference model
# with 3D control points.  Each section's coordinate reference system can be any
# projected CRS or a nonprojected geographic (latitude/longitude) CRS.

# Animation uses a separate placemark for each step for each section, each referencing the 
# cross-section model with different start/end times and altitudes. Works with
# all versions of Google Earth.  

# Randy Smith, MicroImages, Inc.
# 29 July 2010

class RVC_OBJITEM tempObjItemList[];	# hash of ObjItems for the selected cross-section rasters 
class RVC_OBJITEM xsecObjItemList[];	# hash of ObjItems for valid cross-section rasters
class RVC_RASTER XSEC;
class RVC_GEOREFERENCE xsecgeoref, tempgeoref;
numeric i, j, k; # counters
numeric err; # for error checking
numeric numSections;
numeric minZm, maxZm;	# min and max elevations among all section control points (always stored in meters)
minZm = 30000; # starting value higher than any terrestrial elevation in meters;
maxZm = -1000;	# starting value lower than any terrestrial elevation in meters

numeric numPts;

# number of steps in extrusion animation 
# (= number of cross-section placemarks to write to KML file) 
numeric numExtrudeSteps = 10;

clear();

##############   FUNCTIONS AND PROCEDURES   #########################

# error checking procedure
proc ReportError(numeric linenum, numeric err) {
	printf("FAILED -line: %d, error: %d\n", linenum - 1, err);
	PopupError(err); 
	}

# function to determine if three points are in counter-clockwise order
# using image coordinates (origin in upper left, positive downward);
# returns value < 0 if counter-clockwise, value > 0 if clockwise.
func ccwImageCoords(class POINT2D p1, class POINT2d p2, class POINT2d p3)
	{
	return (p2.x - p1.x)*(p3.y - p1.y) - (p2.y - p1.y)*(p3.x - p1.x);
	}




###############################################################################
##############################   MAIN PROGRAM   ###############################

############### Get input cross-sections

class STRING prompt$ = "Choose cross-section rasters to export to COLLADA:";

DlgGetObjects(prompt$, "Raster", tempObjItemList, "ExistingOnly", 1);
numSections = tempObjItemList.GetNumItems();
printf("%d rasters selected.\n", numSections );
array numeric statusArray[numSections]; # array to record whether section rasters are valid

# check for manifold georeference and delete objitems of invalid objects from xsecObjItemList
for i = 1 to numSections
	{
	XSEC.Open(tempObjItemList[i], "Read");

	# get georeference subobject and check that it is present and manifold type with valid CRS
	if (XSEC.GetDefaultGeoref(tempgeoref) )	# true if object has a georeference subobject
		{
		if (tempgeoref.GetGeorefType() == "Manifold" )
			{
			printf("Section object %d has manifold georeference.\n", i);
			if (tempgeoref.GetCoordRefSys().IsLocal() )
				{
				print("Section raster %d has arbitrary local georeference. Removing from input list.\n", i);
				}
			else
				{
				statusArray[i] = 1;
				xsecgeoref.OpenLastUsed(XSEC);
				numPts = xsecgeoref.GetNumCtrlPoints();
				class CTRLPOINT TctrlPtList[numPts];
				xsecgeoref.GetPointList(TctrlPtList);

				for j = 1 to numPts
					{ 
					if (TctrlPtList[j].dest.z < minZm) then minZm = TctrlPtList[j].dest.z;
					if (TctrlPtList[j].dest.z > maxZm) then maxZm = TctrlPtList[j].dest.z;
					}
				xsecgeoref.Close();
				}
			}
		else
			{
			printf("Section object %d does not have manifold georeference. Removing from input list.\n", i);
			statusArray[i] = 0;
			}
		}
	else
		{
		printf("Section object %d is not georeferenced. Removing from input list.\n", i);
		statusArray[i] = 0;
		}

	XSEC.Close();
	}

# loop through statusArray to copy objItems for sections with valid manifold georeference
# to the final objItemList
numeric validcount;		# count of valid sections to be moved to the new objItemList

for i = 1 to numSections
	{
	if (statusArray[i] == 1 ) 
		{
		++validcount;
		xsecObjItemList[validcount] = tempObjItemList[i];
		}
	}

if (xsecObjItemList.GetNumItems() == 0 )
	{
	print("No selected cross-section raster has valid manifold georeference. Exiting now.");
	Exit();
	}

printf("Number of valid cross-sections to process = %d\n", xsecObjItemList.GetNumItems() );
printf("Minimum elevation in meters among all section control points = %.1f\n", minZm);
printf("Maximum elevation in meters among all section control points = %.1f\n", maxZm); 

# prompt user for the name of the output KMZ file
class STRING kmzfilename$;
prompt$ = "Choose and output KMZ file to package the cross-section data \n (spaces will be replaced with underscores):"; 
kmzfilename$ = GetOutputFileName(_context.ScriptDir, prompt$, ".kmz");

for i = 0 to kmzfilename$.length - 1
	{
	if (kmzfilename$.charAt(i) == " ") then
		kmzfilename$.replace(i, i, "_");
	}

class FILEPATH kmzFilepath(kmzfilename$);
class FILEPATH dirFilepath(kmzfilename$);	# filepath to parent directory of KMZ file
dirFilepath.StripToFolder();	# path to directory containing the KMZ
class STRING dir$;
dir$ = dirFilepath;

class STRING dataName$;
dataName$ = kmzFilepath.GetNameOnly();

# set up a WGS84 / Geographic CRS for lat/lon coordinates to position sections
class SR_COORDREFSYS geoCRS;
geoCRS.Assign("Geographic2D_WGS84_Deg");
printf("Geographic CRS = %s\n", geoCRS.Name);

####################################################################
# string with COLLADA file template for a cross-section 
class STRING dae$ = '
   
   
      
         TNTmips SML Script
      
      
      
      
      Z_UP
   
   
      
         
      
   
   
      
         
      
   
   
      
         
            
               
                  Xsect-image
               
            
            
               
                  Xsect-image-surface
               
            
            
               
                  
                     0.000000 0.000000 0.000000 1
                  
                  
                     0.000000 0.000000 0.000000 1
                  
                  
                     
                  
                  
                     0.330000 0.330000 0.330000 1
                  
                  
                     20.000000
                  
                  
                     0.100000
                  
                  
                     1 1 1 1
                  
                  
                     0.000000
                  
               
            
            
               
                  1
               
            
         
      
   
   
      
         
            
               
               		 
               
                  
                     
                     
                     
                  
               
            
            
               -0.000000 -1.000000 -0.000000 
               
                  
                     
                     
                     
                  
               
            
            
               
               		 
               
                  
                     
                     
                  
               
            
            
               
            
            
               
               
               
               
 
            
         
      
   
   
      
         
            
               
                  
                     
                        
                           
                        
                     
                  
               
            
         
      
   
   
      
   

';

# variables for modifying the COLLADA template for each section
class DATETIME dt;
class STRING dt$, format$;
class XMLDOC daedoc;
class XMLNODE colladaRootNode;
class XMLNODE node, node1, node2, node3;
numeric arraycount;

###############################################################
# string with template for KML file
class STRING kmldoc$;
class STRING kmlfile$ = dir$ + "/" + "doc.kml";
class FILEPATH kmlFilepath(kmlfile$);

kmldoc$ = '


	
	1

';

# variables for modifying the KML file template
class XMLDOC kmldoc;
class XMLNODE kmlRootNode;
class XMLNODE folder, document, placemark, model;
numeric startDate;
numeric endDate;
numeric altitudeIncrement;



##############################
# parse KML file template to reference and position the COLLADA models
err = kmldoc.Parse(kmldoc$);
if (err < 0) {	PopupError(err); Exit(); }

kmlRootNode = kmldoc.GetRootNode();	# get root node 
folder = kmlRootNode.FindChild("Folder");		# get folder node that will hold a Document for each cross-section
node = folder.FindChild("name");
node.SetText(dataName$);		# set name of data collection



##################################################################
################ Create KMZ file (ZIP file with .kmz file extension
class ZIPFILE kmz;


###########################################################################
# variables for processing the cross-section rasters and their control points
class RVC_DESCRIPTOR xsecDescriptor;
class STRING xsecName$;
class SR_COORDREFSYS xsecCRS;
class STRING unit$;
numeric geogCRSflag;
class SR_COORDAXIS axisX;
numeric unitConvToMeters;		# scale factor to convert CRS map unit to meters
numeric metersConvToUnit;		# scale factor to convert meters to CRS map unit

class BOUNDARYLIST bndryList;		# container for boundary point arrays
array numeric bndryArray[1];	# array indices for the current boundary
array numeric xyz_array[1];	# holds 3D map coordinates for all points: X Y Z X Y Z etc.
array numeric uv_array[1];	# holds normalized texture coordinates: S T S T etc.

# arrays to hold map coordinates (destination) used to determine coordinate
# ranges and compute local 3D position coordinates for the COLLADA model 
array numeric xarray[1];
array numeric yarray[1];
array numeric zarray[1];

numeric minX, minY, minZ;		# minimum values for each of x, y, and z coordinates for the section control points in CRS units
numeric maxZ;						# maximum z value from control points in CRS units

# variables for processing control point coordinates
class STRING xyzString, uvString;		# string arrays of coordinates to write to COLLADA
numeric col, lin, s, t;
class POINT3D destPt;

class POINT3D originPt;		# origin point for COLLADA model to place it correctly in Google Earth
numeric angleToNorth;
altitudeIncrement = (maxZm - originPt.z) / numExtrudeSteps;

# variables for creating and processing triangulation from the control points
class TRIANGULATOR triangulator; # Triangulator class to create a triangulation in image space for the control points 
class TRIANGLENODES triNodes;
class CTRLEDGE hardedge;
numeric nodeNum, nodeIndex, diffNode;
numeric numTriangles;
class POINT2D pt, nodePt;

array numeric triNodeIndexArray[3];
array numeric tempArray[3];
numeric testCCW;
class STRING pString;	# string array of indices into coordinate arrays to represent their triangulations in COLLADA file



#########################################################################
#########################################################################
#### loop through objItemList of valid sections to process them 

for k = 1 to validcount
	{
	printf("\nProcessing section %d of %d\n", k, validcount);
	XSEC.Open(xsecObjItemList[k], "Read");
	xsecgeoref.OpenLastUsed(XSEC);
	xsecDescriptor = XSEC.GetDescriptor();
	xsecName$ = xsecDescriptor.GetShortName();
	printf("Name of cross-section raster %d is: %s\n", k, xsecName$);

	# get Coordinate Reference System from section
	xsecCRS = xsecgeoref.GetCoordRefSys();
	printf("Section %d CRS = %s\n", k, xsecCRS.Name);

	if (xsecCRS.IsProjected() )	# projected CRS with cartesian map coordinates
		{
		printf("Section %d has projected CRS.\n", k);
		axisX = xsecCRS.CoordSys.GetAxis();
		unit$ = axisX.Unit.GetName();

		# get scale factors for converting between units in the projected CRS and meters.
		# Scale factor from unit to meters is required as a parameter in the COLLADA file.
		# Elevations in manifold georeference are always stored as meters, so need to convert
		# elevation values to units in projected CRS as well.

		unitConvToMeters = GetUnitConvDist(unit$, "meter"); 
		printf("Map unit = %s, conversion factor to meters = %.4f\n", unit$, unitConvToMeters);
		metersConvToUnit = GetUnitConvDist("meter", unit$);
		}
	else 
		{					# nonprojected = geographic (latitude longitude)
		printf("Section %d has geographic CRS.\n", k);
		geogCRSflag = 1;	# set flag indicating CRS is geographic 
		unit$ = "meters";
		metersConvToUnit = 1;
		unitConvToMeters = 1;

		# find center point of section in lat/lon coordinates to determine a UTM
		# zone to use for a temporary projected coordinate system for position coordinates
		# in the COLLADA file
		class TRANS2D_MAPGEN transObjToMap;
		xsecgeoref.GetTransParm(transObjToMap, 0, xsecgeoref.GetCalibModel() );

		class RECT3D extents;
		XSEC.GetExtents(extents);

		class POINT2D center;
		center = extents.center;

		printf("center object coordinates: x = %.1f, y = %.1f\n", center.x, center.y);

		center = transObjToMap.ConvertPoint2DFwd(center);

		printf("center map coordinates: x = %.6f, y = %.6f\n", center.x, center.y);

		# variables for setting up the UTM CRS
		class SR_COORDREFSYS utmCRS;
		class SR_COORDSYS utmCoordSys;		# coordinate system
		class SR_COORDOPDEF projectDef;		# projection parameters and method
		class SR_COORDOPMETHOD method;
		class SR_DATUM datum;

		datum = xsecCRS.Datum;
		utmCoordSys.Assign("Projected2D_EN_m");	# assign projected coordinate system
		method.Assign("1909");					# assign using numeric ID for Transverse Mercator projection
		projectDef.Method = method;		# assign method to projection

		# assign projection parameters that are the same for all UTM zones
		projectDef.SetParmValueNum("LatitudeOfNaturalOrigin", 0);
		projectDef.SetParmValueNum("ScaleFactorAtNaturalOrigin", 0.9996);
		projectDef.SetParmValueNum("FalseEasting", 500000);

		# assign false northing based on whether north or south of equator
		if ( center.y < 0) then
			projectDef.SetParmValueNum("FalseNorthing", 10000);
		else
			projectDef.SetParmValueNum("FalseNorthing", 0);

		# assign longitude of natural origin (central meridian of UTM zone) based on longitude
		numeric zoneNum, centMerid;		# UTM zone and longitude of central meridian
		zoneNum = ceil( (180 + center.x) / 6 );
		centMerid = (zoneNum * 6) - 183;
		projectDef.SetParmValueNum("LongitudeOfNaturalOrigin", centMerid);

		# create defined UTM coordinate reference system
		utmCRS.Create(utmCoordSys, datum, projectDef);

		# set up coordinate transformation from section's geographic CRS to the defined UTM CRS
		class TRANS2D_MAPGEN transGeogUTM;
		transGeogUTM.InputCoordRefSys = xsecCRS;
		transGeogUTM.OutputCoordRefSys = utmCRS;
		}

	# set up coordinate transformation from section CRS to geoCRS
	class TRANS2D_MAPGEN transToGeo;
	transToGeo.InputCoordRefSys = xsecCRS;
	transToGeo.OutputCoordRefSys = geoCRS;

	# report dimensions of the cross-section raster
	printf("Cross section %d dimensions: %d lines, %d columns\n", k, XSEC.$Info.NumLins, XSEC.$Info.NumCols);


	###############################################################
	#################  make a PNG file from the cross-section raster 

	# PIPELINE SOURCE
	class IMAGE_PIPELINE_SOURCE_RVC sourceXS(xsecObjItemList[k]);
	err = sourceXS.Initialize();
	if ( err < 0) ReportError(_context.CurrentLineNum, err);

	# PIPELINE FILTER NULL TO ALPHA
	# converts null mask to alpha channel for output PNG
	class IMAGE_PIPELINE_FILTER_NULLTOALPHA alpha(sourceXS);
	err = alpha.Initialize();
	if ( err < 0) ReportError(_context.CurrentLineNum, err);

	# PIPELINE TARGET
	class STRING png$ = dir$ + "/" + xsecName$ + ".png";
	class FILEPATH pngFilepath(png$);
	printf("PNG filepath = %s\n", png$);

	class IMAGE_PIPELINE_TARGET_PNG targetPNG(alpha, pngFilepath, "none");
	err = targetPNG.Initialize();
	if ( err < 0) ReportError(_context.CurrentLineNum, err); 

	# Execute the pipeline
	err = targetPNG.Process();
	if ( err < 0) ReportError(_context.CurrentLineNum, err);

	#####################################################################
	############### process manifold georeference control points

	# get the list of control points as an array of class CTRLPOINT
	printf("Getting list of georeference control points for section %d.\n", k);
	numeric numPts = xsecgeoref.GetNumCtrlPoints(); 
	class CTRLPOINT ctrlPtList[numPts];
	xsecgeoref.GetPointList(ctrlPtList);
	numeric numBndry = xsecgeoref.GetNumBoundaries();
	numeric numEdges = xsecgeoref.GetNumCtrlEdges();

	printf("Number of control points = %d\n", numPts);
	printf("Number of control edges = %d\n",  numEdges);
	printf("Number of triangulation boundaries = %d\n\n", numBndry );

	###### get boundary of triangulated control points if used (may be more than one boundary)
	print("Get control point boundary if present and list indices:");
	numeric numBndryPts;		# number of points in current boundary

	if (numBndry > 0)
		{
		if (bndryList.GetNumBoundaries() > 0) then
			bndryList.ClearBoundaries();

		for j = 1 to numBndry
			{
			numBndryPts = xsecgeoref.ReadBoundary(j-1, bndryArray);	# indexing of boundaries is 0-based
			printf("Number of points in boundary %d = %d\n", j, numBndryPts);

			for i = 1 to numBndryPts {
				printf("Boundary index %d = %d\n", i, bndryArray[i] + 1);
				}

			printf("\nAdd control point boundary to boundary list.\n");
			bndryList.AddBoundary(bndryArray, numBndryPts);
			ResizeArrayClear(bndryArray, 1);
			}
		}

	printf("Check boundary list for section %d:\n", k);
	printf("Number of boundaries in boundary list = %d\n", bndryList.GetNumBoundaries() );

	# resize and clear array holding all map coordinates (from destination points) and
	# texture coordinates (source points).
	ResizeArrayClear(xyz_array, numPts * 3);
	ResizeArrayClear(uv_array, numPts * 2);

	# resize and clear arrays for x and y coordinates used to determine coordinate ranges and
	# compute local 3D position coordinates for the COLLADA model
	ResizeArrayClear(xarray, numPts);
	ResizeArrayClear(yarray, numPts);
	ResizeArrayClear(zarray, numPts);

	minX = 100000000;	minY = 100000000;		# starting values greater than any anticipated real values


	# loop through points to read the destination (map) coordinates to arrays and determine
	# minimum values for x and y to use to create offset coordinates for each axis of
	# the COLLADA model; 
	for i = 1 to numPts
		{
		printf("Control Pt %d:\n", i);
		# get map point from current control point
		destPt = ctrlPtList[i].dest;

		if (geogCRSflag)  	# if input CRS is geographic, convert to our virtual UTM coordinates
			{
			destPt = transGeogUTM.ConvertPoint3DFwd(destPt);
			destPt.z = ctrlPtList[i].dest.z;
			}

		xarray[i] = destPt.x;
		yarray[i] = destPt.y;
		zarray[i] = destPt.z * metersConvToUnit;  # apply scaling from meters to coordinate system unit

		printf("x = %.6f, y = %.6f\n", xarray[i], yarray[i]);

		if (xarray[i] < minX) then minX = xarray[i];
		if (yarray[i] < minY) then minY = yarray[i];
		}

	printf("\nMinimum coordinate values for offset: x = %.4f, y = %.4f\n", minX, minY);

	# set min x, min y as origin point for the model; find its coordinates in WGS84 Geographic CRS and
	# find angle to north at this location in original CRS to set position and orientation of model in Google Earth
	originPt.x = minX;  originPt.y = minY;		# coordinates in CRS

	# if input geographic CRS, convert UTM origin point back to original CRS
	if (geogCRSflag) then originPt = transGeogUTM.ConvertPoint3DInv(originPt);

	angleToNorth = xsecCRS.ComputeAngleToNorth(originPt);
	originPt = transToGeo.ConvertPoint3DFwd(originPt);		# reproject to WGS 84 / Geographic 

	# set Z value for origin point in meters to provide correct vertical position for model in Google Earth
	originPt.z = minZm;	# minZ in meters for all sections

	printf("origin point: angle to north = %.1f, lon = %.6f, lat = %.6f z = %.1f m\n\n", angleToNorth, originPt.x, originPt.y, originPt.z);

	# loop through points again to write offset XYZ values to xyzString, read image coordinates and
	# write to uvString, and add original points to triangulator to determine mapping of vertices to
	# triangles for position mesh in COLLADA file
	xyzString.Clear();
	uvString.Clear();
	pString.Clear();
	triangulator.ClearAll();
	minZ = minZm * metersConvToUnit;		# origin elevation in CRS units as offset for z

	print("Create strings with vertex coordinates for position and texture meshes in COLLADA file:"); 
	for i = 1 to numPts
		{
		# add offset 3D map coordinates to the xyzString
		xyzString += sprintf("%.6f %.6f %.6f ", xarray[i] - minX, yarray[i] - minY, zarray[i] - minZ); 
		printf("xyzString = %s\n", xyzString);

		# read image coordinates
		col = ctrlPtList[i].source.x;		lin = ctrlPtList[i].source.y;
		# convert to texture coordinates (normalized to 0 to 1 range, origin in lower left) 
		s = col / XSEC.$Info.NumCols;
		t = 1 - (lin / XSEC.$Info.NumLins);		# inverted due to lower left origin

		# add to uvString
		printf("Image coordinates: col = %.6f, lin = %.6f\n", col, lin);
		printf("Texture coordinates: s = %.6f, t = %.6f\n", s, t);
		uvString += sprintf("%.6f %.6f ", s, t);
		printf("uvString = %s\n", uvString);

		triangulator.AddPoint(ctrlPtList[i].source);
		}

	# add hard edges from georeference (if any) to the triangulator
	printf("\nAdd hard edges to triangulator:\n");

	if (numEdges > 0)
		{
		for i = 1 to xsecgeoref.GetNumCtrlEdges()
			{ 
			xsecgeoref.ReadCtrlEdge(i-1, hardedge);	# internal indexing of hard edges is 0-based, so use i-1

			# internal indexing of hard edge points is 0-based, so add 1 to correspond to 1-based control point numbering
			printf("Hardedge %d: start = %d, end = %d\n", i-1, hardedge.start + 1, hardedge.end + 1);
			triangulator.AddHardEdgeAsIndices(hardedge.start + 1, hardedge.end + 1);
			}
		}

	# triangulate the control points
	err = triangulator.Triangulate(0, "ShortestEdge");
	printf("\nTriangulating, triangulation returned %d\n", err); 

	printf("Number of nodes = %d\n", triangulator.GetNumNodes() );
	printf("Number of edges = %d\n", triangulator.GetNumEdges() );
	printf("Number of triangles = %d\n", triangulator.GetNumTriangles() );

	## add control point boundary read from the georeference object to the triangulator.
	## This must be done AFTER the points have been triangulated; points and triangles outside
	## the boundary are then discarded.
	if (numBndry > 0) 
		{
		printf("\nAdd control point boundary to triangulator.\n");
		triangulator.SetBoundaries(bndryList);

		# check number of elements again
		printf("Number of nodes = %d\n", triangulator.GetNumNodes() );
		printf("Number of edges = %d\n", triangulator.GetNumEdges() );
		printf("Number of triangles = %d\n\n", triangulator.GetNumTriangles() );
		}

	# The triangulator creates a new set of Node Numbers but provides a GetIndex() method
	# to find the node index (1-based) which is set by the order in which the points were added
	# to the triangulator.  Node indices therefore map to the georeference control point number.  
	# Hard edges added to the set result in additional nodes (each at the
	# same location as one of the original points). 
	print("Loop through nodes to report indices and coordinates.");
	for i = 1 to triangulator.GetNumNodes()
		{
		triangulator.GetPoint(i, pt);
		printf("Node %d, index = %d, x = %.6f, y = %.6f\n", i, triangulator.GetIndex(i), pt.x, pt.y);
		}

	numTriangles = triangulator.GetNumTriangles();

	# loop through the triangles to get the nodes for each and find the corresponding georeference control points
	for i = 1 to numTriangles
		{
		err = triangulator.GetTriangleNodes(i, triNodes);
		printf("\nTriangle %d, GetTriangleNodes returned %d\n", i, err);

		for j = 1 to 3 	# loop through nodes of triangle to get the control point number (node index) for each node
			{
			nodeNum = triNodes.GetNode(j);
			nodeIndex = triangulator.GetIndex(nodeNum);

			if (nodeIndex > numPts)	# check if node was added with a hard edge and if so find the
													# corresponding node from the original control point set
				{
				triangulator.GetPoint(nodeNum, nodePt);
				nodeNum = triangulator.FindDifferentNode(nodePt, nodeNum);
				nodeIndex = triangulator.GetIndex(nodeNum);
				}

			printf("Node %d = %d, Point Index = %d\n", j, nodeNum, nodeIndex);

			tempArray[j] = nodeIndex;	# add node index to the array of triangle node indices
			}

		printf("Temp array indices: %d, %d, %d\n", tempArray[1], tempArray[2], tempArray[3]);

		testCCW = ccwImageCoords(ctrlPtList[ tempArray[1] ].source,	ctrlPtList[ tempArray[2] ].source, ctrlPtList[ tempArray[3] ].source);
		printf("counter-clockwise test returned %.4f\n", testCCW);

		triNodeIndexArray[1] = tempArray[1];

		if (testCCW > 0) 	# if nodes not in counterclockwise order, reorder nodes 2 and 3
			{
			triNodeIndexArray[2] = tempArray[3];
			triNodeIndexArray[3] = tempArray[2];
			}
		else
			{
			triNodeIndexArray[2] = tempArray[2];
			triNodeIndexArray[3] = tempArray[3];
			}

		printf("Triangle node indices in counter-clockwise order: %d, %d, %d\n", triNodeIndexArray[1], triNodeIndexArray[2], triNodeIndexArray[3]);

		# loop through the nodes for the current triangle to write index positions to the "primitive" string (pString)
		# for the COLLADA file.  For each triangle need three values, in order: index for vertex, index for normals,
		# index for texture. All indices are 0-based; index for vertex and texture are equal to each other and equal
		# node index - 1.
		for j = 1 to 3
			{
			pString += sprintf("%d %d %d ", triNodeIndexArray[j] - 1, 0, triNodeIndexArray[j] - 1);
			}

		printf("pString = %s\n\n", pString);
		}

	# close the georeference object and its parent raster object
	xsecgeoref.Close();
	XSEC.Close();

	# get current date and time for inclusion in the COLLADA file
	dt.SetCurrent();
	format$ = "%Y-%m-%dT%H:%M:%SZ";
	dt$ = dt.GetString(format$);
	printf("formatted date and time: %s\n", dt$);


	#######################################################################
	# parse the onesection.dae COLLADA file template and get the root node
	print("Making COLLADA file.");
	err = daedoc.Parse(dae$);
	if (err < 0) {
		PopupError(err);
		break; 
		}

	colladaRootNode = daedoc.GetRootNode();

	# set created and modified dates, units, and scale factor to meters
	node = colladaRootNode.FindChild("asset");
	node1 = node.FindChild("created");
	node1.SetText(dt$);
	node1 = node.FindChild("modified");
	node1.SetText(dt$);

	node1 = node.FindChild("unit");
	node1.SetAttribute("name", unit$);
	node1.SetAttribute("meter", sprintf("%.4f", unitConvToMeters) );

	### set the output PNG file as the library image
	node = daedoc.GetElementByID("Xsect-image");
	node1 = node.FindChild("init_from");
	node1.SetText(pngFilepath.GetName() );

	### set the position mesh count and array
	arraycount = numPts * 3;
	node = daedoc.GetElementByID("mesh1-geometry-position-array");
	node.SetAttribute("count", NumToStr(arraycount) );
	node.SetText(xyzString);

	### set the count for the position array accessor
	node1 = node.GetParent();
	node2 = node1.FindChild("technique_common");
	node3 = node2.FindChild("accessor");
	node3.SetAttribute("count", NumToStr(numPts) );

	### set the texture mesh count and array
	arraycount = numPts * 2;
	node = daedoc.GetElementByID("mesh1-geometry-uv-array");
	node.SetAttribute("count", NumToStr(arraycount) );
	node.SetText(uvString);

	### set the count for the uv array accessor
	node1 = node.GetParent();
	node2 = node1.FindChild("technique_common");
	node3 = node2.FindChild("accessor");
	node3.SetAttribute("count", NumToStr(numPts) );

	### set the material and primitive for the triangles element
	node = daedoc.GetElementByID("mesh1-geometry");
	node1 = node.FindChild("mesh");
	node2 = node1.FindChild("triangles");
	node2.SetAttribute("count", NumToStr(numTriangles) );

	node3 = node2.FindChild("p");
	node3.SetText(pString);

	############## write out the COLLADA file
	class STRING daefile$ = dir$ + "/" + xsecName$ + ".dae";
	printf("\nDAE filepath = %s\n", daefile$);
	daedoc.Write(daefile$);

	class FILEPATH daeFilepath(daefile$);


	##################################################################
	#############
	############# add the section PNG file and COLLADA file to the KMZ
	kmz.AddFile(kmzFilepath, pngFilepath, "files/" + xsecName$ + ".png");
	kmz.AddFile(kmzFilepath, daeFilepath, "files/" + xsecName$ + ".dae");

	# delete the source files
	DeleteFile(png$);
	DeleteFile(daefile$);


	##################################################################
	############
	############  add info to the KML structure for this section
	startDate = 2000;

	document = folder.NewChild("Document");
	node = document.NewChild("name");
	node.SetText(xsecName$);

	# create cross-section placemark under document for each extrusion step
	for i = 1 to numExtrudeSteps + 1
		{
		placemark = document.NewChild("Placemark");

		# create name element for placemark
		placemark.NewChild("name", sprintf("%s %d", xsecName$, i) );

		# create time span element for placemark
		node = placemark.NewChild("TimeSpan");

		endDate = startDate + 100;
		node.NewChild("begin", NumToStr(startDate) );
		node.NewChild("end", NumToStr(endDate) );
		startDate = endDate; 	# increment start date for next loop

		# create style element for placemark
		node = placemark.NewChild("Style");
		node.SetAttribute("id", "default");

		# create model element for placemark
		model = placemark.NewChild("Model");
		model.SetAttribute("id", "model");

		# create altitude mode element for model
		model.NewChild("altitudeMode", "absolute");

		# create location element for model
		node = model.NewChild("Location");
		node.NewChild("longitude", NumToStr(originPt.x) );
		node.NewChild("latitude", NumToStr(originPt.y) );
		node.NewChild("altitude", NumToStr(originPt.z + (i-1) * altitudeIncrement ) );

		printf("Extrude step %d: altitude = %.1f m\n", i - 1, originPt.z + (i-1) * altitudeIncrement);

		# create orientation element for model
		node = model.NewChild("Orientation");
		node.NewChild("heading", NumToStr(angleToNorth) );
		node.NewChild("tilt", "0");
		node.NewChild("roll", "0");

		# create scale element for model
		node = model.NewChild("Scale");
		node.NewChild("x", "1.0");
		node.NewChild("y", "1.0");
		node.NewChild("z", "1.0");

		# create link element for model
		node = model.NewChild("Link");
		node.NewChild("href", "files/" + xsecName$ + ".dae");

		# create ResourceMap element for model
		node = model.NewChild("ResourceMap");
		node1 = node.NewChild("Alias");

		node1.NewChild("targetHref", "files/" + xsecName$ +".png");
		node1.NewChild("sourceHref", "files/" + xsecName$ +".png");
		}
	}	# end:  for k = 1 to validcount



############## Write the KML structure to a file and add to KMZ
kmldoc.Write(kmlfile$);
kmz.AddFile(kmzFilepath, kmlFilepath, "doc.kml");
DeleteFile(kmlfile$);


print("End.");