# 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$ = '<?xml version="1.0" encoding="utf-8"?>
<COLLADA xmlns="http://www.collada.org/2005/11/COLLADASchema" version="1.4.1">
<asset>
<contributor>
<authoring_tool>TNTmips SML Script</authoring_tool>
</contributor>
<created></created>
<modified></modified>
<unit name="" meter="1.0"/>
<up_axis>Z_UP</up_axis>
</asset>
<library_images>
<image id="Xsect-image" name="Xsect-image">
<init_from></init_from>
</image>
</library_images>
<library_materials>
<material id="XsectID" name="Xsect">
<instance_effect url="#Xsect-effect"/>
</material>
</library_materials>
<library_effects>
<effect id="Xsect-effect" name="Xsect-effect">
<profile_COMMON>
<newparam sid="Xsect-image-surface">
<surface type="2D">
<init_from>Xsect-image</init_from>
</surface>
</newparam>
<newparam sid="Xsect-image-sampler">
<sampler2D>
<source>Xsect-image-surface</source>
</sampler2D>
</newparam>
<technique sid="COMMON">
<phong>
<emission>
<color>0.000000 0.000000 0.000000 1</color>
</emission>
<ambient>
<color>0.000000 0.000000 0.000000 1</color>
</ambient>
<diffuse>
<texture texture="Xsect-image-sampler" texcoord="UVSET0"/>
</diffuse>
<specular>
<color>0.330000 0.330000 0.330000 1</color>
</specular>
<shininess>
<float>20.000000</float>
</shininess>
<reflectivity>
<float>0.100000</float>
</reflectivity>
<transparent>
<color>1 1 1 1</color>
</transparent>
<transparency>
<float>0.000000</float>
</transparency>
</phong>
</technique>
<extra>
<technique profile="GOOGLEEARTH">
<double_sided>1</double_sided>
</technique>
</extra>
</profile_COMMON>
</effect>
</library_effects>
<library_geometries>
<geometry id="mesh1-geometry" name="mesh1-geometry">
<mesh>
<source id="mesh1-geometry-position">
<float_array id="mesh1-geometry-position-array" count="">
</float_array>
<technique_common>
<accessor source="#mesh1-geometry-position-array" count="" stride="3">
<param name="X" type="float"/>
<param name="Y" type="float"/>
<param name="Z" type="float"/>
</accessor>
</technique_common>
</source>
<source id="mesh1-geometry-normal">
<float_array id="mesh1-geometry-normal-array" count="3">-0.000000 -1.000000 -0.000000 </float_array>
<technique_common>
<accessor source="#mesh1-geometry-normal-array" count="1" stride="3">
<param name="X" type="float"/>
<param name="Y" type="float"/>
<param name="Z" type="float"/>
</accessor>
</technique_common>
</source>
<source id="mesh1-geometry-uv">
<float_array id="mesh1-geometry-uv-array" count="">
</float_array>
<technique_common>
<accessor source="#mesh1-geometry-uv-array" count="" stride="2">
<param name="S" type="float"/>
<param name="T" type="float"/>
</accessor>
</technique_common>
</source>
<vertices id="mesh1-geometry-vertex">
<input semantic="POSITION" source="#mesh1-geometry-position"/>
</vertices>
<triangles material="Xsect" count="">
<input semantic="VERTEX" source="#mesh1-geometry-vertex" offset="0"/>
<input semantic="NORMAL" source="#mesh1-geometry-normal" offset="1"/>
<input semantic="TEXCOORD" source="#mesh1-geometry-uv" offset="2" set="0"/>
<p> </p>
</triangles>
</mesh>
</geometry>
</library_geometries>
<library_visual_scenes>
<visual_scene id="CrossSectionScene" name="CrossSectionScene">
<node id="Model" name="Model">
<node id="mesh1" name="mesh1">
<instance_geometry url="#mesh1-geometry">
<bind_material>
<technique_common>
<instance_material symbol="Xsect" target="#XsectID">
<bind_vertex_input semantic="UVSET0" input_semantic="TEXCOORD" input_set="0"/>
</instance_material>
</technique_common>
</bind_material>
</instance_geometry>
</node>
</node>
</visual_scene>
</library_visual_scenes>
<scene>
<instance_visual_scene url="#CrossSectionScene"/>
</scene>
</COLLADA>
';
# 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$ = '<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom">
<Folder>
<name></name>
<open>1</open>
</Folder>
</kml>';
# 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.");
Download