VSHEDMOV.sml

#########################
#
# VSHEDMOV.SML
#
# Demonstration of movie generation from 2D view showing
# series of viewsheds computed from points derived by moving
# along a 2D path vector.
#
# Requires TNTmips Version 6.5
# 
##########################

clear();

#### Set movie format, frame rate, and recording time (seconds)

# Movie format (Possible values : "MPEG"(All platform) or "AVI"(Windows only) 
string format$;
format$ = "AVI";

# Frame rate 
# Possible values : 
#   "MOVIE_FRAMERATE_23_976"  23.976 (24000/1001) fps - NTSC encapsulated film rate
#   "MOVIE_FRAMERATE_24"		24 fps - Standard international cinema film rate
#   "MOVIE_FRAMERATE_25"   	25 fps - PAL (625/50) video frame rate
#   "MOVIE_FRAMERATE_29_970"  29.97 (30000/1001) fps - NTSC video frame rate
#   "MOVIE_FRAMERATE_30"   	30 fps - NTSC drop-frame (525/60) video frame rate
#   "MOVIE_FRAMERATE_50"   	50 fps - Double frame rate / progressive PAL
#   "MOVIE_FRAMERATE_59_940"  59.94 (60000/1001) fps - Double frame rate NTSC
#   "MOVIE_FRAMERATE_60"   	60 fps - Double frame rate drop-frame NTSC

string framerate$;
framerate$ = "MOVIE_FRAMERATE_24";

# Recording time
numeric time;
time = 3;

##### Get RVC objects to load

print("Select raster to use for elevation");
raster DEM, Viewshed;
vector VecLine;
GetInputRaster(DEM);
GetInputVector(VecLine);

numeric numLins;
numeric numCols;

numLins = NumLins(DEM);
numCols = NumCols(DEM);

CreateTempRaster(Viewshed,numLins,numCols,"binary",1);

##### Create color palette subobject for binary Viewshed raster
class RVC_COLORMAP vColorMap;

# Define one palette color with percent transparency
class Color vcolor;
vcolor.red = 100;
vcolor.green = 100;
vcolor.blue = 0;
vcolor.transp = 40;

# Write color palette to temporary viewshed raster and set index color 1
ColorMapWriteToRastVar(Viewshed,vColorMap,"ColorMap","");
ColorMapSetColor(vColorMap,1,vcolor);

##### RVC style object to draw viewer point

string styleFilename$;
string styleObjectname$;
GetInputObject("Style","Select style object for center and viewer point symbols:", 
	styleFilename$, styleObjectname$);
string viewer$;
viewer$ = "VIEWER";

######## Create display group and view and add DEM

print("START");

# Size of squared 2D view
# Should be evenly divisible by 8
numeric size;
size = 400;

# Zoom out factor for 2D view
numeric zoomfactor;
zoomfactor = 1.0;

# Create group
print("Creating Group");
class GRE_GROUP group;
group = GroupCreate();

# Create flags to create view without iconbar, scrollbars, status line and scale/position line
# This is important to maintain fixed window size for movie generation
string flags$;
flags$ = "NoScalePosLine,NoIconBar,NoScrollbars,NoStatusLine";

# Create dialog and 2D view
print("Creating dialog and 2D view");
class XmForm dialog2d;
class GRE_VIEW view2d;
dialog2d = CreateFormDialog("VIEW 2D");
view2d = GroupCreateView(group,dialog2d,"",size,size,flags$);
view2d.BackgroundColor.red = 67;
view2d.BackgroundColor.green = 100;
view2d.BackgroundColor.blue = 100;

# Add DEM to group
GroupQuickAddRasterVar(group,DEM);

# Open view and redraw full
DialogOpen(dialog2d);
ViewRedrawFull(view2d);
ViewZoomOut(view2d,zoomfactor,1);

####### Set up parameters for movie frame

# Destination of view in frame for movie
numeric x2d, y2d, w, h;
x2d = 0;
y2d = 0;
w = size;
h = size;

# Create frame
print("Creating frame for movie");
class Frame frame;
frame = FrameCreate(w,h);
numeric fontsize;

# Create graphics context (GC) for frame
print("Creating GC for frame and activate it");
ActivateGC(FrameCreateGC(frame));
DrawTextSetHeightPixels(fontsize);
DrawUseStyleObject(styleFilename$,styleObjectname$);

######## Set some more movie parameters

# Initialize Movie
print("Initializing Movie");
class Movie movie;
movie = MovieInit();

# Check framerate and force it to "MOVIE_FRAMERATE_24" if it is invalid
numeric rate;
rate = 24;
if (framerate$ == "MOVIE_FRAMERATE_23_976") rate = 23.976;
if (framerate$ == "MOVIE_FRAMERATE_25") rate = 25.0;
if (framerate$ == "MOVIE_FRAMERATE_29_970") rate = 29.970;
if (framerate$ == "MOVIE_FRAMERATE_30") rate = 30.0;
if (framerate$ == "MOVIE_FRAMERATE_50") rate = 50.0; 
if (framerate$ == "MOVIE_FRAMERATE_59_940") rate = 59.940; 
if (framerate$ == "MOVIE_FRAMERATE_60") rate = 60.0;
if (rate == 24.0) framerate$ = "MOVIE_FRAMERATE_24";

# Set Movie Parameters
print("Setting Movie Parameters");
MovieSetFormat(movie,format$);
MovieSetFrameRate(movie,framerate$);
MovieSetFrameWidth(movie,w);
MovieSetFrameHeight(movie,h);

# Make Output File
string ext$;
ext$ = MovieGetFileExt(movie);
string filename$;
filename$ = GetOutputFileName("","Make filename for movie",ext$);
printf("Filename = %s\n",filename$);

# Check recording time
if (time <= 1.0) time = 1.0;

# Calculate number of frames
numeric numFrames;
numFrames = time * rate;

####### Get georeference parameters for layers and reset to group projection
####### (defined by first raster layer)

class Georef vgeoref;
vgeoref = GetLastUsedGeorefObject(VecLine);
GeorefSetProjection(vgeoref,group.Projection);

class Georef rgeoref;
rgeoref = GetLastUsedGeorefObject(DEM);
GeorefSetProjection(rgeoref,group.Projection);

##### Compute resampled path for traverse from 2D vector line

# Create arrays to hold vertex coordinates for input traverse line
# Arrays will be automatically resized by GetVectorLinePointList()
numeric numvert;	# number of vertices in input line

array numeric xarrayL [1];	# array to hold x-coordinates
array numeric yarrayL [1];	# array to hold y-coordinates

# Get x and y coordinates of line vertices
numvert = GetVectorLinePointList(VecLine,xarrayL,yarrayL,1);

# Transform vertex coordinates for traverse line to map coordinates for group

array numeric zarrayL [numvert];	# array to hold z-coordinates of traverse line

numeric obj_x;
numeric obj_y;
numeric map_x;
numeric map_y;

numeric i;
for i = 1 to numvert {
	obj_x = xarrayL[i];
	obj_y = yarrayL[i];
	ObjectToMap(VecLine,obj_x,obj_y,vgeoref,map_x,map_y);
	xarrayL[i] = map_x;
	yarrayL[i] = map_y;
	} 

# Compute length of 2D traverse path and movement increment
numeric dx;
numeric dy;
numeric dl;
numeric length;		# path length in meters
numeric inclength;	# movement increment along traverse

length = 0;

for i = 1 to (numvert - 1) {
	dx = xarrayL[i+1] - xarrayL[i];
	dy = yarrayL[i+1] - yarrayL[i];
	dl = sqrt( dx*dx + dy*dy );
	length = length + dl;
}

inclength = length / (numFrames - 1);

# Compute equal-distance vertices for traverse path

array numeric xarray_eq [numFrames]; # array to hold x-coordinates of line vertices 
array numeric yarray_eq [numFrames]; # array to hold y-coordinates of line vertices

xarray_eq[1] = xarrayL[1];
yarray_eq[1] = yarrayL[1];

numeric numvert_eq;
numeric indexL;
numeric coef;
numeric sum;
numeric rest;

numvert_eq = 1;
indexL = 1;
rest = 0.0;

while (numvert_eq < (numFrames - 1)) {
	dx = xarrayL[indexL+1] - xarrayL[indexL];
	dy = yarrayL[indexL+1] - yarrayL[indexL];
	dl = sqrt( dx*dx + dy*dy );
	rest = rest + dl;
	sum = 0.0;
	while ((rest > inclength) and ( numvert_eq < (numFrames - 1) ) ) {
		rest = rest - inclength;
		sum = sum + inclength;
		numvert_eq = numvert_eq + 1;
		coef = sum / dl;
		xarray_eq[numvert_eq] = xarrayL[indexL] + coef * dx;
		yarray_eq[numvert_eq] = yarrayL[indexL] + coef * dy;
		}
	indexL = indexL + 1;
	}

xarray_eq[numFrames] = xarrayL[numvert];
yarray_eq[numFrames] = yarrayL[numvert];

#### Transform equidistant vertex map coordinates to raster and screen coordinates

array numeric xarrayR [numFrames];
array numeric yarrayR [numFrames];

array numeric xarrayS [numFrames];
array numeric yarrayS [numFrames];

class POINT2D point;

for i = 1 to numFrames {
	map_x = xarray_eq[i];
	map_y = yarray_eq[i];
	MapToObject(rgeoref,map_x,map_y,DEM,obj_x,obj_y);
	xarrayR[i] = obj_x;
	yarrayR[i] = obj_y;

	point.x = map_x;
	point.y = map_y;
	point = TransPoint2D(point,ViewGetTransMapToView(view2d,group.Projection));
	point = TransPoint2D(point,ViewGetTransViewToScreen(view2d));

	xarrayS[i] = point.x;
	yarrayS[i] = point.y;
	} 

##### Set up initial parameters for viewshed computation
numeric nPoints;		# number of view points
numeric percent;		# percentage of view points that must see a location point
numeric height;		# height of view point above surface
numeric zScale;		# vertical scale factor for elevation

nPoints = 1;			# 1 view point per frame
percent = 100;
height = 2;				# height above surface for each location
zScale = 1;

# Create matrices to hold coordinates for view point(s)
class MATRIX matX, matY;
matX = CreateMatrix(1,1); 
matY = CreateMatrix(1,1); 

numeric rastX;			# x-coordinate in DEM in object coordinates
numeric rastY;			# y-coordinate in DEM in object coordinates

##### Start recording movie to file

MovieStart(movie,filename$);

# Loop for each frame
for i = 1 to numFrames {
	SetStatusMessage(sprintf("Processing frame %d of %d",i,numFrames));

	### Compute viewshed

	rastX = xarrayR[i]; 	# Get X and Y coordinates of current view point
	rastY = yarrayR[i];

	SetMatrixItem(matX,0,0,rastX);	# Assign coordinates to matrices
	SetMatrixItem(matY,0,0,rastY);

	RasterToBinaryViewshed(DEM,Viewshed,nPoints,matX,matY,percent,height,zScale);
	CopySubobjects(DEM,Viewshed,"georef");

	### Add viewshed to group
	class GRE_LAYER_RASTER vslayer;	 # define raster layer class variable for viewshed raster

	vslayer = GroupQuickAddRasterVar(group,Viewshed);
	vslayer.NullCellsTransparent = 1;
	vslayer.ColorMap = vColorMap;			# use previously-created color palette

	vslayer.UsesTransparency = 1;			# turn on transparency effects for layer

	# ViewRedrawDirect(view2d,"NoBlankScreen")
	# This new function	added after release of TNTmips 6.5
	# can redraw the view without blanking it first.  Use of
	# this function eliminates "flashing" of the view as the
	# movie is initially rendered.  It has no effect on the
	# output movie file.  For 6.5 release version, use the
	# function in the next statement.

	ViewRedraw(view2d);

	# Copy 2d view to destination frame
	FrameCopyFromView(frame,view2d,0,0,size,size,x2d,y2d);

	# Draw viewer position point in 2d frame
	DrawSetPointStyle(viewer$);
	DrawPoint(xarrayS[i],yarrayS[i]);

	# Add frame to movie
	MovieAddFrame(movie,frame);

	# Destroy viewshed layer
	LayerDestroy(vslayer);
}									# End main processing loop

# Stop and Exit movie
MovieStop(movie);
MovieExit(movie);

# Close dialogs
DialogClose(dialog2d);

print("END");