Given a terrain, how to draw the stream flow path?

There are different possible implementations, but most procedures will start from a grid and not from a TIN.

The simplest one is most likely the D8 procedure: you calculate the direction where water would be flowing. There are 8 possibilities, the 8 cells that are next to a central grid cell. You can first calculate these directions, than how cells are connected, and finally you can draw the lines). An easy implementation is found in SAGA, it almost reads as pseudocode: http://saga-gis.svn.sourceforge.net/viewvc/saga-gis/trunk/saga-gis/src/modules_terrain_analysis/terrain_analysis/ta_channels/D8_Flow_Analysis.cpp?revision=911&view=markup

Although very easy, this is not very realistic: you will not have a stream starting in every cell. More advanced algorithms usually first close the pits (especially if you have a detailed DEM), then calculate the catchment area per cell, that is the number of cells that contribute water to a particular cell, and then use a threshold to determine whether a stream is present.

SAGA GIS implements a lot of these catchment area methods, you can find a description of them in this manual http://sourceforge.net/settings/mirror_choices?projectname=saga-gis&filename=SAGA%20-%20Documentation/SAGA%20Documents/SagaManual.pdf

It was written for an older version of SAGA GIS, but the description of the algorithms is still rather accurate, and I'll copy it here for quick reference (this is around page 120), since it is open source, you can check the implementation details by looking at the code.

• Deterministic 8 (D8): The classical. Flow goes from the center of a cell to the center of one(and only one) of the surrounding cells. Flow directions are, therefore, restricted to multiples of 45o , which is the main reason for most of the drawbacks of the method. (O’Callaghan & Mark 1984).
• Rho8: Same as above but with an stochastic component that should improve it. The flow direction is determined by a random argument that is dependent on the difference between aspect and the direction of the two adjacent neighbor cells. Not very useful. . . (Fairfield & Leymarie 1991).
• Deterministic infinity (D∞): Flow goes from one cell to two contiguous surrounding cells, thus considering a bidimensional flow and overcoming the drawbacks of the D8 method. (Tarboton 1998).
• Braunschweiger Digitales Reliefmodell: Another Multiple Flow Direction algo- rithm. The flow is split between the surrounding cell whose orientation is nearest to the aspect of the center cell and its two adjacent cells. (Bauer, Bork & Rohdenburg 1985).
• FD8 (found in SAGA simply as Multiple Flow Direction): A D8–derived bidimensional flow routing algorithm. (Quinn et al 1991).
• Kinematic Routing Algorithm (KRA). An unidimensional flow tracing algorithm. Flow behaves like a ball rolling down the DEM, without restricting its position to the center of cells. (Lea 1992).
• Digital Elevation Model Network (DEMON): The most complex one. A bidimensional flow tracing algorithm. Rather time consuming. (Costa-Cabral & Burgess 1994).

Even more models have been added recently:

• Triangular Multiple Flow Direction- Seibert, J. / McGlynn, B. (2007): 'A new triangular multiple flow direction algorithm for computing upslope areas from gridded digital elevation models', Water Resources Research, Vol. 43, W04501. This may be interesting for you because it might also work directly on a TIN
• The Mass-Flux Method (MFM) for the DEM based calculation of flow accumulation as proposed by Gruber and Peckham (2008). Gruber, S., Peckham, S. (2008): Land-Surface Parameters and Objects in Hydrology. In: Hengl, T. and Reuter, H.I. [Eds.]: Geomorphometry: Concepts, Software, Applications. Developments in Soil Science, Elsevier, Bd.33, S.293-308.
• The side algorithm: http://watershed.montana.edu/Hydrology/Home_files/2010WR009296.pdf and his code is on his website as well: http://thomasgrabs.com/side-algorithm/

An original approach is the one proposed in this paper:

Fisher, P., J. Wood, and T. Cheng (2004). Where is Helvellyn? Fuzziness of multiscale landscape morphometry. Transactions of the Institute of British Geographers 29, 106-128.

It proposes a method based on fuzzy and multi scale representation. I am not sure but this method may be the one implemented in LandSerf.

If you have access to Spatial Analyst in ArcGIS, then you have a series of tools to compute stream paths. A full workflow is provided in the ESRI reference, but the typical workflow includes:

1. Convert your TIN to an elevation raster.
2. Calculate the Flow Direction.
3. Fill small Sinks.
4. Calculate Flow Accumulation
5. Using a threshold, choose only the cells with a given amount of flow.
6. Use the Stream to Feature tool to export the streams to a vector shapefile.

Of course, there are numerous academic papers describing different methods, but this method is easy for all who have access to Spatial Analyst.