StrikeDipTool76.sml

  Download

More scripts: Display Toolbar

Syntax Highlighing:

comments, key words, predefined symbols, class members & methods, functions & classes
            
##################################################################################
#
# StrikeDipTool76.sml (Tool Script)
#
##################################################################################
#
# Started: 3 November 2003
# Updated: 30 April 2013.
#
# Authors:
#   Dave Breitwisch
#   Randy Smith
#   Dan Glasser
#
#	Version 29 April 2013
#	   Added code to the OnModeChange, OnAccept, and OnCancel functions to properly  
#		manage the "ScriptGadget" and the polygon tool. Fixes a problem with
#		selecting an existing point in Edit and Delete modes.
#
#	Version 19 July 2011
#		Changed the function to create the polyline tool to ViewCreatePolyLineGadget()
#		to be consistent with internal class changes.  Also assigns initial extents
#		to the new vector objects created from the group extents.
#
#	Version 29 April 2008
#		The CartoScript for the Strike/Dip vector is now read from a string variable
#			included within this ToolScript; the display parameters for this vector
#			are also now automatically saved so that the CartoScript remains associated
#			with the saved Strike/Dip vector.
#		If the DEM Coordinate Reference System has a projected coordinate system with
#			units not in meters, elevation values from the DEM are now automatically scaled
#			to that unit to ensure that the proper dip value is computed.
#		The Bedding table now records the Z value for each of the 3 points used to
#			compute the strike and dip of that station.  The unit of measure for the recorded
#			coordinates is included in the description of the table.
#
#	Version 18 April 2008
#		Fixes a problem with reversing the strike direction when the overturned
#		toggle button is changed.
#
#	Version 24 March 2008
#		Fixes a problem with the user-defined function to return angle to north
#		that is used to correct the strike azimuth.
#
# 	Version Nov. 2007
#	Revisions:
#		Produces correct attitude values for DEMs with geographic coordinates by
#			using appropriate UTM zone coordinates for attitude computations.
#		Automatically applies raster scale value to allow scaled DEMs with values in feet or decimeters.
#		Strike azimuth and dip direction are automatically adjusted for the local
#			grid angle to north at each point.
##################################################################################
##################################################################################
#
# Purpose:
#
#     The purpose of this Tool Script is to give a geologist a streamlined method
#     for determining the Strike and Dip values from an accurate DEM raster and
#     overlaid arial imagery.  The accompanying CartoScript displays the added
#     points using the appropriate attitude symbols.  The script creates two
#     vector objects: a StrikeDipVector to store point elements with attached
#     attitude values in the associated point database, and an OutcropTrace
#     object to store outcrop traces for each point (see below).  The script creates
#     and displays a dialog to show the computed attitude values and allow
#     designation of overturned beds.
#
# Assumptions:
#
#		The DEM is the bottom layer in the active group.  The DEM elevation cell values
#		are in meters or the DEM raster's elevation scale is set so that the scaled
#		elevation values are meters.
#
# Description:
#
#     This Tool Script has four different modes that the geologist can use:
#
#        ADD - While viewing the reference image, place 3 points (three vertices of a
#              triangle) on a bedding plane by left clicking on the screen using the
#              provided polygon tool.  Adjust the vertices if needed using the
#              Line/Polygon Edit Controls.  Right-click or press [Apply] on the
#              Edit Controls:  x = easting and y = northing values are registered for
#              each point from the group map coordinates and z = elevation is registered
#              for each point from the corresponding cell of the elevation raster that
#              is the first layer in the group.  If the xyz values pass a validity check,
#              these three noncolinear points are used to compute the strike azimuth,
#              dip angle, and dip azimuth of the plane.  These values are shown in
#              the Strike-Dip Toolscript Controls window, which also includes a
#              toggle button to indicate if the bed is overturned.
#
#              The Tool Script also creates and displays an outcrop trace in the vicinity
#              of the attitude measurement.  This trace is the computed line of intersection
#              between the measured plane and the surface represented by the DEM.
#              Verify that the attitude values appear reasonable and that the computed
#              outcrop trace parallels the actual outcrop trace in the image.  (The
#              Line/Polygon tool remains active to allow you to adjust the points and
#              reapply to recompute a new attitude measurement and outcrop trace if
#              needed.)
#
#              When you are satisfied, press the Accept button on the Tool Script
#              Controls window.  A point is added to the StrikeDip vector object at the
#              center of the triangle and the attitude values are written to an attached
#              record in the associated point element database.  The point is styled
#              automatically with the appropriately-oriented and labeled attitude symbol
#              using the accompanying CartoScript. The CartoScript must be in the
#              same directory as this Tool Script (or, if the Tool Script is saved with
#              a Group or Layout, in the same directory as the Project File
#              containing that Group or Layout).  The outcrop trace line is also saved
#              in the separate OutcropTrace vector object along with a specified color.
#              The original three points are also saved in this color with the trace line.
#
#     DELETE - Select a previously-added Strike/Dip point and press [Accept] to delete it
#              along with the outcrop trace pattern that is associated with it.
#
#       EDIT - Select a previously-added Strike/Dip point and press [Accept] to
#              edit it.  The original triangle and outcrop trace are reloaded to
#              enable vertex locations to be moved using the Line/Polygon tool as
#              in Add mode; Apply or right-click to recompute a new attitude and
#              outcrop tract.  Press [Accept] on the Tool Script Control window to
#              save the edited point and outcrop trace.
#
#       VIEW - Select a Strike/Dip point and use the Accept button to toggle whether
#              the outcrop trace pattern associated with that point is displayed.
#              Both the trace pattern and the three points are displayed in their
#              saved color.
#
# Definitions:
#
#     Strike Angle = azimuth of the line of intersection of the plane with the horizontal.
#          Note: The Strike Angle value conforms to the right hand rule, where the strike
#                angle is in the direction along the strike line for which the dip
#                direction is to the right.
#
#     Dip Angle = maximum vertical angle of the plane measured downward from the horizontal.
#
#     Dip Direction = azimuth of the dip vector, by definition is perpendicular to Strike.
#
#     Linear equation of the plane = a*x + b*y c*z + d = 0
#
#          Note: The coordinates of the three selected points are used to compute the
#                coefficient values of a, b and c.  d is then solved mathematically.
#
##################################################################################
##################################################################################
#
# Standard ToolScript Comments:
#
# The following symbols are predefined
#    class GRE_VIEW View            {use to access the view the tool script is attached to}
#    class GRE_GROUP Group          {use to access the group being viewed if the script is run from a group view}
#    class GRE_LAYOUT Layout        {use to access the layout being viewed if the script is run from a layout view}
#    numeric ToolIsActive           Will be 0 if tool is inactive or 1 if tool is active
#
# The following values are also predefined and are valid when the various On...()
# functions are called which deal with pointer and keyboard events.
#    numeric PointerX               Pointer X coordinate within view in pixels
#    numeric PointerY               Pointer Y coordinate within view in pixels
#    numeric ShiftPressed           1 if <shift> key being pressed or 0 if not
#    numeric CtrlPressed            1 if <ctrl> key being pressed or 0 if not
#    numeric LeftButtonPressed      1 if left pointer button pressed or 0 if not
#    numeric RightButtonPressed     1 if right pointer button pressed or 0 if not
#    numeric MiddleButtonPressed    1 if middle pointer button pressed or 0 if not
#
##################################################################################
##################################################################################
#
# User Defined Global Variable
#
# Increase this variable value to expand the size of the Outcrop Trace Pattern that
# is created when using the tool.
# Example:
# 	Final Left  Extents = Original Left  - (width of box * stretchClipBoundry)
# 	Final Right Extents = Original Right + (width of box * stretchClipBoundry)
# 	Final Upper Extents = Original Upper + (height of box * stretchClipBoundry)
# 	Final Lower Extents = Original Lower - (height of box * stretchClipBoundry)
#
##################################################################################
numeric stretchClipBoundry = 0.5;
##################################################################################
##################################################################################
# ToolScript Global Variables
##################################################################################
class XMLDOC dlgdoc;
class XMLNODE node;
class GUI_DLG dlg;
class GUI_CTRL_LABEL promptString;
class GUI_CTRL_TOGGLEBUTTON overturnedToggle;
class GUI_CTRL_EDIT_NUMBER strikeFld, dipFld, dipDirFld;
class GUI_CTRL_PUSHBUTTON btnACCEPT, btnCANCEL, btnDELETE, btnCLOSE;
class GUI_FORM_RADIOGROUP rGpMODE;
class GUI_CTRL_LABEL aLINECOLOR_LABEL;
class GUI_CTRL_COLORBUTTON btnLINECOLOR;
class GUI_GADGET_POLYLINE tool;
class POLYLINE poly;
class RVC_RASTER DEM, PLANE, clipRast;
class VECTOR StrikeDipVector, OutcropPVector, intersectV, clipV;
class GEOREF SDgeoref, demGeoref;
class RVC_GEOREFERENCE rvcGeorefDEM;
class TRANSMODEL demCalibModel;
class TRANS2D_MAPGEN transDEMobjToMap;
class SR_COORDREFSYS demCRS, utmCRS;
class TRANS2D_MAPGEN transGeogUTM;	# coordinate transformation from geographic to UTM coordinates
class GRE_GROUP activegroup, firstGroup;
class GRE_LAYER activeLayer;
class GRE_LAYER_RASTER rasterLayer, intersectRlayer;
class GRE_LAYER_VECTOR SDvectorLayer, OPvectorLayer, intersectVlayer, clipVlayer;
class GRE_VECTOR_POINTS SDpoints, OPpoints;
class GRE_VECTOR_LINES OPlines;
class GRE_VECTOR_POLYS OPpolys;
class DATABASE SDdb, OPlinedb, OPpointdb, iVpointdb, iVlinedb;
class DBTABLEINFO SDtable, OPlinetable, OPpointtable, iVpointtable, iVlinetable;
class DBFIELDINFO SDstationfield;
numeric demZscale;				# z-scale value for DEM
numeric unitZscale;				# additional scale factor for z if units for projected CRS are not meters
numeric a=0, b=0, c=0, d=0;	# linear coefficients for the equation of the plane
numeric dipX, dipY;           # partial derivatives of dip (negative slope) in x and y directions
numeric dipAng;               # magnitude of dip angle
numeric dipDir;               # azimuth of dip direction (0 to 360 degrees clockwise from north)
numeric strikeAng;            # azimuth of strike line
numeric station, overturned;
numeric setDefaultWhenClose;
numeric selectedElementNum;
numeric currentlyEditing = 0;
numeric waitingForAddConfirmation = 0;
numeric deletedStation;
numeric transGeogCoords;		# flag to indicate map coordinates need to be transformed from
										# geographic to UTM
array numeric records[1];
class POINT3D ctrPt;
class STRING coordUnitName$;		# coordinate unit name for x, and y values recorded in Strike/Dip table
class STRING zCoordUnit$;			# coordinate unit for z values recorded in Strike/Dip table
class STRING Current_Mode = "ADD";
class STRING genericQuery = "(Outcrop.Station == ";
class STRING currentQuery;
##################################################################################
# CartoScript used to draw the strike/dip symbols; script is automatically saved
# with the StrikeDipVector's display parameters.
class STRING sdQry$ = '# beddingStrikeDip.qry
# This is a sample script for drawing geologic point
# symbols for strike and dip of bedding, including
# special cases of horizontal, vertical, overturned bedding.
# The dip value is shown as a text label, except in the
# cases of horizontal and vertical bedding.
# Information required for each point includes the strike
# value (degrees), dip value (degrees), and whether the bed is
# overturned.  Each of these values is read from a particular
# field in an attached database table.  To use the script you
# will need to edit the statements that include database references
# to conform to your table name and field names.
# Strike angle must be specified as azimuth (0 to 360 degrees
# clockwise from north).  Dip direction is defined by the
# strike azimuth using the right hand rule: of the two possible
# azimuths for a given strike line, choose the azimuth the
# strike line points toward with the dip direction on the right hand
# side of the strike line. (For overturned beds, substitute facing
# direction [stratigraphic top] for dip direction in the rule).
# The script assumes that the database field for overturned
# bedding is a Logical field, with Yes indicating overturned bed
# and No (default value) indicating upright bedding.
# Modified for Legend samples August 2002.
# Modified to declare all variables October 2005;
# Modified to correct for angle between grid north and true north May 2006.
# Version Nov. 2007
# Requires TNTmips 2007:73 or later
# Modified to adjust scale based on georeference map units.
numeric azStrike, dip1, overturned1, red, green, blue;
numeric scale, strikeLengthMap, lineWidthMap, strikeLength, lineWidth;
numeric halfLength, tickLength, doubTick, halfTick, heightMap, height, offset;
string fontName$, label$;
numeric direction, oppStrike, dipDir, oppDip;
numeric next_x ,next_y, labelLength, shift1, shift2;
class RVC_GEOREFERENCE vGeoref;
class SR_COORDREFSYS vectCRS;
numeric northAng;
class POINT2D point;
string coordUnit$;
###################### Set Parameters ##############################
# Read strike azimuth and dip value from table.field
azStrike = Bedding.StrikeAngle			#### change to fit your data ####
dip1 = Bedding.DipAngle;					#### change to fit your data ####
# This variable defines the denominator of the intended map scale.
# It is used as the basis for defining line width and symbol size.
# Example: for 1:24,000 map scale, Scale = 24000
scale = CurrentMapScale * 1.8;
# Check if vector has geographic coordinates (units of degrees instead of meters)
# and if so adjust scale factor to draw symbols of appropriate size.
Vect.GetDefaultGeoref(vGeoref);
vectCRS = vGeoref.GetCoordRefSys();
if (vectCRS.IsProjected() <> 1) {
	if (vectCRS.IsLocal() <> 1) {
		scale = scale * 0.00001;
		}
	}
else {	# CRS is projected; check coordinate units to adjust scale
	# get coordinate unit from the first axis of the planar coordinate system
	coordUnit$ = vectCRS.Coordsys.GetAxis(1).Unit.GetSymbol();
	scale = scale * GetUnitConvDist("m", coordUnit$);
	}
# find angle to north at the current point\'s map coordinates
class TRANS2D_MAPGEN transObjToMap;
class TRANSMODEL model;
model = vGeoref.GetCalibModel();
vGeoref.GetTransParm(transObjToMap, 0, model);
point.x = Vect.Point.Internal.x;		# object coordinates of current point
point.y = Vect.Point.Internal.y;
point = transObjToMap.ConvertPoint2DFwd(point);		# convert to map coordinates
northAng = vectCRS.ComputeAngleToNorth(point);
# subtract angle to north from stored strike value for display
azStrike = azStrike - northAng;
# red, green, blue variables define the color of the symbols and label
red = 0;			green = 0;			blue = 0;
# These variables define the dimensions and line widths of the
# symbol.  StrikeLengthMap is the desired length of the symbol
# strike line in mm, assuming vector coordinates are in meters.
# LineWidthMap is the desired line width in mm.
strikeLengthMap = 6;
lineWidthMap = 0.3;
if (DrawingLegendView == 1) {
	strikeLength = 0.5 * SampleRect.GetWidth();
	lineWidth = 0.08 * strikeLength;
	}
else {
	strikeLength = strikeLengthMap * scale / 1000;
	lineWidth = lineWidthMap * scale / 1000;
	}
halfLength = 0.5 * strikeLength;
tickLength = halfLength * 0.5;
halfTick = tickLength * 0.5;
doubTick = tickLength * 2;
# These variables control the drawing of the label text string.
# HeightMap is the desired height of the letters in mm,
# assuming vector coordinates are in meters.
heightMap = 2.5;
height = heightMap * scale / 1000;
fontName$ = "ARIALBD.TTF";
offset = 0.5 * tickLength; 		# offset of label from symbol
######################### Process ############################
# Convert strike azimuth to internal coordinate system.
direction = -(azStrike - 90);
if (direction < 0)
	direction = direction + 360;
oppStrike = direction -180;
dipDir = direction - 90;
oppDip = dipDir - 180;
# Convert dip value to a string and assign it
# to a string variable for use as label
label$ = sprintf("%2d", dip1);
# Find length of label text for label positioning
LineStyleTextNextPosition(label$, height, 0, 1, next_x ,next_y, labelLength);
# Set line color, width, and end type
LineStyleSetColor(red, green, blue);		# set symbol color
LineStyleSetLineWidth(lineWidth);
LineStyleSetCapJoinType(1,1); 	# square ends of lines
########### Draw symbol
# Special symbol for horizontal bedding (cross in circle)
if (dip1 == 0){
	LineStyleDropAnchor(0);
	LineStyleDrawCircle(halfLength);
	LineStyleMoveTo(90, halfLength);
	LineStyleLineTo(-90, strikeLength);
	LineStyleMoveToAnchor(0);
	LineStyleMoveTo(0, halfLength);
	LineStyleLineTo(180, strikeLength);
	}
else {
	# For nonzero dip, draw strike line with center at point
	LineStyleDropAnchor(0);
	LineStyleMoveTo(direction, halfLength);
	LineStyleLineTo(oppStrike, strikeLength);
	LineStyleMoveToAnchor(0);
	# Draw appropriate symbol for dip direction
	if (dip1 == 90) {				# crossbar for vertical bed
		LineStyleMoveTo(dipDir, tickLength);
		LineStyleLineTo(oppDip, doubTick);
		}
	else {
		if (Bedding.Overturned) {		# dip symbol for overturned beds
			LineStyleDrawArc(0, 0, halfTick, halfTick, direction, -180, 0);
			LineStyleMoveTo(direction, halfTick);
			LineStyleLineTo(oppDip, doubTick);
			}
		else {						# dip direction tick mark
			LineStyleLineTo(dipDir, tickLength);
			}
		}
	}
########## Add dip label if bedding is not horizontal or vertical
if (dip1 > 0 and dip1 < 90) {
	# Set font name and text color
	LineStyleSetFont(fontName$);
	LineStyleSetTextColor(red, green, blue);
	# Compute label shift parallel to strike to center label
	shift1 = 0.5 * ( labelLength * cosd(direction) + height * sind(direction) );
	if (shift1 <0) {
		shift1 = abs(shift1);
		LineStyleMoveTo(direction, shift1);
		}
		else {
		LineStyleMoveTo(oppStrike, shift1);
		}
	# Compute label shift in dip direction to avoid overwriting symbol
	if (Bedding.Overturned) {
		if (direction >= 0 and direction < 90)
			shift2 = offset + 0.8 * labelLength * cosd(dipDir);
		if (direction >= 90 and direction < 180)
			shift2 = offset + 0.8 * (labelLength * cosd(dipDir) + height * sind(dipDir));
		if (direction >= 180 and direction < 270)
			shift2 = offset + 0.8 * height * sind(dipDir);
		if (direction >= 270 and direction <= 360)
			shift2 = offset;
		LineStyleMoveTo(direction + 90, shift2);
		LineStyleDrawText(label$, height, 0, 1);
		}
	else {			# place label for upright bedding
		if (direction >= 0 and direction< 90)
			shift2 = offset - 0.8 * height * sind(dipDir);
		if (direction >= 90 and direction < 180)
			shift2 = offset;
		if (direction >= 180 and direction < 270)
			shift2 = offset - 0.8 * labelLength * cosd(dipDir);
		if (direction >= 270 and direction <= 360)
			shift2 = offset - 0.8 * ((labelLength * cosd(dipDir)) + (height * sind(dipDir)));
		LineStyleMoveTo(dipDir, shift2);
		LineStyleDrawText(label$, height, 0, 1);
		}
	}';
##################################################################################
# Function Prototypes
##################################################################################
# Function prototypes with a class as a return value are unallowed at this time
##################################################################################
#func class POINT3D setXY(numeric num);
#func class POINT3D setZ(class POINT3D pt3d, numeric num);
#func class POINT3D computeDifference(class POINT3D pt1, class POINT3D pt2);
#func class POINT3D computeCenterPoint(class POINT3D pt1, class POINT3D pt2, class POINT3D pt3);
#func class POINT3D computePlanarCoefficients();
# Functions used to compute the Strike and Dip values
##################################################################################
proc createPlaneRaster(class POINT3D testPoint, class POINT3D point2, class POINT3D point3);
func Determinant(a1, b1, a2, b2);
func isZValid(numeric z);
proc computeVerticalPlane();
proc computeNonVerticalPlane();
proc computeStrikeDip();
# Functions used to manage the Strike and Dip data
##################################################################################
func addPointToVector(class POINT3D pt);
proc addRecordToDatabase( numeric pointElemNumAdded );
proc UpdateQuery( );
proc deleteStation();
# Functions used to initiliaze or maintain the display layers and databases
##################################################################################
func checkLayer();
proc checkGeoref();
proc createDestVectors();
func createVectorGeoreference();
func vectorLayerExists(string name);
proc addVectorsToDisplay();
proc initializeBeddingDatabaseTable();
proc createBeddingFields();
proc initializeOutcropDatabaseTable();
# Functions used to manage callbacks
##################################################################################
proc cbToolApply();
proc cbLayer();
proc cbGroup();
proc cbClose();
# Functions called from interaction with the dialog
##################################################################################
proc OnModeChange( );
proc OnAccept( );
proc OnCancel( );
proc OnDelete();
proc OnClose();
proc OnOverturnedChanged();
# Basic ToolScript functions
##################################################################################
proc OnRightButtonPress();
proc OnLeftButtonPress();
proc OnInitialize ();
proc OnDestroy ();
proc OnActivate ();
proc OnDeactivate ();
##################################################################################
##################################################################################
# Function Definitions
##################################################################################
##################################################################################
#### Get the x and y screen coordinates from a polyline vertex
func class POINT2D setXY(numeric num) {
	local class POINT2D tmp = poly.GetVertex(num - 1);
	return tmp;
}
##################################################################################
#### Convert polyline vertex coordinates to DEM map coordinates. Get the z value
#### from the elevation raster at that location.
func class POINT3D setZ(class POINT2D pt) {
	# convert from screen to view coords
	pt = TransPoint2D(pt, ViewGetTransViewToScreen(View, 1));
	# convert from view to map coords
	pt = TransPoint2D(pt, ViewGetTransMapToView(View, rasterLayer.MapRegion.CoordRefSys, 1));
	# convert from map coords to object
	local class Georef georef;
	local class POINT2D objcoord;
	georef = GetLastUsedGeorefObject(DEM);
	objcoord = MapToObject(georef, pt.x, pt.y, DEM);
	# with the object coordinates we can get the cell value
	numeric cellval = DEM[objcoord.y, objcoord.x];
	# return a 3D point with map coordinates and elevation
	local class POINT3D retPnt;
	retPnt.x = pt.x;
	retPnt.y = pt.y;
	retPnt.z = cellval * demZscale * unitZscale; # a raster cell value with Z-scaling applied
	return retPnt;
}
##################################################################################
#### Calculate the difference between two points
func class POINT3D computeDifference(class POINT3D pt1, class POINT3D pt2) {
	local class POINT3D diff;
	diff.x = pt1.x - pt2.x;
	diff.y = pt1.y - pt2.y;
	diff.z = pt1.z - pt2.z;
	return diff;
}
##################################################################################
#### Compute the center point of the three points
func class POINT3D computeCenterPoint(class POINT3D pt1, class POINT3D pt2, class POINT3D pt3) {
	local class POINT3D ctr;
	ctr.x = (pt1.x + pt2.x + pt3.x) / 3;
	ctr.y = (pt1.y + pt2.y + pt3.y) / 3;
	ctr.z = (pt1.z + pt2.z + pt3.z) / 3;
	return ctr;
}
##################################################################################
#### Convert geographic coordinates to planar coordinates in the designated UTM CRS
func class POINT3D convertGeogToUTM(class POINT3D geogPt) {
	local class POINT2D mapCoords;	# 2D map coordinates
	local class POINT3D utmPt;			# 3D point to return
	mapCoords.x = geogPt.x;		# read geographic coordinates from 3D point to 2D point
	mapCoords.y = geogPt.y;
	mapCoords = TransPoint2D(mapCoords, transGeogUTM);		# convert 2D point to UTM
	utmPt.x = mapCoords.x;	# assign coordinate values to 3D point to return
	utmPt.y = mapCoords.y;
	utmPt.z = geogPt.z;
	return utmPt;
	}
##################################################################################
#### Convert planar coordinates in the designated UTM CRS back to geographic
func class POINT3D convertUTMtoGeog(class POINT3D utmPt) {
	local class POINT2D mapCoords;	# 2D map coordinates
	local class POINT3D geogPt;			# 3D point to return
	mapCoords.x = utmPt.x;
	mapCoords.y = utmPt.y;
	mapCoords = TransPoint2D(mapCoords, transGeogUTM, 1);	# inverse transformation
	geogPt.x = mapCoords.x;
	geogPt.y = mapCoords.y;
	geogPt.z = utmPt.z;
	return geogPt;
	}
##################################################################################
#### Find angle between grid north and true north for station location (center point)
func findAngleToNorth (class POINT3D ctr) {
	local class POINT2D mapCoords;
	local class SR_COORDREFSYS ptCRS;
	if (transGeogCoords == 1) {
		ctr = convertGeogToUTM(ctrPt);
		ptCRS = utmCRS;
		}
	else {
		ptCRS = demCRS;
		}
	return ptCRS.ComputeAngleToNorth(ctr);
	}
##################################################################################
#### Compute linear coefficients of the equation of the plane from the three
#### 3D points in the polyline tool
func class POINT3D computePlanarCoefficients( ) {
	local class POINT3D pt1, pt2, pt3;	# three points on a plane
	local class POINT3D pt21, pt31;		# differences between point positions
	local class POINT3D ctrPt;
	pt1 = setXY(1);	# get screen coordinates for polyline vertex
	pt1 = setZ(pt1);	# convert x,y to map coordinates and set Z value from DEM
	pt2 = setXY(2);
	pt2 = setZ(pt2);
	pt3 = setXY(3);
	pt3 = setZ(pt3);
	# if one of the z values is bad, then stop - return invalid ctrPt
	if (!(isZValid(pt1.z) && isZValid(pt2.z) && isZValid(pt3.z)))
	{
		string s$ = "Error: One (or more) of the points has an invalid elevation.";
		PopupMessage(s$);
		ctrPt.x = NullValue(DEM);
		ctrPt.y = NullValue(DEM);
		ctrPt.z = NullValue(DEM);
		return ctrPt;
	}
	else
	{
		# compute the center point of the triangle, this is the point that is created
		# and attributed with the information about strike and dip.
		ctrPt = computeCenterPoint(pt1, pt2, pt3);
		# if DEM map coordinates are geographic, call function to convert each vertex point
		# to planar UTM coordinates to compute equation of the plane
		if (transGeogCoords == 1) {
			pt1 = convertGeogToUTM(pt1);
			pt2 = convertGeogToUTM(pt2);
			pt3 = convertGeogToUTM(pt3);
			}
		# Compute linear coefficients for the equation of the plane
		pt21 = computeDifference(pt2, pt1);
		pt31 = computeDifference(pt3, pt1);
		a = Determinant( pt21.y, pt21.z, pt31.y, pt31.z );
		b = Determinant( pt21.x, pt21.z, pt31.x, pt31.z ) * -1;
		c = Determinant( pt21.x, pt21.y, pt31.x, pt31.y );
		# call function to compute outcrop trace
		createPlaneRaster(pt1, pt2, pt3);
		return ctrPt;
	}
}
##################################################################################
#### Use the computed planar coefficients to determine the outcrop trace pattern
#### based on the plane intersecting the DEM raster
proc createPlaneRaster(class POINT3D testPoint, class POINT3D point2, class POINT3D point3) {
	local numeric i, j, k, l, tempVal;
	local class POINT2D tempPoint;
	local numeric numLins, numCols;
	local class RVC_RASTER RastPlane;
	local class RVC_GEOREFERENCE georefPlane, georefClip;
	local class POINT3D ulPt, urPt, lrPt, llPt;
	local class CTRLPOINT ctrlPointList[4];		# array of control points
	local class RVC_VECTOR tempVect, ClipVect;
	# compute constant d (z-intercept) of plane using test point
	d = (-1) * ((a*testPoint.x) + (b* testPoint.y) + (c*testPoint.z));
	# Determine the extents of the polytool
	local class RECT polyRect = poly.ComputeExtents();
	local class TRANS2D_MAPGEN transScrnToLayer = View.GetTransLayerToScreen(rasterLayer, true);
	ulPt = polyRect.pt1;		# screen coordinates for upper left and lower right corners of polytool extents
	lrPt = polyRect.pt2;
	ulPt = transScrnToLayer.ConvertPoint2DFwd(ulPt);	# extents in DEM object coordinates
	lrPt = transScrnToLayer.ConvertPoint2DFwd(lrPt);
	# Use the user-defined variable to expand the extents of the outcrop trace
	numLins = (lrPt.y - ulPt.y) * stretchClipBoundry;
	numCols = (lrPt.x - ulPt.x) * stretchClipBoundry;
	ulPt.x = floor(ulPt.x - numCols);
	ulPt.y = floor(ulPt.y - numLins);
	lrPt.x = floor(lrPt.x + numCols);
	lrPt.y = floor(lrPt.y + numLins);
	# Make sure extents are within the bounds of the DEM raster
	if (ulPt.x < 1) then ulPt.x = 1;
	if (ulPt.y < 1) then ulPt.y = 1;
	if (lrPt.x > DEM.$Info.NumCols) then lrPt.x = DEM.$Info.NumCols;
	if (lrPt.y > DEM.$Info.NumLins) then lrPt.y = DEM.$Info.NumLins;
	# Get number of DEM lines and columns in the outcrop trace area
	numLins = lrPt.y - ulPt.y + 1;
	numCols = lrPt.x - ulPt.x + 1;
	# set raster object coordinates for upper right and lower left corners of trace area
	urPt.x = lrPt.x;	urPt.y = ulPt.y;
	llPt.x = ulPt.x;	llPt.y = lrPt.y;
	# Create temporary binary raster for outcrop trace area with georefence using same CRS as DEM
	CreateTempRaster(RastPlane, numLins, numCols, "binary");
	georefPlane.SetCoordRefSys(demCRS);
	georefPlane.SetCalibModel(demCalibModel);
	# create georeference control points for plane raster
	ctrlPointList[1].source.x = 1;	# upper left control point
	ctrlPointList[1].source.y = 1;
	ctrlPointList[1].dest = transDEMobjToMap.ConvertPoint2DFwd(ulPt);
	ctrlPointList[2].source.x = numCols;	# upper right control point
	ctrlPointList[2].source.y = 1;
	ctrlPointList[2].dest = transDEMobjToMap.ConvertPoint2DFwd(urPt);
	ctrlPointList[3].source.x = numCols;	# lower right control point
	ctrlPointList[3].source.y = numLins;
	ctrlPointList[3].dest = transDEMobjToMap.ConvertPoint2DFwd(lrPt);
	ctrlPointList[4].source.x = 1;			# lower left control point
	ctrlPointList[4].source.y = numLins;
	ctrlPointList[4].dest = transDEMobjToMap.ConvertPoint2DFwd(llPt);
	# set control point list for georeference and assign georeference to RastPlane
	georefPlane.SetPointList(ctrlPointList);
	georefPlane.Make(RastPlane);
	# Set up a status message to inform user of progress
	local class STATUSDIALOG statusD;
	statusD.Create();
	local class STATUSCONTEXT statusC = statusD.CreateContext();
	StatusSetMessage(statusC, "Computing Strike and Dip Values and Outcrop Pattern");
	StatusSetBar(statusC, i, DEM.$Info.NumLins);
	# Loop through cells in RastPlane; get Z value for corresponding cell in DEM and
	# use its map coordinates to compute Z value from the equation for the strike/dip plane.
	# If Plane is above ground, set cell value in RastPlane to 1.  After all cells are
	# processed, boundary between 1's and 0's is the outcrop trace for the plane.
	for i = 1 to RastPlane.$Info.NumLins {
		StatusSetBar(statusC, i, DEM.$Info.NumLins);
		for j = 1 to RastPlane.$Info.NumCols {
			k = ulPt.y + i - 1; 		l = ulPt.x + j - 1;		# matching object coordinates in DEM
			tempPoint.x = l;
			tempPoint.y = k;
			tempPoint = transDEMobjToMap.ConvertPoint2DFwd(tempPoint); # cell location in DEM map coordinates
			if (transGeogCoords == 1) {	# if DEM has geographic CRS, convert map coords to UTM to match plane coordinates
				tempPoint = TransPoint2D(tempPoint, transGeogUTM);
				}
			tempVal = (-1) * (a * tempPoint.x + b * tempPoint.y + d) / c;
			if (tempVal > DEM[k,l] * demZscale * unitZscale) then
				RastPlane[i, j] = 1;
			}
		}
	# Create a Vector containing the outcrop trace line.  This is created by finding the boundaries
	# between the 0's and 1's in the raster that compared the DEM to the computed plane.
	CreateTempVector(tempVect);
	tempVect = RasterToVectorBound(RastPlane, "DoImpliedGeoref");
	# Create clip area to extract outcrop trace line and leave border behind;
	# set corners of clip area in map coordinates, equivalent to 1 DEM cell smaller than trace area
	local class POINT2D rastScale;
	DEM.ComputeScaleFromGeoref(rvcGeorefDEM, rastScale, 0);	# x and y map dimensions of a DEM cell
	local class POINT2D vertexPt;
	local class POLYLINE polyL;
	vertexPt.x = ctrlPointList[1].dest.x;	# upper left vertex
	vertexPt.y = ctrlPointList[1].dest.y;
	polyL.AppendVertex(vertexPt);
	vertexPt.x = ctrlPointList[2].dest.x - rastScale.x;	# upper right vertex
	vertexPt.y = ctrlPointList[2].dest.y;
	polyL.AppendVertex(vertexPt);
	vertexPt.x = ctrlPointList[3].dest.x - rastScale.x;	# lower right vertex
	vertexPt.y = ctrlPointList[3].dest.y + rastScale.y;
	polyL.AppendVertex(vertexPt);
	vertexPt.x = ctrlPointList[4].dest.x;	# lower left vertex
	vertexPt.y = ctrlPointList[4].dest.y + rastScale.y;
	polyL.AppendVertex(vertexPt);
	polyL.SetClosed(1);		# close the polyline
	polyRect = polyL.ComputeExtents();
	CreateTempVector(ClipVect, "VectorTookit,Polygonal","NoDatabase", polyRect);
	local class RVC_GEOREFERENCE georefClipVect;
	georefClipVect.SetCoordRefSys(demCRS);
	georefClipVect.SetImplied();
	georefClipVect.Make(ClipVect);
	VectorAddPolyLine(ClipVect, polyL);
	VectorValidate(ClipVect);
	# Extract outcrop trace line to new vector using the ClipVect to leave behind the border created
	# by the autobounds procedure.
	CreateTempVector(intersectV);
	intersectV = VectorExtract(ClipVect, tempVect, "InsideClip", "", "1", "1", "0");
	# if DEM has geographic coordinates, reproject UTM points passed to this function back to geographic
	# so they can be added to the outcrop trace vector
	if (transGeogCoords == 1) {
		testPoint = convertUTMtoGeog(testPoint);
		point2 = convertUTMtoGeog(point2);
		point3 = convertUTMtoGeog(point3);
		}
	# Add the three points selected by the user
	VectorAddPoint(intersectV, testPoint.x, testPoint.y, testPoint.z);
	VectorAddPoint(intersectV, point2.x, point2.y, point2.z);
	VectorAddPoint(intersectV, point3.x, point3.y, point3.z);
	# Create a Point Database table and record with the station number and color values
	iVpointdb = OpenVectorPointDatabase(intersectV);
	iVpointtable = TableCreate(iVpointdb, "Outcrop", "Intersecting Outcrop Pattern");
	TableAddFieldInteger(iVpointtable, "Station", 4);
	TableAddFieldInteger(iVpointtable, "RedComp", 3);
	TableAddFieldInteger(iVpointtable, "GreenComp", 3);
	TableAddFieldInteger(iVpointtable, "BlueComp", 3);
	local numeric recNumPoints = TableWriteRecord(iVpointtable, 0, -999, 0, 0, 0);
	# Create a Line Database table and with the station number and color values
	iVlinedb = OpenVectorLineDatabase(intersectV);
	iVlinetable = TableCreate(iVlinedb, "Outcrop", "Intersecting Outcrop Pattern");
	TableAddFieldInteger(iVlinetable, "Station", 4);
	TableAddFieldInteger(iVlinetable, "RedComp", 3);
	TableAddFieldInteger(iVlinetable, "GreenComp", 3);
	TableAddFieldInteger(iVlinetable, "BlueComp", 3);
	local numeric recNumLines = TableWriteRecord(iVlinetable, 0, -999, 0, 0, 0);
	# Attach all the points to the station
	local numeric numPoints = NumVectorPoints(intersectV);
	local numeric k;
	local array numeric rec[1];
	rec[1] = recNumPoints;
	for k = 1 to numPoints {
		TableWriteAttachment(iVpointtable, k, rec, 1, "point");
	}
	# Attach all the lines to the station
	local numeric numLines = NumVectorLines(intersectV);
	rec[1] = recNumLines;
	for k = 1 to numLines {
		TableWriteAttachment(iVlinetable, k, rec, 1, "line");
	}
	# Add the outcrop trace line to the view for the user to evaluate
	intersectVlayer = GroupQuickAddVectorVar(firstGroup, intersectV);
	intersectVlayer.Point.NormalStyle.DrawCircle = 1;
	intersectVlayer.Point.NormalStyle.Fill = 1;
	intersectVlayer.Line.StyleMode = "AllSame";
	intersectVlayer.Line.NormalStyle.Color = btnLINECOLOR.GetColor();
	intersectVlayer.Line.NormalStyle.Width = 1.05;
	ViewRedrawLayer(View, intersectVlayer);
	# Clean up the resources used
	statusC.Clear();
	statusD.Destroy();
	CloseVector(tempVect);
	CloseVector(ClipVect);
	CloseRaster(RastPlane);
#	CloseRaster(clipRast);
#	CloseRaster(PLANE);
}
##################################################################################
#### Function to return the determinant of a 2 x 2 matrix formed by the four
#### numbers in the argument.  The numbers can be in either row or column order;
#### the result is the same in either order.
func Determinant(a1, b1, a2, b2) {
	local det;
	det =  a1 * b2 - a2 * b1 ;
	return det;
}
##################################################################################
#### Check that the z value is not null (either null or beyond rast extents)
func isZValid(numeric z) {
	local numeric nullval = NullValue(DEM);
	return z != nullval;
}
##################################################################################
#### Compute the strike and dip for a vertical plane
proc computeVerticalPlane( ) {
	local numeric northAng;
	northAng = findAngleToNorth(ctrPt);
	dipAng = 90;
	if ( a == 0 )
	{
		# If a = b = c = 0, points are colinear, can't find plane;
		# tool script should trap this as an error and alert the user
		# to repick points -- tool shouldn't let this happen -djg
		if ( b == 0 )
		{
			strikeAng = null;
			dipDir = null;
			dipAng = null;
			PopupMessage("Points are colinear and so don't define a unique plane.");
		}
		# Plane is also parallel to x-axis, raw strike = 90 and arctan(-b/a) is undefined
		else
		{
			strikeAng = 90 + northAng;
		}
	}
	else
	{
		# Raw strike angle = arctan of the partial derivative of x versus y = (-b / a)
		# (this is slope of line of intersection of plane with horizontal = strike line
		# in plot of x vs y, which gives the angle equivalent to azimuth in conventional
		# y vs x plot)
		strikeAng = atand( -b / a ) + northAng;
		if ( strikeAng < 0 ) then
			strikeAng = 180 + strikeAng;
 	}
}
##################################################################################
#### Compute the strike and dip for a nonvertical plane
proc computeNonVerticalPlane( ) {
	local numeric northAng;
	northAng = findAngleToNorth(ctrPt);
	# Partial derivative of dip in x-direction (opposite the slope component)
	dipX = a / c;
	# Plane is parallel to y-axis (north-south), dipY is 0;
	if ( b == 0 ) {
		# Plane is parallel to both x and y and thus is horizontal
		if ( a == 0 ) {
			dipAng = 0;
			strikeAng = null;	# strikeAng is undefined
			dipDir = null;
		# Plane is not horizontal, raw strike is north-south
		} else {
			dipAng = atand( dipX );
			if ( dipAng < 0 ) {
				dipAng = -1 * dipAng;
			}
			# Dip direction is due east
			if ( dipX > 0 ) {
				dipDir = 90 + northAng;
				strikeAng = 0 + northAng;
			} else {
				dipDir = 270 + northAng;
				strikeAng = 180 + northAng;
			}
		}
	# Strike not parallel to y-axis, dipY is nonzero
	} else {
		# Partial derivative of dip in y-direction (opposite the slope component)
		dipY = b / c;
		# Compute dip angle from dipX and dipY
		dipAng = atand( sqrt( sqr(dipX) + sqr(dipY) ) );
		if ( dipAng < 0 ) {
			dipAng = -1 * dipAng;
		}
		# Compute azimuth of dip direction from dipX and dipY; value returned by atand() is between -90 and +90,
		# so must adjust result to get positive azimuth value between 0 and 360 depending on which quadrant dip
		# vector is in.
		dipDir = atand( dipX / dipY ) + northAng;
		# Dip direction is due north or south
		if ( dipDir == 0 ) {
			# Dip direction is south
			if ( dipY < 0 ) {
				dipDir = 180;
			}
		}
		# Dip direction in SE or NW quadrants, must adjust value
		if ( dipDir < 0 ) {
			# Dip direction in SE quadrant
			if ( dipX > 0 ) {
				dipDir = 180 + dipDir;
			# Dip direction in NW quadrant
			} else {
				dipDir = 360 + dipDir;
			}
		# Dip direction in SW quadrant, must adjust value
		} else if ( dipX < 0 ) {
			dipDir = 180 + dipDir;
		}
		# Compute strike angle from dip direction and convert negative values to
		# positive azimuth in range 0 to 360
		# Strike direction from right-hand rule
		strikeAng = dipDir - 90;
		if ( strikeAng < 0 ) {
			strikeAng = 360 + strikeAng;
		}
	}
}
##################################################################################
#### Compute strike and dip values based on the user input
proc computeStrikeDip( ) {
	local numeric isVertical = (c==0);
	if (isVertical) {
		computeVerticalPlane();
	} else {
		computeNonVerticalPlane();
	}
}
##################################################################################
#### Add the point to the StrikeDipVector and return its element number
func addPointToVector(class POINT3D pt) {
	local numeric x = pt.x;
	local numeric y = pt.y;
	local numeric z = pt.z;
	VectorAddPoint(StrikeDipVector, x, y, z);
	VectorValidate(StrikeDipVector);
	return FindClosestPoint(StrikeDipVector, x, y, SDgeoref);
}
##################################################################################
#### Add a record to the StrikeDipVector point database with appropriate values
proc addRecordToDatabase( numeric pointElemNumAdded ) {
	local class POINT3D pt1, pt2, pt3;
	# Get the original three points from the polytool
	pt1 = setXY(1);
	pt1 = setZ(pt1);
	pt2 = setXY(2);
	pt2 = setZ(pt2);
	pt3 = setXY(3);
	pt3 = setZ(pt3);
	# If editing, use the same station number
	if (Current_Mode == "EDIT") {
		station = deletedStation;
	# If adding a new one, set as the highest station incremented by 1
	} else {
		station = 0;
		local numeric tempVal;
		local numeric numRecords = NumRecords(SDtable);
		local numeric i;
		for i = 1 to numRecords {
			tempVal = TableReadFieldNum(SDtable, "Station", i);
			if (tempVal > station) {
				station = tempVal;
			}
		}
		station = station + 1;
	}
	# Round and save Strike Dip Values as Integers
	strikeAng = round(strikeAng);
	dipAng = round(dipAng);
	dipDir = round(dipDir);
	# Write the record for the Strike Dip point and attach it
	records[1] = TableWriteRecord(SDtable, 0, station, strikeAng, dipAng, dipDir, overturned,
	                                            pt1.x, pt1.y, pt1.z, pt2.x, pt2.y, pt2.z, pt3.x, pt3.y, pt3.z, 0);
	TableWriteAttachment(SDtable, pointElemNumAdded, records, 1, "point");
	# Write records into the temporary outcrop trace pattern vector.  They will be transfered
	# to the actual outcrop trace vector in the next step.
	local class COLOR tempColor = btnLINECOLOR.GetColor();
	TableWriteRecord(iVpointtable, 1, station, tempColor.red, tempColor.green, tempColor.blue);
	TableWriteRecord(iVlinetable, 1, station, tempColor.red, tempColor.green, tempColor.blue);
	# Vector Copy Elements doesn't like copying to an object that already has elements, so if
	# it does, we introduce a temporary vector before using the VectMerge process.
	if ( NumVectorPoints(OutcropPVector) + NumVectorLines(OutcropPVector) > 0 ) {
		# Close the databases and remove the layers
		CloseDatabase(SDdb);
		CloseDatabase(OPpointdb);
		CloseDatabase(OPlinedb);
		GroupRemoveLayer(firstGroup, SDvectorLayer);
		GroupRemoveLayer(firstGroup, OPvectorLayer);
		# Copy the temporary outcrop trace vector to a temporary vector, then merge it with the
		# final destination outcrop trace pattern vector
		local VECTOR tempMergeVector;
		CreateTempVector(tempMergeVector);
		VectorCopyElements(OutcropPVector, tempMergeVector);
		OutcropPVector = VectMerge(tempMergeVector, intersectV);
		# Add the layers back to the display;
		addVectorsToDisplay();
		# Reinitialize and clean up the vector databases
		initializeBeddingDatabaseTable();
		initializeOutcropDatabaseTable();
		TableRemoveDuplicateRecords(OPpointtable);
		TableRemoveDuplicateRecords(OPlinetable);
		CloseVector(tempMergeVector);
	# If the outcrop trace vector has no elements, just copy straight to it.
	} else {
		# Close the databases and remove the layers
		CloseDatabase(SDdb);
		CloseDatabase(OPpointdb);
		CloseDatabase(OPlinedb);
		GroupRemoveLayer(firstGroup, SDvectorLayer);
		GroupRemoveLayer(firstGroup, OPvectorLayer);
		VectorCopyElements(intersectV, OutcropPVector);
		# Add the layers back to the display;
		addVectorsToDisplay();
		# Reinitialize and clean up the vector databases
		initializeBeddingDatabaseTable();
		initializeOutcropDatabaseTable();
		TableRemoveDuplicateRecords(OPpointtable);
		TableRemoveDuplicateRecords(OPlinetable);
	}
}
##################################################################################
#### This function creates the Query used to display the appropriate lines attached
#### to specific stations.  It also manages the color they appear in.
proc UpdateQuery( ) {
	currentQuery = genericQuery + "-1)";
	# Scan through every record in the StrikeDipVector Point database
	local numeric numPoints = NumVectorPoints(StrikeDipVector);
	local array numeric records[10];
	local numeric i;
	for i = 1 to numPoints {
		TableReadAttachment(SDtable, i, records, "point");
		# If the record value states that it should be displayed, add an extra OR statement
		# to the current query so it is displayed.
		if (TableReadFieldNum(SDtable, "DisplayOPdata", records[1]) != 0) {
			currentQuery = currentQuery + " || " + genericQuery
	                       + NumToStr(TableReadFieldNum(SDtable, "Station", records[1])) + ")";
		}
	}
	# Set the current query and redraw the layer
	OPpoints.Select.Query = currentQuery;
	OPlines.Select.Query = currentQuery;
	ViewRedrawLayer(View, OPvectorLayer);
}
##################################################################################
#### This function is designed to delete all lines, points and database records that
#### are related to the station that is being deleted.
proc deleteStation( ) {
	local numeric numToDelete = 0;
	local numeric i;
	local array numeric tempRecords[10];
	local numeric tempStation;
	local numeric numRecords;
	# Turn off the highlighted point
	SDpoints.HighlightSingle(selectedElementNum, "Toggle");
	# Get the station number from the selected point
	deletedStation = StrikeDipVector.point[selectedElementNum].Bedding.Station;
	# Scan through and mark all lines attached to the deleted station to be deleted.
	local numeric numLines = NumVectorLines(OutcropPVector);
	local array numeric linesToDelete[numLines];
	for i = 1 to numLines {
		tempStation = OutcropPVector.line[i].Outcrop.Station;
		if (tempStation == deletedStation) {
			# Unattach the records if the line is being marked to delete
			numRecords = TableReadAttachment(OPlinetable, i, tempRecords, "line");
			TableRemoveAttachment(OPlinetable, i, tempRecords, numRecords, "line");
			# Save line number to delete
			numToDelete = numToDelete + 1;
			linesToDelete[numToDelete] = i;
		}
	}
	# Delete all the marked lines
	if (numToDelete > 0) {
		VectorDeleteLines(OutcropPVector, linesToDelete, numToDelete);
	}
	numToDelete = 0;
	# Scan through and mark all points attached to the deleted station to be deleted.
	local numeric numPoints = NumVectorPoints(OutcropPVector);
	local array numeric pointsToDelete[numPoints];
	for i = 1 to numPoints {
		tempStation = OutcropPVector.point[i].Outcrop.Station;
		if (tempStation == deletedStation) {
			# Unattach the records if the point is being marked to delete
			numRecords = TableReadAttachment(OPpointtable, i, tempRecords, "point");
			TableRemoveAttachment(OPpointtable, i, tempRecords, numRecords, "point");
			# Save point number to delete
			numToDelete = numToDelete + 1;
			pointsToDelete[numToDelete] = i;
		}
	}
	# Delete all the marked points
	if (numToDelete > 0) {
		VectorDeletePoints(OutcropPVector, pointsToDelete, numToDelete);
	}
	numToDelete = 0;
	# Scan through and delete the Strike and Dip station
	numPoints = NumVectorPoints(StrikeDipVector);
	local numeric pointToDelete;
	for i = 1 to numPoints {
		tempStation = StrikeDipVector.point[i].Bedding.Station;
		if (tempStation == deletedStation) {
			# Unnattach the records if the point is being marked to delete
			numRecords = TableReadAttachment(SDtable, i, tempRecords, "point");
			TableRemoveAttachment(SDtable, i, tempRecords, numRecords, "point");
			# Save point number to delete
			numToDelete = numToDelete + 1;
			pointsToDelete[numToDelete] = i;
		}
	}
	# Delete all the marked points
	if (numToDelete > 0) {
		VectorDeletePoints(StrikeDipVector, pointsToDelete, numToDelete);
	}
	# Remove all the records that are no longer attached
	TableRemoveUnattachedRecords(OPlinetable);
	TableRemoveUnattachedRecords(OPpointtable);
	TableRemoveUnattachedRecords(SDtable);
	VectorValidate(OutcropPVector);
	VectorValidate(StrikeDipVector);
}
##################################################################################
#### Creates the vector objects to store the resultant points and outcrop trace lines
proc createDestVectors( ) {
	local string flags$ = "Polygonal,3DVector";
	GetOutputVector(StrikeDipVector, flags$, "", Group.Extents);
	GetOutputVector(OutcropPVector, flags$, "", Group.Extents);
	createVectorGeoreference();
}
##################################################################################
#### Creates the vector georeference from the raster if necessary
func createVectorGeoreference( ) {
	local class Georef rgeoref;
	# determine if any georef objects exist
	local string file$ = StrikeDipVector.$INFO.Filename;
	local string type$ = "GEOREF";
	local numeric parent = GetObjectNumber(StrikeDipVector);
	local class STRINGLIST sl = GetAllObjectNames(file$, type$, parent);
	local numeric hasGeoref = (sl.GetNumItems()!=0);
	# If it has a Georef, use it.  If it doesn't, create one from the DEM
	if (hasGeoref)	{
		SDgeoref = GetLastUsedGeorefObject(StrikeDipVector);
	} else {
		rgeoref = GetLastUsedGeorefObject(DEM);
		CreateImpliedGeoref(StrikeDipVector, rgeoref.CoordRefSys);
		SDgeoref = GetLastUsedGeorefObject(StrikeDipVector);
		CreateImpliedGeoref(OutcropPVector, rgeoref.CoordRefSys);
	}
}
##################################################################################
#### This function is used to check if the vector is already in a layer in the View
func vectorLayerExists( string name ) {
	local class VECTOR V;
	# Set the active layer to the top layer in the first group
	firstGroup.SetActiveLayer(firstGroup.LastLayer);
	# Cycle through all the layers in the group
	while (firstGroup.ActiveLayer != 0) {
		if (firstGroup.ActiveLayer.TypeID == "Vector") {
			DispGetVectorFromLayer(V, firstGroup.ActiveLayer);
			# If the vector layer has the same name, return true
			if (V.$Info.Name == name) {
				CloseVector(V);
				return true;
			}
			CloseVector(V);
		}
		# If there are no more layers to check, return false
		if (firstGroup.ActiveLayer.PrevLayer == 0) {
			return false;
		}
		# Set the active layer to the next one down
		firstGroup.SetActiveLayer(firstGroup.ActiveLayer.PrevLayer);
	}
	return false;
}
##################################################################################
#### Add the vector to the display group if it is not already there
proc addVectorsToDisplay( ) {
	# Check if the StrikeDipVector layer does not already exist
	if(!vectorLayerExists(StrikeDipVector.$Info.Name))
	{
		# Add the vector to the View
		SDvectorLayer = GroupQuickAddVectorVar(firstGroup, StrikeDipVector);
		SDvectorLayer.IgnoreExtents = 1;
	} else {
		SDvectorLayer = firstGroup.ActiveLayer;
		CloseVector(StrikeDipVector);
		DispGetVectorFromLayer(StrikeDipVector, SDvectorLayer);
		SDvectorLayer.IgnoreExtents = 1;
	}
	# Check if the OutcropPVector layer does not already exist
	if(!vectorLayerExists(OutcropPVector.$Info.Name))
	{
		# Add the vector to the View
		OPvectorLayer = GroupQuickAddVectorVar(firstGroup, OutcropPVector);
		OPvectorLayer.IgnoreExtents = 1;
	} else {
		OPvectorLayer = firstGroup.ActiveLayer;
		CloseVector(OutcropPVector);
		DispGetVectorFromLayer(OutcropPVector, OPvectorLayer);
		OPvectorLayer.IgnoreExtents = 1;
	}
	# Set up how the layers display points, lines and polygons
	SDpoints = SDvectorLayer.Point;
	OPpoints = OPvectorLayer.Point;
	OPlines = OPvectorLayer.Line;
	OPpolys = OPvectorLayer.Poly;
	# Use the modified bedding.qry CartoScript developed by MicroImages to
	# be used by this ToolScript if it exists in the same directory
	SDpoints.Select.Mode = "All";
	SDvectorLayer.Point.StyleMode = "ByScript";
	SDvectorLayer.Point.Script = sdQry$;
	SDvectorLayer.SaveDispParmSubObject();
	# Create a script that will style the points by color
	local string script = "DrawColor$ = NumToStr(Outcrop.RedComp) + \" \" + NumToStr(Outcrop.GreenComp)"
	                          + " + \" \" + NumToStr(Outcrop.BlueComp);\n";
	# Set up the points be selected by query and styled by script
	OPpoints.Select.Mode = "ByQuery";
	OPpoints.Select.Query = "Outcrop.Station == -1";
	OPpoints.StyleMode = "ByScript";
	OPpoints.Script = script;
	# Add to the script to make the lines thicker
	script = script + "LineScale = 1;";
	# Set up the lines be selected by query and styled by script
	OPlines.Select.Mode = "ByQuery";
	OPlines.Select.Query = "Outcrop.Station == -1";
	OPlines.StyleMode = "ByScript";
	OPlines.Script = script;
	# Turn off polygons
	OPpolys.Select.Mode = "None";
}
##################################################################################
#### Create the vector point database tables for the StrikeDipVector
proc initializeBeddingDatabaseTable( ) {
	local string name$ = "Bedding";
	local string desc$ = sprintf("Strike and dip of bedding; map coordinates in %s, elevation in %s", coordUnitName$, zCoordUnit$);
	SDdb = OpenVectorPointDatabase(StrikeDipVector);
	# If the table doesn't exist it needs to be created
	if( !TableExists(SDdb, name$) ) {
		SDtable = TableCreate(SDdb, name$, desc$);
		SDtable.OneRecordPerElement = 1;
		# Create the fields
		createBeddingFields();
	# Else the table does exist, get the info for it
	} else {
		SDtable = DatabaseGetTableInfo(SDdb, name$);
	}
	SDstationfield = FieldGetInfoByName(SDtable, "Station");
}
##################################################################################
#### Create the fields in the bedding table
proc createBeddingFields( ) {
	TableAddFieldInteger(SDtable, "Station", 8);
	TableAddFieldInteger(SDtable, "StrikeAngle", 11);
	TableAddFieldInteger(SDtable, "DipAngle", 8);
	TableAddFieldInteger(SDtable, "DipDirection", 12);
	TableAddFieldInteger(SDtable, "Overturned", 10);
	TableAddFieldFloat(SDtable, "ToolPos1_x", 12, 4);
	TableAddFieldFloat(SDtable, "ToolPos1_y", 12, 4);
	TableAddFieldFloat(SDtable, "ToolPos1_z", 10, 1);
	TableAddFieldFloat(SDtable, "ToolPos2_x", 12, 4);
	TableAddFieldFloat(SDtable, "ToolPos2_y", 12, 4);
	TableAddFieldFloat(SDtable, "ToolPos2_z", 10, 1);
	TableAddFieldFloat(SDtable, "ToolPos3_x", 12, 4);
	TableAddFieldFloat(SDtable, "ToolPos3_y", 12, 4);
	TableAddFieldFloat(SDtable, "ToolPos3_z", 10, 1);
	TableAddFieldInteger(SDtable, "DisplayOPdata", 13);
}
##################################################################################
#### Create the vector point database tables for the OutcropPVector
proc initializeOutcropDatabaseTable( ) {
	local string name$ = "Outcrop";
	local string desc$ = "Intersecting Outcrop Pattern";
	OPpointdb = OpenVectorPointDatabase(OutcropPVector);
	OPlinedb = OpenVectorLineDatabase(OutcropPVector);
	# If the table doesn't exist it needs to be created
	if( !TableExists(OPpointdb, name$) ) {
		OPpointtable = TableCreate(OPpointdb, name$, desc$);
		OPpointtable.OneRecordPerElement = 1;
		# Add the fields in the table
		TableAddFieldInteger(OPpointtable, "Station", 4);
		TableAddFieldInteger(OPpointtable, "RedComp", 3);
		TableAddFieldInteger(OPpointtable, "GreenComp", 3);
		TableAddFieldInteger(OPpointtable, "BlueComp", 3);
	# Else the table does exist, get the info for it
	} else {
		OPpointtable = DatabaseGetTableInfo(OPpointdb, name$);
	}
	# If the table doesn't exist it needs to be created
	if( !TableExists(OPlinedb, name$) ) {
		OPlinetable = TableCreate(OPlinedb, name$, desc$);
		OPlinetable.OneRecordPerElement = 1;
		# Add the fields in the table
		TableAddFieldInteger(OPlinetable, "Station", 4);
		TableAddFieldInteger(OPlinetable, "RedComp", 3);
		TableAddFieldInteger(OPlinetable, "GreenComp", 3);
		TableAddFieldInteger(OPlinetable, "BlueComp", 3);
	# else the table does exist, get info for it
	} else {
		OPlinetable = DatabaseGetTableInfo(OPlinedb, name$);
	}
}
##################################################################################
#### Checks the first layer to see if it is valid DEM raster
func checkLayer( ) {
	local numeric valid = false;
	local class GRE_LAYER layer = firstGroup.FirstLayer;
	if (layer.TypeID == "Surface") layer = layer.NextLayer;
	# Get name of active layer if it is usable.  If not output an error message.
	if (layer.TypeID != "")
	{
		if (layer.TypeID == "Raster")
		{
			DispGetRasterFromLayer(DEM, layer);
			if (DEM.$Info.Type == "binary" or
				DEM.$Info.Type == "4-bit unsigned" or
				DEM.$Info.Type == "8-bit signed" or
				DEM.$Info.Type == "8-bit unsigned" or
				DEM.$Info.Type == "16-bit signed" or
				DEM.$Info.Type == "16-bit unsigned" or
				DEM.$Info.Type == "32-bit signed" or
				DEM.$Info.Type == "32-bit unsigned" or
				DEM.$Info.Type == "32-bit floating-point" or
				DEM.$Info.Type == "64-bit signed" or
				DEM.$Info.Type == "64-bit unsigned" or
				DEM.$Info.Type == "64-bit floating-point")
			{
				rasterLayer = layer;
				valid = true;
			}
		}
	}
	return valid;
}
###########################################################
### Procedure to check georeference of DEM; if it uses a geographic coordinate
### system, set up a coordinate transformation to the appropriate UTM zone
### to transform map coordinate used in determine strike/dip values.
### Called by onInitialize() and by cbLayer()
proc checkGeoref () {
	local class RECT3D extentsObj;
	local class POINT2D center;
	local class STRING demCRSunit$;
	transGeogCoords = 0;
	unitZscale = 1;		# default value for CRS coordinates in meters; assumes elevations in meters
	demZscale = DEM.$Info.DataScale;
	rvcGeorefDEM.OpenLastUsed(DEM);
	demCRS = rvcGeorefDEM.GetCoordRefSys();
	demCalibModel = rvcGeorefDEM.GetCalibModel();
	rvcGeorefDEM.GetTransParm(transDEMobjToMap, 0, demCalibModel);
	demCRSunit$ = demCRS.CoordSys.GetAxis(1).Unit.GetSymbol();
	coordUnitName$ = demCRS.CoordSys.GetAxis(1).Unit.GetName(1);
	zCoordUnit$ = "meters";
	if (demCRS.IsProjected() <> 1) {		# DEM CRS is not projected
		if (demCRS.IsLocal() <> 1) {		# DEM CRS also is not local, then must be geographic (lat/lon)
			transGeogCoords = 1;
			DEM.GetExtents(extentsObj);
			center = extentsObj.center;
			center = transDEMobjToMap.ConvertPoint2DFwd(center);  # DEM center point in map coordinates, used to set UTM zone
			# variables for setting up UTM CRS for transforming DEM map coordinates in tool
			local class SR_COORDSYS utmCoordSys;		# coordinate system
			local class SR_COORDOPDEF projectDef;		# projection parameters and method
			local class SR_COORDOPMETHOD method;
			local class SR_DATUM datum;
			datum = demCRS.Datum;		# get datum from the DEM raster
			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
			local 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);
			# setup coordinate transformation from DEM's geographic CRS to the defined UTM CRS.
			transGeogUTM.InputCoordRefSys = demCRS;
			transGeogUTM.OutputCoordRefSys = utmCRS;
			}
		}
	else if (demCRSunit$ <> "m")
		{
		unitZscale = GetUnitConvDist("m", demCRSunit$);		# scale elevations to same units as CRS coordinate system
		zCoordUnit$ = coordUnitName$;
		}
	}
##################################################################################
#### Callback for Apply button press for polyline (right click as well).
proc cbToolApply()
{
	# Make sure the bottom layer is still the DEM
	if ( checkLayer() )	{
		promptString.SetLabel("Computing Data...");
		poly = tool.GetPolygon();
		numeric numPoints = poly.GetNumPoints();
		# If there are not exactly 3 points in accepted polyline
		if ( numPoints != 3 ) {
			string message$ = "Polyline must have exactly 3 points: "+ NumToStr(numPoints) + " found.";
			PopupMessage(message$);
		# Else there are 3 good points so do computations
		} else {
			# If the three points have been modified during the add process and are being
			# recomputed, remove the original results
			if (waitingForAddConfirmation) {
				GroupRemoveLayer(firstGroup, intersectVlayer);
				CloseVector(intersectV);
			}
			# Compute the coefficients from the 3 points
			ctrPt = computePlanarCoefficients();
			if (isZValid(ctrPt.z)) {
				computeStrikeDip();
				waitingForAddConfirmation = 1;
				# Manage the Overturned toggle button
				if (Current_Mode == "ADD" or Current_Mode == "EDIT")
					overturnedToggle.SetValue(0,0);		# set toggle to OFF without calling the change callback
				# Display the resulting Strike and Dip values in the dialog
				strikeAng = round(strikeAng);
				strikeFld.SetValue(strikeAng, 0);
				dipFld.SetValue(round(dipAng), 0);
				dipDirFld.SetValue(round(dipDir), 0);
				overturnedToggle.SetEnabled(1);
				# Manage the other controls in the display
				aLINECOLOR_LABEL.SetEnabled(1);
				btnLINECOLOR.SetEnabled(1);
				btnACCEPT.SetEnabled(1);
				btnCANCEL.SetEnabled(1);
				if (Current_Mode == "ADD") {
					rGpMODE.SetEnabled(0);
					btnCLOSE.SetEnabled(0);
					promptString.SetLabel("Accept this New point?");
				} else if (Current_Mode == "EDIT") {
					promptString.SetLabel("Accept this Edited point?");
					rGpMODE.SetEnabled(0);
					btnCLOSE.SetEnabled(0);
				}
			}
		}
	# Else the bottom layer is not a valid DEM raster
	} else {
		string message$ = "The First Layer: "+activeLayer.Name+" is invalid.";
		PopupMessage(message$);
	}
}
##################################################################################
#### Callback for when the active layer changes
proc cbLayer( ) {
	checkLayer();
	checkGeoref();
}
##################################################################################
#### Callback for when the active group changes
proc cbGroup( ) {
	activegroup = Layout.ActiveGroup;
	WidgetAddCallback(activegroup.LayerSelectedCallback, cbLayer);
	cbLayer();
}
##################################################################################
#### Close the window, switching to default tool
proc cbClose( ) {
	tool.Managed = 0;
	if (setDefaultWhenClose) {
		setDefaultWhenClose = false;
		View.SetDefaultTool();
	}
}
##################################################################################
#### Manages the control dialog when the active mode has been changed
proc OnModeChange( ) {
	# Get the current Mode
	Current_Mode = rGpMODE.GetSelected();
	# Set controls to default values
	strikeFld.SetValue(0,0);
	dipFld.SetValue(0, 0);
	dipDirFld.SetValue(0, 0);
	overturnedToggle.SetEnabled(1);
	overturnedToggle.SetValue(0,0);
	overturnedToggle.SetEnabled(0);
	aLINECOLOR_LABEL.SetEnabled(0);
	btnLINECOLOR.SetEnabled(0);
	btnACCEPT.SetEnabled(0);
	btnCANCEL.SetEnabled(0);
	btnDELETE.SetEnabled(0);
	# Based on Mode, set the controls and tool to appropriate settings
	if ( Current_Mode == "VIEW") {
		promptString.SetLabel("Select point to Toggle Data ON/OFF");
		tool.managed = 0;
		ScriptGadget.Managed = 1;
	} else if ( Current_Mode == "ADD" ) {
		aLINECOLOR_LABEL.SetEnabled(1);
		btnLINECOLOR.SetEnabled(1);
		promptString.SetLabel("Select 3 new points");
		poly.Clear();
		tool.SetPolygon(poly);
		ScriptGadget.Managed = 0;
		tool.managed = 1;
	} else if ( Current_Mode == "EDIT" ) {
		promptString.SetLabel("Select Point to Edit");
		tool.managed = 0;
		ScriptGadget.Managed = 1;
	} else if ( Current_Mode == "DELETE" ) {
		promptString.SetLabel("Select Point to Delete");
		tool.managed = 0;
		ScriptGadget.Managed = 1;
	}
}
##################################################################################
#### Called when the user clicks on the "Set as Overturned" togglebutton.  It
#### adjusts the value of the Strike angle by 180 degrees in appropriate modes.
proc OnOverturnedChanged( ) {
	# If strike and dip fields do not have default 0 values
	if ( !(strikeFld.GetValueNum() == 0 && dipFld.GetValueNum() == 0 && dipDirFld.GetValueNum() == 0) )
		{
		if (Current_Mode == "ADD" or (Current_Mode == "EDIT" && currentlyEditing) )
			{
			# If the angle is less then 180 degrees and is being set to overturned
			if ( (strikeAng < 180) )
				strikeAng = strikeAng + 180;
			# If the agnle is already overturned and being set as not overturned
			else if ( (strikeAng >= 180) )
				strikeAng = strikeAng - 180;
			strikeFld.SetValue(strikeAng,0);
			}
		}
	}
##################################################################################
#### Called when the user clicks on the Accept button in the dialog.  Actions are
#### based on  the current Mode of the ToolScript
proc OnAccept( ) {
	# If in VIEW mode, toggle the display values of the selected station and update
	if (Current_Mode == "VIEW") {
		promptString.SetLabel("Resetting Query...");
		btnACCEPT.SetEnabled(0);
		btnCANCEL.SetEnabled(0);
		# Toggle the values in the database
		local array numeric records[10];
		TableReadAttachment(SDtable, selectedElementNum, records, "point");
		if (TableReadFieldNum(SDtable, "DisplayOPdata", records[1]) == 0) {
			TableWriteField(SDtable, records[1], "DisplayOPdata", 1);
		} else {
			TableWriteField(SDtable, records[1], "DisplayOPdata", 0);
		}
		UpdateQuery();
		SDpoints.HighlightSingle(selectedElementNum, "Toggle");
		selectedElementNum = 0;
		rGpMODE.SetEnabled(1);
		btnCLOSE.SetEnabled(1);
		promptString.SetLabel("Select point to Toggle Data ON/OFF");
	# If in ADD Mode, add the computed Strike and Dip data to the Destination Vectors
	} else if (Current_Mode == "ADD") {
		promptString.SetLabel("Saving Data...");
		overturned = overturnedToggle.GetValue();
		overturnedToggle.SetEnabled(0);
		aLINECOLOR_LABEL.SetEnabled(0);
		btnLINECOLOR.SetEnabled(0);
		btnACCEPT.SetEnabled(0);
		btnCANCEL.SetEnabled(0);
		# Remove the temporary vectors and save the Strike and Dip data
		GroupRemoveLayer(firstGroup, intersectVlayer);
		addRecordToDatabase( addPointToVector(ctrPt) );
		CloseVector(intersectV);
		waitingForAddConfirmation = 0;
		# Reset the polytool
		poly.Clear();
		tool.SetPolygon(poly);
		# Update the display
		UpdateQuery();
		ViewRedrawLayer(View, SDvectorLayer);
		selectedElementNum = 0;
		# Manage the dialog controls
		aLINECOLOR_LABEL.SetEnabled(1);
		btnLINECOLOR.SetEnabled(1);
		rGpMODE.SetEnabled(1);
		btnCLOSE.SetEnabled(1);
		strikeFld.SetValue(0, 0);
		dipFld.SetValue(0, 0);
		dipDirFld.SetValue(0, 0);
		promptString.SetLabel("Select 3 new points");
	# If in EDIT Mode, reload the data or save the new edited data
	} else if (Current_Mode == "EDIT") {
		# Check if the point has already been selected and in modifying mode
		if ( currentlyEditing ) {
			promptString.SetLabel("Editing Data...");
			overturned = overturnedToggle.GetValue();
			overturnedToggle.SetEnabled(0);
			aLINECOLOR_LABEL.SetEnabled(0);
			btnLINECOLOR.SetEnabled(0);
			btnACCEPT.SetEnabled(0);
			btnCANCEL.SetEnabled(0);
			GroupRemoveLayer(firstGroup, intersectVlayer);
			# Delete the old station information and add the new station information
			deleteStation();
			addRecordToDatabase( addPointToVector(ctrPt) );
			CloseVector(intersectV);
			waitingForAddConfirmation = 0;
			# Update the display
			UpdateQuery();
			ViewRedrawLayer(View, SDvectorLayer);
			# Reset the polytool
			tool.Managed = 0;
			ScriptGadget.Managed = 1;
			currentlyEditing = 0;
			selectedElementNum = 0;
			# Manage the dialog controls
			strikeFld.SetValue(0, 0);
			dipFld.SetValue(0, 0);
			dipDirFld.SetValue(0, 0);
			overturnedToggle.SetValue(0, 0);
			rGpMODE.SetEnabled(1);
			btnCLOSE.SetEnabled(1);
			promptString.SetLabel("Select point to Edit");
		# Else the point needs to be loaded and allowed to be modified
		} else {
			promptString.SetLabel("Loading Data...");
			btnACCEPT.SetEnabled(0);
			btnCANCEL.SetEnabled(0);
			overturnedToggle.SetEnabled(1);
			aLINECOLOR_LABEL.SetEnabled(1);
			btnLINECOLOR.SetEnabled(1);
			# Get the original 3 points
			local class POINT2D p1, p2, p3;
			p1.x = StrikeDipVector.point[selectedElementNum].Bedding.ToolPos1_x;
			p1.y = StrikeDipVector.point[selectedElementNum].Bedding.ToolPos1_y;
			p2.x = StrikeDipVector.point[selectedElementNum].Bedding.ToolPos2_x;
			p2.y = StrikeDipVector.point[selectedElementNum].Bedding.ToolPos2_y;
			p3.x = StrikeDipVector.point[selectedElementNum].Bedding.ToolPos3_x;
			p3.y = StrikeDipVector.point[selectedElementNum].Bedding.ToolPos3_y;
			# Map them to the screen coordinates
			local class TRANSPARM trans;
			trans = ViewGetTransMapToView(View, rasterLayer.MapRegion.CoordRefSys , 0);
			p1 = TransPoint2D(p1, trans);
			p2 = TransPoint2D(p2, trans);
			p3 = TransPoint2D(p3, trans);
			trans = ViewGetTransViewToScreen(View, 0);
			p1 = TransPoint2D(p1, trans);
			p2 = TransPoint2D(p2, trans);
			p3 = TransPoint2D(p3, trans);
			# Reset the polytool with the 3 points
			poly.Clear();
			poly.AppendVertex(p1);
			poly.AppendVertex(p2);
			poly.AppendVertex(p3);
			# Set the polytool for editing
			currentlyEditing = 1;
			tool.SetPolygon(poly);
			tool.Managed = 1;
			promptString.SetLabel("Modify the 3 points");
		}
	}
}
##################################################################################
#### Called when the user clicks on the Cancel button in the dialog.  Actions are
#### based on  the current Mode of the ToolScript
proc OnCancel( ) {
	strikeFld.SetValue(0, 0);
	dipFld.SetValue(0, 0);
	dipDirFld.SetValue(0, 0);
	overturnedToggle.SetValue(0, 0);
	# If in VIEW mode, unselect the point and reset the control dialog
	if (Current_Mode == "VIEW") {
		btnACCEPT.SetEnabled(0);
		btnCANCEL.SetEnabled(0);
		SDpoints.HighlightSingle(selectedElementNum, "Toggle");
		selectedElementNum = 0;
		rGpMODE.SetEnabled(1);
		btnCLOSE.SetEnabled(1);
		promptString.SetLabel("Select point to Toggle Data ON/OFF");
	# If in ADD mode, remove the temporary outcrop trace pattern and reset the dialog
	} else if (Current_Mode == "ADD") {
		overturnedToggle.SetEnabled(0);
		aLINECOLOR_LABEL.SetEnabled(0);
		btnLINECOLOR.SetEnabled(0);
		btnACCEPT.SetEnabled(0);
		btnCANCEL.SetEnabled(0);
		promptString.SetLabel("Clearing Data");
		# Remove the computed temporary outcrop trace pattern and refresh the display
		GroupRemoveLayer(firstGroup, intersectVlayer);
		CloseVector(intersectV);
		waitingForAddConfirmation = 0;
		ViewRedrawLayer(View, SDvectorLayer);
		ViewRedrawLayer(View, OPvectorLayer);
		# Reset the tool
		poly.Clear();
		tool.SetPolygon(poly);
		# Manage the dialog controls
		rGpMODE.SetEnabled(1);
		btnCLOSE.SetEnabled(1);
		promptString.SetLabel("Select 3 new points");
	# If in EDIT mode, cancel out of selecting a point, or remove the temporary data
	} else if (Current_Mode == "EDIT") {
		btnACCEPT.SetEnabled(0);
		btnCANCEL.SetEnabled(0);
		SDpoints.HighlightSingle(selectedElementNum, "Toggle");
		tool.Managed = 0;
		ScriptGadget.Managed = 1;
		poly.Clear();
		# Remove the temporary outcrop trace pattern data
		if (currentlyEditing) {
			GroupRemoveLayer(firstGroup, intersectVlayer);
			CloseVector(intersectV);
		}
		currentlyEditing = 0;
		waitingForAddConfirmation = 0;
		rGpMODE.SetEnabled(1);
		btnCLOSE.SetEnabled(1);
		promptString.SetLabel("Select point to Edit");
	# If in DELETE mode, deselect the point
	} else if (Current_Mode == "DELETE") {
		btnDELETE.SetEnabled(0);
		btnCANCEL.SetEnabled(0);
		SDpoints.HighlightSingle(selectedElementNum, "Toggle");
		selectedElementNum = 0;
		rGpMODE.SetEnabled(1);
		btnCLOSE.SetEnabled(1);
		promptString.SetLabel("Select point to Delete");
	}
	# Reset the selected element to nothing
	selectedElementNum = 0;
}
##################################################################################
#### Called when the user clicks on the Cancel button in the dialog.  Deletes the
#### selected point and all the records attached to that station.
proc OnDelete( ) {
	promptString.SetLabel("Deleting Data...");
	rGpMODE.SetEnabled(0);
	btnCLOSE.SetEnabled(0);
	btnCANCEL.SetEnabled(0);
	btnDELETE.SetEnabled(0);
	# Delete the station and attached records and update the display
	deleteStation();
	UpdateQuery();
	ViewRedrawLayer(View, SDvectorLayer);
	strikeFld.SetValue(0, 0);
	dipFld.SetValue(0, 0);
	dipDirFld.SetValue(0, 0);
	rGpMODE.SetEnabled(1);
	btnCLOSE.SetEnabled(1);
	promptString.SetLabel("Select point to Delete");
	selectedElementNum = 0;
}
##################################################################################
#### Called when the user clicks on the Close button to close the dialog
proc OnClose( ) {
	tool.Managed = 0;
	setDefaultWhenClose = true;
	View.SetDefaultTool();
}
##################################################################################
#### Called when the Right Button is pressed.  This is an alternative method of
#### clicking "Accept" in certain situations
proc OnRightButtonPress( ) {
	if (!tool.Managed) {
		if (selectedElementNum > 0) {
			if (Current_Mode == "VIEW") {
				OnAccept();
			} else if (Current_Mode == "DELETE") {
				OnDelete();
			}
		}
	}
}
##################################################################################
#### Called when the Left Button is pressed.  This selects the closest point in the
#### StrikeDipVector and displays the information related to that point.  Certain
#### modes use this selected point for their operations
proc OnLeftButtonPress( ) {
	# Make sure the polytool is not currently active or it will conflict with it
	if (!tool.Managed) {
		# Get the point on the screen that was clicked
		local class POINT2D screenPoint, layerPoint;
		screenPoint.x = PointerX;
		screenPoint.y = PointerY;
		# Get the layer coordinates from the screen point and find the closest point
		local class TRANSPARM transparm = ViewGetTransLayerToScreen(View, SDvectorLayer, 1);
		layerPoint = TransPoint2D(screenPoint, transparm);
		selectedElementNum = FindClosestPoint(StrikeDipVector, layerPoint.x, layerPoint.y, SDgeoref, 300);
		# If greater then 0, it has found the closest point
		if (selectedElementNum > 0) {
			SDpoints.HighlightSingle(selectedElementNum, "Replace");
			# Load the values for that point in the dialog
			strikeFld.SetValue(StrikeDipVector.point[selectedElementNum].Bedding.StrikeAngle, 0);
			dipFld.SetValue(StrikeDipVector.point[selectedElementNum].Bedding.DipAngle, 0);
			dipDirFld.SetValue(StrikeDipVector.point[selectedElementNum].Bedding.DipDirection, 0);
			overturnedToggle.SetEnabled(1);
			overturnedToggle.SetValue(StrikeDipVector.point[selectedElementNum].Bedding.Overturned, 0);
			overturnedToggle.SetEnabled(0);
			# Find the color attached to that station and update the color button
			local class COLOR tempColor;
			local numeric selectedStation = StrikeDipVector.point[selectedElementNum].Bedding.Station;
			local numeric numPoints = NumVectorPoints(OutcropPVector);
			local numeric i = 1;
			local numeric done = 0;
			while (!done) {
				# If it can find the station, load the color values
				if (OutcropPVector.point[i].Outcrop.Station == selectedStation) {
					tempColor.red = OutcropPVector.point[i].Outcrop.RedComp;
					tempColor.green = OutcropPVector.point[i].Outcrop.GreenComp;
					tempColor.blue = OutcropPVector.point[i].Outcrop.BlueComp;
					done = 1;
				# If not, load the color as Black and Notify user
				} else if (i == numPoints) {
					PopupMessage("Station Color Not Found");
					tempColor.red = 0;
					tempColor.green = 0;
					tempColor.blue = 0;
					done = 1;
				}
				i++;
			}
			btnLINECOLOR.SetColor(tempColor);
			# Update the control dialog settings based on the current mode
			if (Current_Mode == "VIEW") {
				promptString.SetLabel("CONFIRM: Toggle Data for this point?");
				rGpMODE.SetEnabled(0);
				btnACCEPT.SetEnabled(1);
				btnCANCEL.SetEnabled(1);
				btnCLOSE.SetEnabled(0);
			} else if (Current_Mode == "EDIT") {
				promptString.SetLabel("CONFIRM: Edit this point?");
				rGpMODE.SetEnabled(0);
				btnACCEPT.SetEnabled(1);
				btnCANCEL.SetEnabled(1);
				btnCLOSE.SetEnabled(0);
			} else if (Current_Mode == "DELETE") {
				promptString.SetLabel("CONFIRM: Delete this point?");
				rGpMODE.SetEnabled(0);
				btnCANCEL.SetEnabled(1);
				btnDELETE.SetEnabled(1);
				btnCLOSE.SetEnabled(0);
			}
		}
	}
}
##################################################################################
#### Called the first time the tool is activated.  The dialog is created here, but
#### not yet displayed
proc OnInitialize ( ) {
	# Is this a layout?  Set callbacks in case things change
	if (Layout) {
		WidgetAddCallback(Layout.GroupSelectedCallback, cbGroup);
		firstGroup = Layout.FirstGroup;
	} else {
		firstGroup = Group;
	}
	WidgetAddCallback(firstGroup.LayerSelectedCallback, cbLayer);
	# Is the bottom layer a valid DEM raster?
	if( checkLayer() ) {
		# Initialize tool
		tool = ViewCreatePolyLineGadget(View);
		ToolAddCallback(tool.ActivateCallback, cbToolApply);
		checkGeoref();
		# Create the destination vectors and their tables
		createDestVectors();
		addVectorsToDisplay();
		initializeBeddingDatabaseTable();
		initializeOutcropDatabaseTable();
		View.Redraw();
	} else {
		string message$ = "Active Layer: "+activeLayer.Name+" is invalid.";
		PopupMessage(message$);
	}
	# XML script that specifies the context of the dialog
	string xml$='<?xml version="1.0"?>
	<root>
		<dialog id="id_StrikeDipDialog" Title="StrikeDip Toolscript Control Panel" buttons="">
			<groupbox Name=" Current Mode: " ExtraBorder="5">
				<radiogroup id="id_MODE" orientation="Horizontal" Default="ADD" OnSelection="OnModeChange();">
					<item Value="VIEW" Name="View"/>
					<item Value="ADD" Name="Add"/>
					<item Value="EDIT" Name="Edit"/>
					<item Value="DELETE" Name="Delete"/>
				</radiogroup>
			</groupbox>
			<pane Orientation="Horizontal">
				<label>Strike Azimuth:   </label>
				<editnumber Width="6" id="id_STRIKEANGLE" Precision="0" ReadOnly="true" HorizResize="Fixed"/>
			</pane>
			<pane Orientation="Horizontal">
				<label>Dip Angle:   </label>
				<editnumber Width="6" id="id_DIPANGLE" Precision="0" ReadOnly="true" HorizResize="Fixed"/>
			</pane>
			<pane Orientation="Horizontal">
				<label>Dip Azimuth:   </label>
				<editnumber Width="6" id="id_DIPDIRECTION" Precision="0" ReadOnly="true" HorizResize="Fixed"/>
			</pane>
			<togglebutton id="id_OVERTURNED" name="Set as Overturned" OnChanged="OnOverturnedChanged();" Selected="false" Enabled="false"/>
			<pane Orientation="Horizontal">
				<label id="id_LINECOLOR_LABEL">Set Outcrop Trace Color:   </label>
				<colorbutton id="id_LINECOLOR"/>
			</pane>
			<groupbox Name=" Apply Action: " Orientation="Horizontal" ExtraBorder="5">
				<pushbutton id="id_ACCEPT" Name="Accept" OnPressed="OnAccept();"/>
				<pushbutton id="id_CANCEL" Name="Cancel" OnPressed="OnCancel();"/>
				<pushbutton id="id_DELETE" Name="Delete" OnPressed="OnDelete();"/>
				<pushbutton id="id_CLOSE" Name="Close" OnPressed="OnClose();"/>
			</groupbox>
			<label id="id_PROMPT" HorizResize="Expand">Select Prompt</label>
		</dialog>
	</root>';
	# Parse the document to make sure it is valid
	local numeric err = dlgdoc.Parse(xml$);
	if ( err < 0 ) {
		PopupError( err ); 	# Popup an error dialog. "Details" button shows syntax errors.
		Exit();
	}
	# Find the document
	node = dlgdoc.GetElementByID("id_StrikeDipDialog");
	if ( node == 0 ) {
		PopupMessage("Could not find dialog node in XML document");
		Exit();
	}
	# Set the XML Node
	err = dlg.SetXMLNode(node);
	if ( err < 0 ) {
		PopupMessage("Could not set the XML Node");
		Exit();
	}
	# Create a modeless dialog and assign some commonly changed elements to controls
	dlg.CreateModeless();
	rGpMODE = dlg.GetCtrlByID("id_MODE");
	btnLINECOLOR = dlg.GetCtrlByID("id_LINECOLOR");
	aLINECOLOR_LABEL = dlg.GetCtrlByID("id_LINECOLOR_LABEL");
	overturnedToggle = dlg.GetCtrlByID("id_OVERTURNED");
	promptString = dlg.GetCtrlByID("id_PROMPT");
	strikeFld = dlg.GetCtrlByID("id_STRIKEANGLE");
	dipFld = dlg.GetCtrlByID("id_DIPANGLE");
	dipDirFld = dlg.GetCtrlByID("id_DIPDIRECTION");
	btnACCEPT = dlg.GetCtrlByID("id_ACCEPT");
	btnCANCEL = dlg.GetCtrlByID("id_CANCEL");
	btnDELETE = dlg.GetCtrlByID("id_DELETE");
	btnCLOSE = dlg.GetCtrlByID("id_CLOSE");
	}
##################################################################################
#### Called when the tool is destroyed.  It cleans up resources used.
proc OnDestroy ( ) {
	tool.Managed = 0;
	CloseDatabase(SDdb);
	CloseDatabase(OPpointdb);
	CloseDatabase(OPlinedb);
	CloseVector(StrikeDipVector);
	CloseVector(OutcropPVector);
}
##################################################################################
#### Called when tool is activated.  It opens the dialog.
proc OnActivate ( ) {
	dlg.Open();
	OnModeChange();
	setDefaultWhenClose = true;
}
##################################################################################
#### Called when tool is deactivated.  Closes the dialog.
proc OnDeactivate ( ) {
	dlg.Close(0);
	setDefaultWhenClose = false;
}