When running an OLTP system you may have a transactions table with things like transaction ids, transaction types and states of transaction (open, closed, cancelled, etc). Different transactions tend to have different usage patterns – some might be automatically dealt with, others require user interaction, some typically rejected and still others where a lot of mistakes are common so they get cancelled a lot. There are differences in the in-flight and final states between the various transactions but since they are mixed into the one transaction table the skew in the states can be obscured when looking at only a single column. The skew may be only visible when both the type and state are considered together – a situation for which extended stats is perfect for. To test the effect of multi-column histograms or EXTENDED_STATS in this type of scenario we’ll need a sample table which shows a skew only when multiple columns are considered. I’ve got type A transactions which typically have a final state of closed and type B transactions which typically end up getting cancelled. However a few As (10) do get cancelled and a few Bs (10) make it to a closed state. I have no open or other state transactions.

The test scenario is:

1. Create a transaction table with skewed data, transtype A mostly closed, B mostly cancelled.
2. Take a look at the data to verify that it has the expected skew. Note that there are 10 cancelled A transactions.
3. Explain plan for filter(”TRANSTYPE”=’A’ AND “STATE”=’Cancelled’) and note the expected row count: 500, which is Oracle’s best guess based on information available in normal single column histograms.
4. Try again but with a dynamic_sampling level of 4 and see the correct expected row count (10).
5. Generate multi-column stats on ‘(transtype,state)’
6. Retry the query without dynamic sampling and see that Oracle 11g is able to use the new stats to correctly estimate the row count of 10 and uses an index in preference to a table scan.

Obviously this is a contrived example which deliberately obscures the skew in the data with perfect balance between the two transaction types but with a larger number of transaction types it is easy for such patterns to be muddied. With the aid of a histogram on the two critical columns here Oracle was able to see the situation accurately and modify the execution plan accordingly going from a full table scan to an index lookup. A more accurate picture of the data will result in superior execution plans. Using the dynamic_sampling hint results in the exact same benefits. Oracle’s documentation on dynamic sampling levels says of level 4:

Level 4: Apply dynamic sampling to all tables that meet Level 3 criteria, plus all tables that have single-table predicates that reference 2 or more columns. The number of blocks sampled is the default number of dynamic sampling blocks. For unanalyzed tables, the number of blocks sampled is two times the default number of dynamic sampling blocks.

Note: The “estimate” mentioned in step 6 is the exact number because for the stats gathering parameter estimate_percent I’ve used NULL or compute which may not be practical in a production system. You’ll need to use a percentage which balances accurately modelling the data against stats gathering time. This interesting piece on that topic comes to the conclusion that between 5-20% is the right number but the analysis extends only as far as 10g.

Here are the results of running the SQL script mcol_histogram which you can download from the link up near the start of this post:

SQL> set linesize 100 echo on
SQL> @mcol_histogram
SQL> -- Step 1: Create a transaction table with skewed data, transtype A mostly closed, B mostly cancelled
SQL> CREATE TABLE transactions AS
  2  SELECT tmp.*, (CASE WHEN (tmp.transtype = 'A' AND MOD(transid,100) != 0) OR
  3                           (tmp.transtype = 'B' AND MOD(transid,100)  = 0) THEN 'Closed'
  4                      ELSE 'Cancelled'
  5                  END) state FROM
  6      (
  7      SELECT LEVEL transid, DECODE(FLOOR(LEVEL/1001),0,'A','B') transtype
  8        FROM dual
  9      CONNECT BY LEVEL <= 2000
 10      ) tmp;

Table created.

SQL>
SQL> CREATE INDEX ix1transactions ON transactions(transtype, state);

Index created.

SQL>
SQL> EXEC DBMS_STATS.GATHER_TABLE_STATS(ownname => SYS_CONTEXT('USERENV','CURRENT_USER'),
tabname => 'transactions', estimate_percent => NULL, method_opt => 'FOR ALL COLUMNS SIZE SKEWONLY', cascade => true);

PL/SQL procedure successfully completed.

SQL>
SQL> -- Step 2: Take a look at the data to verify that it has the expected skew.
SQL> SELECT transtype, state, COUNT(*)
  2    FROM transactions
  3  GROUP BY transtype, state
  4  ORDER BY transtype, state
  5    ;

T STATE       COUNT(*)
- --------- ----------
A Cancelled         10
A Closed           990
B Cancelled        990
B Closed            10

SQL>
SQL> -- Step 3:  This only returns 10 rows, but we have a full table scan and high cardinality
SQL> EXPLAIN PLAN FOR
  2  SELECT *
  3    FROM transactions t
  4   WHERE transtype = 'A'
  5     AND state = 'Cancelled';

Explained.

SQL>
SQL> SELECT * FROM TABLE(DBMS_XPLAN.DISPLAY);

PLAN_TABLE_OUTPUT
----------------------------------------------------------------------------------------------------
Plan hash value: 1321366336

----------------------------------------------------------------------------------
| Id  | Operation         | Name         | Rows  | Bytes | Cost (%CPU)| Time     |
----------------------------------------------------------------------------------
|   0 | SELECT STATEMENT  |              |   500 |  7500 |     3   (0)| 00:00:45 |
|*  1 |  TABLE ACCESS FULL| TRANSACTIONS |   500 |  7500 |     3   (0)| 00:00:45 |
----------------------------------------------------------------------------------

Predicate Information (identified by operation id):
---------------------------------------------------

PLAN_TABLE_OUTPUT
----------------------------------------------------------------------------------------------------

   1 - filter("TRANSTYPE"='A' AND "STATE"='Cancelled')

13 rows selected.

SQL>
SQL> -- Step 4:  Try again but use dynamic sampling this time and get much better results
SQL> EXPLAIN PLAN FOR
  2  SELECT /*+ dynamic_sampling (4) */ *
  3    FROM transactions t
  4   WHERE transtype = 'A'
  5     AND state = 'Cancelled';

Explained.

SQL>
SQL> SELECT * FROM TABLE(DBMS_XPLAN.DISPLAY);

PLAN_TABLE_OUTPUT
----------------------------------------------------------------------------------------------------
Plan hash value: 278489313

-----------------------------------------------------------------------------------------------
| Id  | Operation                   | Name            | Rows  | Bytes | Cost (%CPU)| Time     |
-----------------------------------------------------------------------------------------------
|   0 | SELECT STATEMENT            |                 |    10 |   150 |     2   (0)| 00:00:30 |
|   1 |  TABLE ACCESS BY INDEX ROWID| TRANSACTIONS    |    10 |   150 |     2   (0)| 00:00:30 |
|*  2 |   INDEX RANGE SCAN          | IX1TRANSACTIONS |    10 |       |     1   (0)| 00:00:15 |
-----------------------------------------------------------------------------------------------

Predicate Information (identified by operation id):

PLAN_TABLE_OUTPUT
----------------------------------------------------------------------------------------------------
---------------------------------------------------

   2 - access("TRANSTYPE"='A' AND "STATE"='Cancelled')

Note
-----
   - dynamic sampling used for this statement (level=4)

18 rows selected.

SQL>
SQL> -- Step 5:  Generate some multi-column stats
SQL> SELECT DBMS_STATS.CREATE_EXTENDED_STATS (SYS_CONTEXT('USERENV','CURRENT_USER'),'transactions',
'(transtype,state)') FROM dual;

DBMS_STATS.CREATE_EXTENDED_STATS(SYS_CONTEXT('USERENV','CURRENT_USER'),'TRANSACTIONS','(TRANSTYPE,ST
----------------------------------------------------------------------------------------------------
SYS_STUMWLR#03VMK70K_1F3J$2ZMZ

SQL>
SQL> EXEC DBMS_STATS.GATHER_TABLE_STATS(ownname => SYS_CONTEXT('USERENV','CURRENT_USER'),
tabname => 'transactions', estimate_percent => NULL, method_opt => 'FOR ALL COLUMNS SIZE SKEWONLY', cascade => true);

PL/SQL procedure successfully completed.

SQL>
SQL> -- Step 6:  The new extended stats yield a cardinality of 10 and use an index lookup
SQL> EXPLAIN PLAN FOR
  2  SELECT *
  3    FROM transactions t
  4   WHERE transtype = 'A'
  5     AND state = 'Cancelled';

Explained.

SQL>
SQL> SELECT * FROM TABLE(DBMS_XPLAN.DISPLAY);

PLAN_TABLE_OUTPUT
----------------------------------------------------------------------------------------------------
Plan hash value: 278489313

-----------------------------------------------------------------------------------------------
| Id  | Operation                   | Name            | Rows  | Bytes | Cost (%CPU)| Time     |
-----------------------------------------------------------------------------------------------
|   0 | SELECT STATEMENT            |                 |    10 |   150 |     2   (0)| 00:00:30 |
|   1 |  TABLE ACCESS BY INDEX ROWID| TRANSACTIONS    |    10 |   150 |     2   (0)| 00:00:30 |
|*  2 |   INDEX RANGE SCAN          | IX1TRANSACTIONS |    10 |       |     1   (0)| 00:00:15 |
-----------------------------------------------------------------------------------------------

Predicate Information (identified by operation id):

PLAN_TABLE_OUTPUT
----------------------------------------------------------------------------------------------------
---------------------------------------------------

   2 - access("TRANSTYPE"='A' AND "STATE"='Cancelled')

14 rows selected.

SQL>
SQL> DROP TABLE transactions;

Table dropped.

SQL>
SQL> SELECT * FROM v$version;

BANNER
--------------------------------------------------------------------------------
Oracle Database 11g Enterprise Edition Release 11.2.0.1.0 - Production
PL/SQL Release 11.2.0.1.0 - Production
CORE    11.2.0.1.0      Production
TNS for Linux: Version 11.2.0.1.0 - Production
NLSRTL Version 11.2.0.1.0 - Production

SQL>

Final Note: The new stats are difficult to understand in the dba_histogram/user_histogram views with the endpoint_value being encoded in a different format to normal stats. I tried running Tom Kyte’s hexstr to decode it but was unsuccessful. Perhaps someone else will have better luck. Note that only the first 6 characters of endpoint_value can be decoded for normal histogram columns so it is perhaps better to stick to 6 letter transaction states rather than using things like ‘cancelled’.

If you are in a mutating table situation where your trigger needs to read from the same table that it is firing on your only real option before Oracle 11g was a multi trigger solution with a muddle of row and statement triggers communicating via a package. But now Oracle 11g comes through with a great new feature: compound triggers which neatly solve the problem … or do they. Well they purport to do so but it isn’t all smooth sailing: Beware the compound trigger. But I figured – that’s 11.1.0.6 whereas I’m using 11.1.0.7 and those particular failure modes don’t apply to my situation so I’ll be fine. Not exactly. What I have is a normal after statement trigger firing on the “banner” column and modifying data in the “changed” column and I need to fire a second trigger on that second column. What I saw was that the compound trigger worked properly once and then not again until it was recompiled.

trg_after => fires on column “banner” and modifies column “changed”.
trg_compound_row => fires on column “changed”
trg_simple_row => fires on same criteria as compound trigger for comparative purposes

The test scenario is:
0) Create table + triggers
1) Update table to see triggers firing then rollback to reset
2) Do step 1 again
3) Recompile trigger
4) Do step 1 twice
5) Drop table

SQL> set echo on
SQL> @trigger.sql
SQL> SET SERVEROUT ON LINESIZE 200
SQL>
SQL> CREATE table trigger_test AS
  2  SELECT rownum id, banner, 'A' changed FROM v$version;

Table created.

SQL>
SQL> CREATE OR REPLACE TRIGGER trg_after
  2      AFTER UPDATE OF banner ON trigger_test
  3  BEGIN
  4      dbms_output.put_line('Update');
  5      UPDATE trigger_test SET changed = 'B'
  6       WHERE id = 1;
  7  END;
  8  /

Trigger created.

SQL> CREATE OR REPLACE TRIGGER trg_compound_row
  2    FOR UPDATE OF changed ON trigger_test
  3  COMPOUND TRIGGER
  4
  5  AFTER EACH ROW IS BEGIN
  6      dbms_output.put_line('Compound fired <== Should happen every time');
  7  END AFTER EACH ROW;
  8
  9  END trg_compound_row;
 10  /

Trigger created.

SQL> CREATE OR REPLACE TRIGGER trg_simple_row
  2      AFTER UPDATE OF changed ON trigger_test
  3      FOR EACH ROW
  4  BEGIN
  5      dbms_output.put_line('Simple fired');
  6  END;
  7  /

Trigger created.

SQL> SELECT * FROM trigger_test;

        ID BANNER                                                                           C
---------- -------------------------------------------------------------------------------- -
         1 Oracle Database 11g Enterprise Edition Release 11.1.0.7.0 - 64bit Production     A
         2 PL/SQL Release 11.1.0.7.0 - Production                                           A
         3 CORE 11.1.0.7.0      Production                                                       A
         4 TNS for Solaris: Version 11.1.0.7.0 - Production                                 A
         5 NLSRTL Version 11.1.0.7.0 - Production                                           A

SQL>
SQL> UPDATE trigger_test SET banner = banner || ' ' WHERE id = 2;
Update
Compound fired <== Should happen every time
Simple fired

1 row updated.

SQL> ROLLBACK;

Rollback complete.

SQL>
SQL> UPDATE trigger_test SET banner = banner || ' ' WHERE id = 2;
Update
Simple fired

1 row updated.

SQL> ROLLBACK;

Rollback complete.

SQL>
SQL> ALTER TRIGGER trg_compound_row COMPILE;

Trigger altered.

SQL>
SQL> UPDATE trigger_test SET banner = banner || ' ' WHERE id = 2;
Update
Compound fired <== Should happen every time
Simple fired

1 row updated.

SQL> ROLLBACK;

Rollback complete.

SQL>
SQL> UPDATE trigger_test SET banner = banner || ' ' WHERE id = 2;
Update
Simple fired

1 row updated.

SQL> ROLLBACK;

Rollback complete.

SQL>
SQL> DROP TABLE trigger_test;

Table dropped.

SQL>

So the lesson is clear, as Dominic Brooks already stated: Beware the compound trigger.

You can download the source code for this post here 3DViewerBingMaps.zip

In OpenGL you can draw something called a Stippled Line which is basically a dashed line with some extra control over the dash pattern, however in Bing Maps 3D there is no such thing, which is why I wrote my own dashed line mesh class. My class supports the following things:

- Lines are made up of multiple segments
- Dashes and gaps between dashes are independent in size
- You can set an offset each frame allowing for a moving ant effect
- You can set different widths for the start and the end of the line

Here’s 2 pictures of the same line showing the billboard effect of the wide dashes, and the constant dash and gap size:

If you download the sample at the top of this post you can run the application shown in the screenshots.

Here’s the code snippet for setting up the above Dashed Mesh:

DashedLineMesh dashedLine = new DashedLineMesh();

dashedLine.DashLength = 8;
dashedLine.GapLength = 16;
dashedLine.LineStartWidth = 1;
dashedLine.LineEndWidth = 64;
dashedLine.Color = Color.White;
dashedLine.ZBufferEnable = true;

dashedLine.Points.Clear();
GeodeticViewpoint centrePoint = new GeodeticViewpoint();
double alt = 1000000;
double lat = 45.0;
double lon = 180.0;

for(double i = 0; i<(Math.PI*2); i+=0.103)
{
  double x = 8+30 * Math.Sin(i);
  double y = 8+30 * Math.Cos(i);
  centrePoint.Position.Location = LatLonAlt.CreateUsingDegrees(lat+x, lon+y, alt);
  dashedLine.Points.Add(new Vector3F(centrePoint.Position.Vector));
}

and as a bonus here’s the complete source code for the DashedLineMesh class

// Dashed Line Mesh by John Stewien 2010

using System;
using System.Collections.Generic;
using System.Text;
using Microsoft.MapPoint.Rendering3D.GraphicsProxy;
using Microsoft.MapPoint.Geometry.VectorMath;
using Microsoft.MapPoint.Geometry;
using Microsoft.MapPoint.Graphics3D;
using Microsoft.MapPoint.Rendering3D.State;

namespace _3DViewerBingMaps.Plugin {

  /// <summary>
  /// Draws a multi-segment line that can be broken up into dashes
  /// </summary>
  public class DashedLineMesh : MeshGraphicsObject<Vertex.PositionOnly, ushort> {

    // ************************************************************************
    // Private Fields
    // ************************************************************************
    #region Private Fields

    /// <summary>
    /// A list of points that make up the line to be drawn. Points are in
    /// Lat, Lon, Alt format.
    /// </summary>
    private List<Vector3F> points;
    /// <summary>
    /// The starting width of the line in pixels
    /// </summary>
    private float lineStartWidth = 1f;
    /// <summary>
    /// The ending width of the line in pixels
    /// </summary>
    private float lineEndWidth = 1f;
    /// <summary>
    /// The length of the dashes in pixels
    /// </summary>
    private float dashLength = 0f;
    /// <summary>
    /// the length of the gaps between the dashes in pixels
    /// </summary>
    private float gapLength = 0f;
    /// <summary>
    /// The offset of the dahses in pixels. Incrementing this by .5 every frame
    /// gives a moving and effect.
    /// </summary>
    private float dashOffset = 0f;
    /// <summary>
    /// Enable alpha blending for this mesh
    /// </summary>
    private bool alphaEnable = false;
    /// <summary>
    /// Enable using the Depth Buffer for this mesh
    /// </summary>
    private bool zBufferEnable = true;
    /// <summary>
    /// The color to draw this line in
    /// </summary>
    private System.Drawing.Color color = System.Drawing.Color.White;

    #endregion Private Fields

    // ************************************************************************
    // Public Methods
    // ************************************************************************
    #region Public Methods

    /// <summary>
    /// Constructor
    /// </summary>
    /// <param name="radius"></param>
    /// <param name="material"></param>
    /// <param name="materialId"></param>
    public DashedLineMesh()
      : base(GraphicsBufferUsage.Dynamic, GraphicsBufferUsage.Dynamic, 4, 4, true) {

      this.points = new List<Vector3F>();

      if (base.State.Materials.Count == 0) {
        base.State.Materials.Add(new Material());
      }
      base.ProjectionEnabled = false;
      this.UpdateRenderState();
    }

    /// <summary>
    /// Updates the graphic prior to rendering it in a frame. Needs to be
    /// called every frame.
    /// </summary>
    /// <param name="camera"></param>
    /// <param name="transform"></param>
    /// <remarks>
    /// Unfortunately the GraphicsObject.UpdateInternalRenderObject method
    /// can't be overriden because of its "internal" protection, so this
    /// needs to be called every frame before the graphics is rendered
    /// </remarks>
    public void PreRenderUpdate(CameraData camera, Matrix4x4D transform) {
      int width = camera.View.Size.Width;
      int height = camera.View.Size.Height;
      this.UpdateLineGeometry(width, height, transform*camera.Snapshot.TransformMatrix);
    }

    #endregion Public Methods

    // ************************************************************************
    // Private Methods
    // ************************************************************************
    #region Private Methods

    /// <summary>
    /// Calculates the required line width from the user specified value and the
    /// size of the display area
    /// </summary>
    /// <param name="deviceWidth"></param>
    /// <param name="deviceHeight"></param>
    /// <returns></returns>
    private Vector3F GetWidth(float lineWidth, int deviceWidth, int deviceHeight) {
      float x = lineWidth / ((float)deviceWidth);
      float y = lineWidth / ((float)deviceHeight);
      Vector3F retVal = new Vector3F(x, y, 1f);
      return retVal;
    }

    /// <summary>
    /// Apply the position/scale tranform, and the projection transform to the points
    /// to convert them to screen coordinates
    /// </summary>
    /// <param name="projectionTransform"></param>
    /// <returns></returns>
    private IList<Vector3F> TransformedPoints(Matrix4x4D projectionTransform) {
      List<Vector3F> list = new List<Vector3F>(this.points.Count);
      Matrix4x4D matrix = projectionTransform;
      foreach (Vector3F vectorf in this.points) {
        Vector4D vector = new Vector4D((double)vectorf.X, (double)vectorf.Y, (double)vectorf.Z, 1.0);
        vector = Vector4D.TransformCoordinate(ref vector, ref matrix);
        if (vector.W > 0.0) {
          vector.MultiplyBy(1.0 / vector.W);
        } else {
          vector.MultiplyBy(1.0 / -vector.W);
        }
        list.Add(new Vector3F((float)vector.X, (float)vector.Y, (float)vector.Z));
      }
      return list;
    }

    /// <summary>
    /// This method takes all the solid lines and breaks them up into dashes
    /// </summary>
    /// <param name="lineList"></param>
    /// <param name="widthList"></param>
    /// <param name="deviceWidth"></param>
    /// <param name="deviceHeight"></param>
    private bool ConvertToDashed(ref IList<Vector3F> lineList, ref IList<float> widthList, int deviceWidth, int deviceHeight) {
      this.dashOffset = dashOffset % (dashLength + gapLength);

      // Scale everything to viewport space
      float dashStep = 2 * this.dashLength / Math.Min(deviceWidth, deviceHeight);
      float gapStep = 2 * this.gapLength / Math.Min(deviceWidth, deviceHeight);
      float offsetStep = 2 * this.dashOffset / Math.Min(deviceWidth, deviceHeight);

      // Sanity check dashStep and gapStep
      if (dashStep <= 0 || gapStep <= 0) {
        return false;
      }

      // Calculate the initial conditions
      float remainder = offsetStep;
      bool doDash = remainder==0;
      if (remainder > gapStep) {
        remainder -= gapStep;
        doDash = true;
      }

      bool doingJoin = false;

      List<Vector3F> dashedLineList = new List<Vector3F>();
      List<float> dashedWidthList = new List<float>();

      for (int i = 0; i < (lineList.Count - 1); i++) {
        Vector3F fromPoint = lineList[i];
        Vector3F toPoint = lineList[i + 1];

        // Store the original line length. Used for calculating the line width later on.
        float lineLength = Vector3F.Distance(fromPoint, toPoint);

        // Clamp the line to the visible area otherwise a lot more calculation will be done
        // than is required
        if (ClampLineToBox(ref fromPoint, ref toPoint, new Vector2F(-1, -1), new Vector2F(1, 1))) {

          // Need to break up the line in the X-Y direction regardless of Z
          Vector2F segment = new Vector2F(toPoint.X - fromPoint.X, toPoint.Y - fromPoint.Y);
          float segLength = segment.Length();
          if (segLength > 0) {

            float widthDiff = widthList[i + 1] - widthList[i];

            Vector3F step = new Vector3F(segment.X, segment.Y, toPoint.Z - fromPoint.Z) / segLength;

            float distTravelled = 0;
            float nextDash = remainder > 0 ? remainder : dashStep;
            float nextgap = remainder > 0 ? remainder : gapStep;

            // Calcluate the line width for the point we are starting at and
            // Calculate the line width increments for when we iterate over the dashes

            float invLineLength = 1f / lineLength;
            float progressiveWidth = widthList[i] + widthDiff * (remainder + Vector3F.Distance(fromPoint, lineList[i])) * invLineLength;
            float dashDWidth = widthDiff * dashStep * invLineLength;
            float gapDWidth = widthDiff * gapStep * invLineLength;

            // Keep on adding dashes until we reach the end of the line segment
            // Also calculate the line widths for each end of the dashes

            while (distTravelled < segLength) {
              if (doDash) {
                if (doingJoin) {
                  dashedLineList.RemoveAt(dashedLineList.Count - 1);
                  dashedWidthList.RemoveAt(dashedWidthList.Count - 1);
                } else {
                  dashedLineList.Add(distTravelled * step + fromPoint);
                  dashedWidthList.Add(progressiveWidth);
                }

                doingJoin = false;

                distTravelled += nextDash;
                progressiveWidth += dashDWidth;

                dashedLineList.Add(Math.Min(segLength, distTravelled) * step + fromPoint);
                dashedWidthList.Add(progressiveWidth);

                if (distTravelled > segLength && i < (lineList.Count - 2)) {
                  doingJoin = true;
                }
              } else {
                distTravelled += nextgap;
                progressiveWidth += gapDWidth;
              }
              doDash = !doDash;
              nextDash = dashStep;
              nextgap = gapStep;
            }
            remainder = distTravelled - segLength;
            if (remainder > 0) {
              doDash = !doDash;
            }
          }
        } else {
          doingJoin = false;
        }
      }
      lineList = dashedLineList;
      widthList = dashedWidthList;
      return true;
    }

    /// <summary>
    /// Clamps the line between the 2 vectors passed to be within the box limits passed in.
    /// </summary>
    /// <param name="from"></param>
    /// <param name="to"></param>
    /// <param name="min"></param>
    /// <param name="max"></param>
    /// <returns>
    /// true if the line is visible, false otherwise
    /// </returns>
    private bool ClampLineToBox(ref Vector3F from, ref Vector3F to, Vector2F boxMin, Vector2F boxMax) {
      // need to limit everything to within -1 to 1 on all dimensions

      // Test if anything is inside the viewbox
      if (Math.Min(from.X, to.X) > boxMax.X || Math.Max(from.X, to.X) < boxMin.X ||
        Math.Min(from.Y, to.Y) > boxMax.Y || Math.Max(from.Y, to.Y) < boxMin.Y) {
        return false;
      }

      // Test if the line is completely contained in the viewbox
      if (Math.Min(from.X, to.X) > boxMin.X && Math.Max(from.X, to.X) < boxMax.X &&
        Math.Min(from.Y, to.Y) > boxMin.Y && Math.Max(from.Y, to.Y) < boxMax.Y) {
        return true;
      }

      ClampDestToBox(ref from, ref to, boxMin, boxMax);
      ClampDestToBox(ref to, ref from, boxMin, boxMax);

      return true;
    }

    /// <summary>
    /// Clamps the destination point to the box limits passed in
    /// </summary>
    /// <param name="from"></param>
    /// <param name="dest"></param>
    /// <param name="boxMin"></param>
    /// <param name="boxMax"></param>
    private void ClampDestToBox(ref Vector3F from, ref Vector3F dest, Vector2F boxMin, Vector2F boxMax) {
      Vector3F diff = dest - from;

      if (dest.X > boxMax.X) {
        dest = (boxMax.X - from.X) / diff.X * diff + from;
      } else if (dest.X < boxMin.X) {
        dest = (boxMin.X - from.X) / diff.X * diff + from;
      }
      if (dest.Y > boxMax.Y) {
        dest = (boxMax.Y - from.Y) / diff.Y * diff + from;
      } else if (dest.Y < boxMin.Y) {
        dest = (boxMin.Y - from.Y) / diff.Y * diff + from;
      }
    }

    /// <summary>
    /// Updates the vertex and index buffer for rendering the line for the current viewpoint
    /// </summary>
    /// <param name="deviceWidth"></param>
    /// <param name="deviceHeight"></param>
    /// <param name="projectionTransform"></param>
    private void UpdateLineGeometry(int deviceWidth, int deviceHeight, Matrix4x4D projectionTransform) {
      this.Indices.Clear();
      this.Vertices.Clear();

      // Get all the points in screen coordinates
      IList<Vector3F> list = this.TransformedPoints(projectionTransform);
      IList<float> widths = this.GetWidths(list);
      bool isDashed = this.ConvertToDashed(ref list, ref widths, deviceWidth, deviceHeight);
      for (int i = 0; i < (list.Count - 1); i++) {
        Vector3F width1 = this.GetWidth(widths[i], deviceWidth, deviceHeight);
        Vector3F width2 = this.GetWidth(widths[i+1], deviceWidth, deviceHeight);
       
        Vector3F difference = list[i + 1] - list[i];
        Vector2F direction = Vector2F.Normalize(new Vector2F(difference.X, difference.Y));
        Vector3F perpendicular1 = new Vector3F(-direction.Y * width1.X, direction.X * width1.Y, 0f);
        Vector3F perpendicular2 = new Vector3F(-direction.Y * width2.X, direction.X * width2.Y, 0f);
        this.Vertices.AddData(new Vertex.PositionOnly(list[i] + perpendicular1));
        this.Vertices.AddData(new Vertex.PositionOnly(list[i] - perpendicular1));
        this.Vertices.AddData(new Vertex.PositionOnly(list[i+1] + perpendicular2));
        this.Vertices.AddData(new Vertex.PositionOnly(list[i+1] - perpendicular2));
        int num2 = i*2;
        ushort[] data = new ushort[] { (ushort)num2, (ushort)(num2 + 2), (ushort)(num2 + 3), (ushort)num2, (ushort)(num2 + 3), (ushort)(num2 + 1) };
        this.Indices.AddDataUnsafe(data, PrimitiveType.TriangleList, 0);
        if (isDashed)
          i++;
      }
    }

    /// <summary>
    /// Creates a list of widths for the points passed in.
    /// </summary>
    /// <param name="list"></param>
    /// <returns></returns>
    private IList<float> GetWidths(IList<Vector3F> list) {
      float[] progressiveLength = new float[list.Count];
      float totalLength = 0;
      for (int i = 0; i < (list.Count-1); ++i) {
        progressiveLength[i] = totalLength;
        totalLength += Vector3F.Distance(list[i], list[i + 1]);
      }
      progressiveLength[list.Count - 1] = totalLength;

      if (totalLength > 0) {
        float widthDiff = lineEndWidth - lineStartWidth;
        float invTotalLength = 1f/totalLength;
        for (int i = 0; i < progressiveLength.Length; ++i) {
          progressiveLength[i] = lineStartWidth + widthDiff * progressiveLength[i] * invTotalLength;
        }
      }

      return progressiveLength;
    }

    /// <summary>
    /// Updates the colors, depth buffer, and alpha settings for rendering the mesh
    /// </summary>
    private void UpdateRenderState() {
      RenderState renderState = base.State.RenderState;
      if (this.alphaEnable) {
        renderState.Alpha.Enabled = true;
        renderState.Alpha.SourceBlend = Blend.SourceAlpha;
        renderState.Alpha.DestinationBlend = Blend.InvSourceAlpha;
      } else {
        renderState.Alpha.Enabled = false;
      }
      renderState.Lighting.Enabled = false;
      renderState.Cull.Enabled = false;
      renderState.Stages[0].Blending.AlphaArgument1 = TextureArgument.TFactor;
      renderState.Stages[0].Blending.AlphaOperation = TextureOperation.SelectArg1;
      renderState.Stages[0].Blending.ColorArgument1 = TextureArgument.TFactor;
      renderState.Stages[0].Blending.ColorOperation = TextureOperation.SelectArg1;
      renderState.Stages[1].Enabled = false;
      renderState.TextureFactor.Value = this.color;
      renderState.ZBuffer.Enabled = zBufferEnable;
      renderState.ZBuffer.ZBufferFunction = Compare.Less;
    }

    #endregion Private Methods

    // ************************************************************************
    // Properties
    // ************************************************************************
    #region Properties

    /// <summary>
    /// A list of points that make up the line to be drawn. Points are in
    /// Lat, Lon, Alt format.
    /// </summary>
    public List<Vector3F> Points {
      get {
        return points;
      }
    }

    /// <summary>
    /// Gets and Sets the starting width of the line in pixels
    /// </summary>
    public float LineStartWidth {
      get {
        return this.lineStartWidth;
      }
      set {
        if (value <= 0f) {
          throw new ArgumentOutOfRangeException("value");
        }
        this.lineStartWidth = value;
      }
    }

    /// <summary>
    /// Gets and Sets the end width of the line in pixels
    /// </summary>
    public float LineEndWidth {
      get {
        return this.lineEndWidth;
      }
      set {
        if (value <= 0f) {
          throw new ArgumentOutOfRangeException("value");
        }
        this.lineEndWidth = value;
      }
    }

    /// <summary>
    /// Gets and Sets the length of the dashes in pixels
    /// </summary>
    public float DashLength {
      get {
        return dashLength;
      }
      set {
        if (value < 0f) {
          throw new ArgumentOutOfRangeException("value");
        }
        dashLength = value;
      }
    }

    /// <summary>
    /// Gets and Sets the length of the gaps between the dashes in pixels
    /// </summary>
    public float GapLength {
      get {
        return gapLength;
      }
      set {
        if (value < 0f) {
          throw new ArgumentOutOfRangeException("value");
        }
        gapLength = value;
      }
    }

    /// <summary>
    /// Gets and Sets the offset of the dashes in pixels.
    /// </summary>
    /// <remarks>
    ///  Incrementing this by .5 every frame gives a moving ant effect.
    /// </remarks>
    public float DashOffset {
      get {
        return dashOffset;
      }
      set {
        if (value < 0f) {
          throw new ArgumentOutOfRangeException("value");
        }
        dashOffset = value;
      }
    }

    /// <summary>
    /// Gets and Sets if alpha blending is enabled for this line
    /// </summary>
    public bool AlphaEnable {
      get {
        return this.alphaEnable;
      }
      set {
        if (this.alphaEnable != value) {
          this.alphaEnable = value;
          this.UpdateRenderState();
        }
      }
    }

    /// <summary>
    /// Gets and Sets if the Depth Buffer is enabled for this mesh
    /// </summary>
    public bool ZBufferEnable {
      get {
        return this.zBufferEnable;
      }
      set {
        if (this.zBufferEnable != value) {
          this.zBufferEnable = value;
          this.UpdateRenderState();
        }
      }
    }

    /// <summary>
    /// Gets and Sets the color to draw this line in
    /// </summary>
    public System.Drawing.Color Color {
      get {
        return this.color;
      }
      set {
        if (this.color != value) {
          this.color = value;
          this.UpdateRenderState();
        }
      }
    }

    #endregion Properties
  }
}

A co-worker had a problem today and asked for my help. He was using the Chart control from the WPF Toolkit and wanted to render it to a bitmap and save it to a file without actually adding the Chart to the visual tree. There are plenty of samples on how to render a WPF Visual to a bitmap, including This One that I wrote 4 years ago, but they don’t seem to render if the Visual isn’t a child Visual in a WPF visual tree. The reason for this failure to render is because the Visual isn’t told what size container it has to fill and so the ActualWidth and ActualHeight properties are zero.

So how to get around this? We need to perform the same setup that a panel does on the Visual. If we look at the MSDN Documentation For Custom Panels we will find this code:

public class PlotPanel : Panel
{
    // Default public constructor
    public PlotPanel()
        : base()
    {
    }

    // Override the default Measure method of Panel
    protected override Size MeasureOverride(Size availableSize)
    {
        Size panelDesiredSize = new Size();

        // In our example, we just have one child.
        // Report that our panel requires just the size of its only child.
        foreach (UIElement child in InternalChildren)
        {
            child.Measure(availableSize);
            panelDesiredSize = child.DesiredSize;
        }

        return panelDesiredSize ;
    }
    protected override Size ArrangeOverride(Size finalSize)
    {
        foreach (UIElement child in InternalChildren)
        {
            double x = 50;
            double y = 50;

            child.Arrange(new Rect(new Point(x, y), child.DesiredSize));
        }
        return finalSize; // Returns the final Arranged size
    }
}

Looking at the code above we can see that the methods Measure(Size size) and Arrange(Size size) are called on the child Visuals. Here is some sample code inside an event handler to show how it is done:

private void button1_Click(object sender, RoutedEventArgs e) {

  Button button = new Button();
  button.Content = "Parentless Button";
  double width = 200;
  double height = 200;
  button.Width = width;
  button.Height = height;

  // Cure for getting an ActualWidth and ActualHeight of 0
  button.Measure(new Size(Double.PositiveInfinity,
                         Double.PositiveInfinity));
  button.Arrange(new Rect(button.DesiredSize));

  RenderTargetBitmap rtb = new RenderTargetBitmap((int)width, (int)height, 96.0, 96.0, PixelFormats.Pbgra32);
  rtb.Render(button);
  this.image1.Source = rtb;
}

You will also find a link at the top to the source code for a sample application that demonstrates rendering an un-parented Button to a Bitmap, and then displaying that Bitmap. When you run the application and click the “Test” button it will look like this:

The “Parentless Button” isn’t a button, but just a Bitmap displayed in an image control.

This is the 3rd part of my posts on Bing Maps 3D. The 3 parts are:

- Adding Arbitrary Shapes. Creating a Cube mesh and a Sphere mesh.
- Rendering shapes past the Far Clip. How to render shapes whose coordinates are outside the View Port.
- Clipping shapes with an imposter Earth.

In the previous Bing Maps 3D post we had overcome the problem of 3D objects being clipped by the far plane on the horizon, but now we had the problem that they weren’t being clipped by the earth as seen for the yellow cube in the picture below:

The solution to this is to create another sphere to imitate the earth, but smaller and closer to the camera. The sphere should be scaled down by the same factor that the cube is scaled down when it is brought closer to the camera. The sphere rendered for the Earth in Bing Maps 3D isn’t an exact sphere, it is a little bit squashed, so to make sure we get the same shape we measure the earth from within Bing Maps 3D and apply those measurements to our sphere. Here’s where we use the SetToReplaceEarth() method that I included in the Sphere Mesh code in the first post which looks like this:

/// <summary>
/// Scales and positions this sphere such that it replaces the earth
/// </summary>
public void SetToReplaceEarth() {
  Microsoft.MapPoint.Rendering3D.Cameras.GeodeticViewpoint view = new Microsoft.MapPoint.Rendering3D.Cameras.GeodeticViewpoint();
  Vector3D[] vecs = new Vector3D[6];
  view.Position.Location = LatLonAlt.CreateUsingDegrees(0, 0, 0);
  vecs[0] = view.Position.Vector;

  view.Position.Location = LatLonAlt.CreateUsingDegrees(0,180, 0);
  vecs[1] = view.Position.Vector;

  view.Position.Location = LatLonAlt.CreateUsingDegrees(-90, 0, 0);
  vecs[2] = view.Position.Vector;

  view.Position.Location = LatLonAlt.CreateUsingDegrees(90, 0, 0);
  vecs[3] = view.Position.Vector;

  view.Position.Location = LatLonAlt.CreateUsingDegrees(0, 90, 0);
  vecs[4] = view.Position.Vector;

  view.Position.Location = LatLonAlt.CreateUsingDegrees(0, -90, 0);
  vecs[5] = view.Position.Vector;

  Vector3D average = new Vector3D(0, 0, 0);
  Vector3D min = new Vector3D(double.MaxValue, double.MaxValue, double.MaxValue);
  Vector3D max = new Vector3D(double.MinValue, double.MinValue, double.MinValue);

  for (int i = 0; i < vecs.Length; ++i) {
    average.Add(vecs[i]);
    min = min.Min(vecs[i]);
    max = max.Max(vecs[i]);
  }
  average.MultiplyBy(1.0 / 6.0);

  this.Position = average;

  Vector3D diff = max.Subtract(min);
  diff.MultiplyBy(0.5 / radius);
  this.Scale = diff;

}

We put this sphere into the scene at 0,0,0, and to make it visible I will give it a very transparent yellow color. Here’s the result:

So it’s working now right? Well not quite, you see I’ve drawn the sphere above to have 20 stacks, 20 slices, so it’s not a very high resolution mesh, so when you get lower the clip is quite a bit under the surface. In the image below the yellow line should stop at the same place as the red cube does:

Easy fix right? Lets change the imitation Earth to a 200×200 mesh:

What happened? Well the old DirectX cards could only handle a 16 bit Index Buffer , which is what we set up our sphere mesh as when we declared it like this:

public class SphereMesh : MeshGraphicsObject<Vertex.PositionNormal, ushort>{
}

Fortunately for me my video card can handle a 32 bit index buffer, so I change the code throughout to have uint instead of ushort and now it works fine and the error in the clip at low altitude is improved:

This wraps up my 3 part series on adding 3D Meshes to Bing Maps 3D. Enjoy playing with the sample application attached at the top of all 3 articles.

This is the 2nd part of my posts on Bing Maps 3D. The 3 parts are:

- Adding Arbitrary Shapes. Creating a Cube mesh and a Sphere mesh.
- Rendering shapes past the Far Clip. How to render shapes whose coordinates are outside the View Port.
- Clipping shapes with an imposter Earth.

In the previous part we discovered a problem with rendering large shapes on the Earth in Bing Maps that the far clip plane is at the horizon, and so anything drawn past that is clipped, as is seen happening to far corner of the red cube in the screenshot below:

To fix this we need to move the cube towards the camera, and shrink it so that the whole cube is in front of the far clip plane, but looks the same. Here’s a diagram showing what we need to do:

The red square in the diagram shows the normal cube being clipped by the far clip plane, and the yellow square shows the fixed version. The blue circle is the camera, and the blue lines are the rays cast from the camera. As you can see from the cameras point of view the yellow square looks to be in the same position as the red square, but it is closer to the camera.

To implement this I wrote a method in my Bing Maps 3D Actor called MoveTowardCamera. This method sets up a post render transform to apply to the cube to move it towards the camera. Strangely the post transform in Bing Maps gets applied before the object is positioned, as I need to know where the object is in order to move it toward the camera, I first translate it to its set position, scale it down, then move it towards the camera in proportion to it’s scale. Here’s the code to do this:

public static void MoveTowardCamera(WorldGraphicsObject mesh, double scale, SceneState sceneState)
{
  CameraData camera = null;

  if (sceneState.TryGetData<CameraData>(out camera))
  {
    // Calculate a transform that moves the object closer
    // to the camera, but scales it down so it appears
    // to be the same size

    Vector3D cameraPos = camera.Snapshot.Position.Vector;
    Vector3D meshPos = mesh.Position;

    // Transform works about the centre of an object, so to work properly
    // need to move to position, apply scaling towards camera, and then move back

    Matrix4x4D postTransform = Matrix4x4D.Translation(meshPos);
    postTransform.Multiply(Matrix4x4D.Scaling(scale, scale, scale));
    postTransform.Multiply(Matrix4x4D.Translation(cameraPos * (1 - scale)));
    postTransform.Multiply(Matrix4x4D.Translation(-meshPos));

    mesh.PostTransform = postTransform;
  }
}

Here’s how it looks in 3D:

The red cube is being clipped by the far clip plane, and the yellow cube is closer to the camera but overlays the red cube correctly. The problem now is that the yellow frame isn’t being clipped by the Earth as it is much closer than the Earth. In my next post I will show how to get the yellow cube to be clipped by the Earth.

The function to find the end of the last complete character (and thus how many extra bytes are left over) is cntExtraBytes and it takes the last three bytes as parameters to check as that’s the maximum number of bytes that could possible be a piece of a character. The rules are derived from the wikipedia’s description of the UTF-8 standard of which bytes start two, three and four byes sequences. I’ve put in a few ASSERTs to help clarify what the possible values are of each byte.

using System;
using System.Linq;
using System.Text;
using System.IO;

namespace UTF8_to_Unicode
{
    class Program
    {
        // UTF-8 is encoded in up to 4 bytes, need to make sure we don't chop a sequence in half
        // 240+: First in a 4 byte sequence
        // 224+: 3 byte sequence
        // 194+: 2 byte sequence
        // 128+: Part of a multi-byte sequence
        // 128-: Not part of a multi-byte sequence
        private static int cntExtraBytes(byte c3, byte c2, byte c1)
        {
            if (c1 < 128) return 0;

            if (c1 >= 192) return 1;

            System.Diagnostics.Debug.Assert(c1 >= 128);
            System.Diagnostics.Debug.Assert(c2 >= 128);
            System.Diagnostics.Debug.Assert(c3 >= 128);

            System.Diagnostics.Debug.Assert( 192 >= c1 && c1 <= 224);

            if (c2 >= 224) return 2;
            if (c2 >= 192) return 0;

            if (c3 >= 240) return 3;
            if (c3 >= 224) return 0;

            System.Diagnostics.Debug.Assert(c3 <= 192);

            return 0;
        }

        static void Main(string[] args)
        {
            System.Diagnostics.Stopwatch sw = new System.Diagnostics.Stopwatch();
            System.Diagnostics.Stopwatch swRead = new System.Diagnostics.Stopwatch();
            System.Diagnostics.Stopwatch swEncode = new System.Diagnostics.Stopwatch();
            System.Diagnostics.Stopwatch swWrite = new System.Diagnostics.Stopwatch();
            sw.Start();

            if (args.Count() != 3)
            {
                Console.Out.WriteLine("Usage: buffer_size(kb) UTF8_in Unicode_out");
                Console.Out.WriteLine("Doesn't work for files in excess of 2GB");
                System.Environment.Exit(0);
            }

            FileStream fs = new FileStream(args[0], FileMode.Open, FileAccess.Read, FileShare.None, 1, true);
            FileStream fs_out = new FileStream(args[1], FileMode.CreateNew, FileAccess.Write, FileShare.None,1,true);

            int total = (int)fs.Length;
            int cBytesOut = 0, cBytesIn = 0, cExtraBytes = 0, cLength = 0;
            int piece = Int32.Parse(args[2]) * 1024;
            Console.Out.Write(piece.ToString() + ",");

            byte[] bytes_in = new byte[piece+3];
            byte[] bytes_out = null;

            do
            {
                swRead.Start();

                cLength = fs.Read(bytes_in, cExtraBytes, piece) + cExtraBytes;

                cExtraBytes = cntExtraBytes(bytes_in[cLength - 3], bytes_in[cLength - 2], bytes_in[cLength - 1]);

                swRead.Stop();

                swEncode.Start();
                bytes_out = Encoding.Convert(Encoding.UTF8, Encoding.Unicode, bytes_in, 0, cLength - cExtraBytes);
                swEncode.Stop();

                if (cExtraBytes > 0)
                {
                    Buffer.BlockCopy(bytes_in, cLength - cExtraBytes, bytes_in, 0, cExtraBytes);
                }

                swWrite.Start();
                fs_out.Write(bytes_out, 0, bytes_out.Count());
                swWrite.Stop();

                cBytesOut += bytes_out.Count();
                cBytesIn += piece;
            }
            while (total > cBytesIn);

            fs.Close();
            fs_out.Close();

            sw.Stop();

            Console.Out.Write(total.ToString() + "," + cBytesOut.ToString() + ",");

            Console.Out.Write((sw.ElapsedMilliseconds / 1000.0).ToString() + ",");
            Console.Out.Write((swRead.ElapsedMilliseconds / 1000.0).ToString() + ",");
            Console.Out.Write((swEncode.ElapsedMilliseconds / 1000.0).ToString() + ",");
            Console.Out.WriteLine((swWrite.ElapsedMilliseconds / 1000.0).ToString());
        }
    }
}

As you can see there’s a lot of timing code in there and when run it shows that the reading and writing takes a suspiciously low amount of time, exceeding the throughput of the disks which tells us that it benefiting from both read and write caching.

Interestingly caching makes it rather difficult to accurately measure the performance of your disks but you can use the free version of HDTune or if you want something with source code then diskspd from Microsoft Research. Both return comparable results but HDTune is more comprehensive and provides pretty graphics.

The program takes three parameters, the first is buffer size in kb, and changing that parameter can greatly affect the performance with larger buffers performing better. A run with a 256kb buffer size takes only 0.1 secs to read a 6MB file, convert it and write out the resultant 10MB file:

C:\data>UTF8_to_Unicode.exe 256 input.u8 output.u16
262144,6152243,10391224,0.1,0.011,0.051,0.035

So there it is, fast enough, runnable from a build script and uses only a fixed amount of memory.

Warning: There is no error handling or checking of input parameters of any kind. You will need to add those in order to use this code.

Recently I have been implementing a 3D Viewer application based on Bing Maps 3D which you can download from Here. From this work I have generated enough material for a few posts, they will be as follows:

- Adding Arbitrary Shapes. Creating a Cube mesh and a Sphere mesh.
- Rendering shapes past the Far Clip. How to render shapes whose coordinates are outside the View Port.
- Clipping shapes with an imposter Earth.

The work I’ve done started off with 2 samples that I got from Here

Here is a screenshot of a translucent cube rendered onto the surface of the Earth:

The cube is drawn by a plug-in that implements an “Actor” and is implemented as a MeshGraphicsObject. Here is the source code for the Cube Mesh which is fairly straight forward:

using System;
using System.Collections.Generic;
using System.Text;
using Microsoft.MapPoint.Rendering3D.GraphicsProxy;
using Microsoft.MapPoint.Geometry.VectorMath;
using Microsoft.MapPoint.Geometry;

namespace _3DViewerBingMaps.Plugin {

  /// <summary>
  /// This class is a Mesh for rendering a Cube
  /// </summary>
  public class CubeMesh : MeshGraphicsObject<Vertex.PositionNormal, ushort> {

    #region Fields

    /// <summary>
    /// The material id for the material attached to this mesh
    /// </summary>
    private int materialId;

    #endregion Private

    #region Methods

    /// <summary>
    /// Constructor that takes a material
    /// </summary>
    /// <param name="radius"></param>
    /// <param name="material"></param>
    public CubeMesh(float radius, Material material) : this(radius, material, -1) { }

    /// <summary>
    /// Constructor that takes a material id
    /// </summary>
    /// <param name="radius"></param>
    /// <param name="materialId"></param>
    public CubeMesh(float radius, int materialId) : this(radius, null, materialId) { }

    /// <summary>
    /// Constructor that takes a material and a material id. The material
    /// id is ignored if the material isn't equal to null.
    /// </summary>
    /// <param name="radius"></param>
    /// <param name="material"></param>
    /// <param name="materialId"></param>
    public CubeMesh(float radius, Material material, int materialId)
      : base(GraphicsBufferUsage.Static, GraphicsBufferUsage.Static, 8, 36, true) {

      // Get the vertices and indices
      List<Vector3F> vertices = new List<Vector3F>();
      List<ushort> indices = new List<ushort>();
      CalculateVerticesIndicies(vertices, indices, radius);

      // Add the vertices to the mesh
      foreach (Vector3F vector in vertices) {
        this.Vertices.AddData(new Vertex.PositionNormal(vector.X, vector.Y, vector.Z, vector.X, vector.Y, vector.Z));
      }

      // Set up the material
      if (material != null) {
        this.materialId = this.Materials.Add(material);
      }
      else {
        this.materialId = materialId;
      }

      // Add the indices to the mesh as a list of triangles
      this.Indices.AddData(indices.ToArray(), PrimitiveType.TriangleList, this.materialId);

    }

    /// <summary>
    /// Calculates the vertices and the indices for the mesh
    /// </summary>
    /// <param name="vertices"></param>
    /// <param name="indices"></param>
    /// <param name="radius"></param>
    protected virtual void CalculateVerticesIndicies(List<Vector3F> vertices, List<ushort> indices, float radius) {

      // Add all the corners of the cube
      vertices.Add(new Vector3F(-radius, -radius, -radius));
      vertices.Add(new Vector3F(-radius, -radius, radius));
      vertices.Add(new Vector3F(-radius, radius, -radius));
      vertices.Add(new Vector3F(-radius, radius, radius));
      vertices.Add(new Vector3F(radius, -radius, -radius));
      vertices.Add(new Vector3F(radius, -radius, radius));
      vertices.Add(new Vector3F(radius, radius, -radius));
      vertices.Add(new Vector3F(radius, radius, radius));

      // Build a triangle list for the above vertices
      indices.AddRange(new ushort[] {
        0, 2, 1,
        1, 2, 3,
        1, 3, 5,
        5, 3, 7,
        5, 7, 4,
        4, 7, 6,
        4, 6, 0,
        0, 6, 2,
        6, 7, 2,
        2, 7, 3,
        5, 4, 1,
        1, 4, 0
      });
    }

    #endregion Methods

    #region Properties

    /// <summary>
    /// Gets the material id
    /// </summary>
    public int MaterialId {
      get {
        return this.materialId;
      }
    }

    #endregion Properties
  }
}

 

That isn’t terribly exciting. What about a Sphere that looks like this:

The code for the Sphere Mesh is a bit more interesting:

using System;
using System.Collections.Generic;
using System.Text;
using Microsoft.MapPoint.Rendering3D.GraphicsProxy;
using Microsoft.MapPoint.Geometry.VectorMath;
using Microsoft.MapPoint.Geometry;
using Microsoft.MapPoint.Rendering3D;

namespace _3DViewerBingMaps.Plugin {

  /// <summary>
  /// This class is a mesh for rendering a sphere
  /// </summary>
  public class SphereMesh : MeshGraphicsObject<Vertex.PositionNormal, ushort> {

    #region Fields

    /// <summary>
    /// The material id for the material attached to this mesh
    /// </summary>
    private int materialId;
    /// <summary>
    /// The radius of the sphere
    /// </summary>
    private float radius;

    #endregion Fields

    #region Methods

    /// <summary>
    /// Constructor that takes a material
    /// </summary>
    /// <param name="radius"></param>
    /// <param name="n"></param>
    /// <param name="material"></param>
    public SphereMesh(float radius, int n, Material material) : this(radius, n, material, -1) { }

    /// <summary>
    /// Constructor that takes a material id
    /// </summary>
    /// <param name="radius"></param>
    /// <param name="n"></param>
    /// <param name="materialId"></param>
    public SphereMesh(float radius, int n, int materialId) : this(radius, n, null, materialId) { }

    /// <summary>
    /// Constructor that takes a material and a material id. The material id is ignored
    /// if the material isn't null
    /// </summary>
    /// <param name="radius"></param>
    /// <param name="n"></param>
    /// <param name="material"></param>
    /// <param name="materialId"></param>
    protected SphereMesh(float radius, int n, Material material, int materialId)
      : base(GraphicsBufferUsage.Static, GraphicsBufferUsage.Static, 4 * (2 * n * n + 3 * n + 1), 24 * (2 * n * n + 3 * n + 1), true) {

      this.radius = radius;

      // Get the indices and vertices
      List<Vector3F> vertices = new List<Vector3F>();
      List<ushort> indices = new List<ushort>();
      CalculateVerticesIndicies(vertices, indices, radius, n);

      // Add the vertices to the mesh
      foreach (Vector3F vector in vertices) {
        this.Vertices.AddData(new Vertex.PositionNormal(vector.X, vector.Y, vector.Z, vector.X, vector.Y, vector.Z));
      }

      // Set the material
      if (material != null) {
        this.materialId = this.Materials.Add(material);
      }
      else {
        this.materialId = materialId;
      }

      // Add the indices to the mesh as a list of triangles
      this.Indices.AddData(indices.ToArray(), PrimitiveType.TriangleList, this.materialId);

    }

    /// <summary>
    /// Calculates the vertices and indices for rendering this mesh
    /// </summary>
    /// <param name="vertices"></param>
    /// <param name="indices"></param>
    /// <param name="radius"></param>
    /// <param name="n"></param>
    private static void CalculateVerticesIndicies(List<Vector3F> vertices, List<ushort> indices, float radius, int n) {
      double segmentRad = Math.PI / 2 / (n + 1);
      int numberOfSeparators = 4 * n + 4;

      for (int segmentNo = -n; segmentNo <= n; segmentNo++) {
        double r_e = radius * Math.Cos(segmentRad * segmentNo);
        double y_e = radius * Math.Sin(segmentRad * segmentNo);

        for (int s = 0; s <= (numberOfSeparators - 1); s++) {
          double z_s = r_e * Math.Sin(segmentRad * s) * (-1);
          double x_s = r_e * Math.Cos(segmentRad * s);
          vertices.Add(new Vector3F((float)x_s, (float)y_e, (float)z_s));
        }
      }
      vertices.Add(new Vector3F(0, radius, 0));
      vertices.Add(new Vector3F(0, -1 * radius, 0));

      for (int segmentNo = 0; segmentNo < 2 * n; segmentNo++) {
        for (int i = 0; i < numberOfSeparators; i++) {
          indices.Add((ushort)(segmentNo * numberOfSeparators + i));
          indices.Add((ushort)(segmentNo * numberOfSeparators + i +
                              numberOfSeparators));
          indices.Add((ushort)(segmentNo * numberOfSeparators + (i + 1) %
                              numberOfSeparators + numberOfSeparators));

          indices.Add((ushort)(segmentNo * numberOfSeparators + (i + 1) %
                              numberOfSeparators + numberOfSeparators));
          indices.Add((ushort)(segmentNo * numberOfSeparators +
                             (i + 1) % numberOfSeparators));
          indices.Add((ushort)(segmentNo * numberOfSeparators + i));
        }
      }

      for (int i = 0; i < numberOfSeparators; i++) {
        indices.Add((ushort)((2 * n) * numberOfSeparators + i));
        indices.Add((ushort)(numberOfSeparators * (2 * n + 1)));
        indices.Add((ushort)((2 * n) * numberOfSeparators + (i + 1) %
                            numberOfSeparators));
      }

      for (int i = 0; i < numberOfSeparators; i++) {
        indices.Add((ushort)(i));
        indices.Add((ushort)((i + 1) % numberOfSeparators));
        indices.Add((ushort)(numberOfSeparators * (2 * n + 1) + 1));
      }
    }

    /// <summary>
    /// Scales and positions this sphere such that it replaces the earth
    /// </summary>
    public void SetToReplaceEarth() {
      Microsoft.MapPoint.Rendering3D.Cameras.GeodeticViewpoint view = new Microsoft.MapPoint.Rendering3D.Cameras.GeodeticViewpoint();
      Vector3D[] vecs = new Vector3D[6];
      view.Position.Location = LatLonAlt.CreateUsingDegrees(0, 0, 0);
      vecs[0] = view.Position.Vector;

      view.Position.Location = LatLonAlt.CreateUsingDegrees(0,180, 0);
      vecs[1] = view.Position.Vector;

      view.Position.Location = LatLonAlt.CreateUsingDegrees(-90, 0, 0);
      vecs[2] = view.Position.Vector;

      view.Position.Location = LatLonAlt.CreateUsingDegrees(90, 0, 0);
      vecs[3] = view.Position.Vector;

      view.Position.Location = LatLonAlt.CreateUsingDegrees(0, 90, 0);
      vecs[4] = view.Position.Vector;

      view.Position.Location = LatLonAlt.CreateUsingDegrees(0, -90, 0);
      vecs[5] = view.Position.Vector;

      Vector3D average = new Vector3D(0, 0, 0);
      Vector3D min = new Vector3D(double.MaxValue, double.MaxValue, double.MaxValue);
      Vector3D max = new Vector3D(double.MinValue, double.MinValue, double.MinValue);

      for (int i = 0; i < vecs.Length; ++i) {
        average.Add(vecs[i]);
        min = min.Min(vecs[i]);
        max = max.Max(vecs[i]);
      }
      average.MultiplyBy(1.0 / 6.0);

      this.Position = average;

      Vector3D diff = max.Subtract(min);
      diff.MultiplyBy(0.5 / radius);
      this.Scale = diff;

    }

    #endregion Methods

    #region Properties

    /// <summary>
    /// Gets the material id for the material attached to this mesh
    /// </summary>
    public int MaterialId {
      get {
        return materialId;
      }
    }

    #endregion Properties
  }
}

 

You might notice that there is a method in the Sphere mesh class called SetToImitateEarth(). I will use that in 2 posts time.

One of the problems with these shapes though is that as we get closer to the ground the Far-Clip-Plane in Bing Maps 3D gets closer and cuts off the back of our shape as seen in this screenshot:

In my next Bing Maps 3D article I will show how to get around this limitation.

This method of distribution is quite ad-hoc so when a client rings me up and asks me to view their screen to look at something, it’s useful to know what build they are running. To facillitate this I print the link date in the main Window Title so I instantly have an idea about how old the version is that I am looking at. This date is calculated at run time. To do this requires reading in the Portable Executable (PE) header from the EXE or DLL file that is responsible for updating the Window. As many of the apps I write are in C# I’ve written a general purpose PE Header reading class in C# which has a utility method for getting the link date of the Calling Assembly.

To get the date simply do this:

DateTime linkDate = PeHeaderReader.GetCallingAssemblyHeader().TimeStamp;

Here’s is the complete code:

using System;
using System.Collections.Generic;
using System.Runtime.InteropServices;
using System.IO;

namespace _3DViewerControls.Data {

  // Reads in the header information of the Portable Executable format.
  // Provides information such as the date the assembly was compiled.
  public class PeHeaderReader {
    #region File Header Structures

    public struct IMAGE_DOS_HEADER {      // DOS .EXE header
      public UInt16 e_magic;              // Magic number
      public UInt16 e_cblp;               // Bytes on last page of file
      public UInt16 e_cp;                 // Pages in file
      public UInt16 e_crlc;               // Relocations
      public UInt16 e_cparhdr;            // Size of header in paragraphs
      public UInt16 e_minalloc;           // Minimum extra paragraphs needed
      public UInt16 e_maxalloc;           // Maximum extra paragraphs needed
      public UInt16 e_ss;                 // Initial (relative) SS value
      public UInt16 e_sp;                 // Initial SP value
      public UInt16 e_csum;               // Checksum
      public UInt16 e_ip;                 // Initial IP value
      public UInt16 e_cs;                 // Initial (relative) CS value
      public UInt16 e_lfarlc;             // File address of relocation table
      public UInt16 e_ovno;               // Overlay number
      public UInt16 e_res_0;              // Reserved words
      public UInt16 e_res_1;              // Reserved words
      public UInt16 e_res_2;              // Reserved words
      public UInt16 e_res_3;              // Reserved words
      public UInt16 e_oemid;              // OEM identifier (for e_oeminfo)
      public UInt16 e_oeminfo;            // OEM information; e_oemid specific
      public UInt16 e_res2_0;             // Reserved words
      public UInt16 e_res2_1;             // Reserved words
      public UInt16 e_res2_2;             // Reserved words
      public UInt16 e_res2_3;             // Reserved words
      public UInt16 e_res2_4;             // Reserved words
      public UInt16 e_res2_5;             // Reserved words
      public UInt16 e_res2_6;             // Reserved words
      public UInt16 e_res2_7;             // Reserved words
      public UInt16 e_res2_8;             // Reserved words
      public UInt16 e_res2_9;             // Reserved words
      public UInt32 e_lfanew;             // File address of new exe header
    }

    [StructLayout(LayoutKind.Sequential, Pack = 1)]
    public struct IMAGE_OPTIONAL_HEADER32 {
      public UInt16 Magic;
      public Byte MajorLinkerVersion;
      public Byte MinorLinkerVersion;
      public UInt32 SizeOfCode;
      public UInt32 SizeOfInitializedData;
      public UInt32 SizeOfUninitializedData;
      public UInt32 AddressOfEntryPoint;
      public UInt32 BaseOfCode;
      public UInt32 BaseOfData;
      public UInt32 ImageBase;
      public UInt32 SectionAlignment;
      public UInt32 FileAlignment;
      public UInt16 MajorOperatingSystemVersion;
      public UInt16 MinorOperatingSystemVersion;
      public UInt16 MajorImageVersion;
      public UInt16 MinorImageVersion;
      public UInt16 MajorSubsystemVersion;
      public UInt16 MinorSubsystemVersion;
      public UInt32 Win32VersionValue;
      public UInt32 SizeOfImage;
      public UInt32 SizeOfHeaders;
      public UInt32 CheckSum;
      public UInt16 Subsystem;
      public UInt16 DllCharacteristics;
      public UInt32 SizeOfStackReserve;
      public UInt32 SizeOfStackCommit;
      public UInt32 SizeOfHeapReserve;
      public UInt32 SizeOfHeapCommit;
      public UInt32 LoaderFlags;
      public UInt32 NumberOfRvaAndSizes;
    }

    [StructLayout(LayoutKind.Sequential, Pack = 1)]
    public struct IMAGE_OPTIONAL_HEADER64 {
      public UInt16 Magic;
      public Byte MajorLinkerVersion;
      public Byte MinorLinkerVersion;
      public UInt32 SizeOfCode;
      public UInt32 SizeOfInitializedData;
      public UInt32 SizeOfUninitializedData;
      public UInt32 AddressOfEntryPoint;
      public UInt32 BaseOfCode;
      public UInt64 ImageBase;
      public UInt32 SectionAlignment;
      public UInt32 FileAlignment;
      public UInt16 MajorOperatingSystemVersion;
      public UInt16 MinorOperatingSystemVersion;
      public UInt16 MajorImageVersion;
      public UInt16 MinorImageVersion;
      public UInt16 MajorSubsystemVersion;
      public UInt16 MinorSubsystemVersion;
      public UInt32 Win32VersionValue;
      public UInt32 SizeOfImage;
      public UInt32 SizeOfHeaders;
      public UInt32 CheckSum;
      public UInt16 Subsystem;
      public UInt16 DllCharacteristics;
      public UInt64 SizeOfStackReserve;
      public UInt64 SizeOfStackCommit;
      public UInt64 SizeOfHeapReserve;
      public UInt64 SizeOfHeapCommit;
      public UInt32 LoaderFlags;
      public UInt32 NumberOfRvaAndSizes;
    }

    [StructLayout(LayoutKind.Sequential, Pack = 1)]
    public struct IMAGE_FILE_HEADER {
      public UInt16 Machine;
      public UInt16 NumberOfSections;
      public UInt32 TimeDateStamp;
      public UInt32 PointerToSymbolTable;
      public UInt32 NumberOfSymbols;
      public UInt16 SizeOfOptionalHeader;
      public UInt16 Characteristics;
    }

    #endregion File Header Structures

    #region Private Fields

    // The DOS header
    private IMAGE_DOS_HEADER dosHeader;
    // The file header
    private IMAGE_FILE_HEADER fileHeader;
    // Optional 32 bit file header
    private IMAGE_OPTIONAL_HEADER32 optionalHeader32;
    // Optional 64 bit file header
    private IMAGE_OPTIONAL_HEADER64 optionalHeader64;

    #endregion Private Fields

    #region Public Methods

    public PeHeaderReader(string filePath) {
      // Read in the DLL or EXE and get the timestamp
      using (FileStream stream = new FileStream(filePath, System.IO.FileMode.Open, System.IO.FileAccess.Read)) {
        BinaryReader reader = new BinaryReader(stream);
        dosHeader = FromBinaryReader<IMAGE_DOS_HEADER>(reader);

        // Add 4 bytes to the offset
        stream.Seek(dosHeader.e_lfanew, SeekOrigin.Begin);

        UInt32 ntHeadersSignature = reader.ReadUInt32();
        fileHeader = FromBinaryReader<IMAGE_FILE_HEADER>(reader);
        if (this.Is32BitHeader) {
          optionalHeader32 = FromBinaryReader<IMAGE_OPTIONAL_HEADER32>(reader);
        }
        else {
          optionalHeader64 = FromBinaryReader<IMAGE_OPTIONAL_HEADER64>(reader);
        }
      }
    }

    // Gets the header of the .NET assembly that called this function
    public static PeHeaderReader GetCallingAssemblyHeader() {
      string pathCallingAssembly = System.Reflection.Assembly.GetCallingAssembly().Location;

      // Get the path to the calling assembly, which is the path to the
      // DLL or EXE that we want the time of
      string filePath = System.Reflection.Assembly.GetCallingAssembly().Location;

      // Get and return the timestamp
      return new PeHeaderReader(filePath);
    }

    // Reads in a block from a file and converts it to the struct
    // type specified by the template parameter
    public static T FromBinaryReader<T>(BinaryReader reader) {
      // Read in a byte array
      byte[] bytes = reader.ReadBytes(Marshal.SizeOf(typeof(T)));

      // Pin the managed memory while, copy it out the data, then unpin it
      GCHandle handle = GCHandle.Alloc(bytes, GCHandleType.Pinned);
      T theStructure = (T)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(T));
      handle.Free();

      return theStructure;
    }

    #endregion Public Methods

    #region Properties

    // Gets if the file header is 32 bit or not
    public bool Is32BitHeader {
      get {
        UInt16 IMAGE_FILE_32BIT_MACHINE = 0x0100;
        return (IMAGE_FILE_32BIT_MACHINE & FileHeader.Characteristics) == IMAGE_FILE_32BIT_MACHINE;
      }
    }

    // Gets the file header
    public IMAGE_FILE_HEADER FileHeader {
      get {
        return fileHeader;
      }
    }

    // Gets the optional header
    public IMAGE_OPTIONAL_HEADER32 OptionalHeader32 {
      get {
        return optionalHeader32;
      }
    }

    // Gets the optional header
    public IMAGE_OPTIONAL_HEADER64 OptionalHeader64 {
      get {
        return optionalHeader64;
      }
    }

    // Gets the timestamp from the file header
    public DateTime TimeStamp {
      get {
        // Timestamp is a date offset from 1970
        DateTime returnValue = new DateTime(1970, 1, 1, 0, 0, 0);

        // Add in the number of seconds since 1970/1/1
        returnValue = returnValue.AddSeconds(fileHeader.TimeDateStamp);
        // Adjust to local timezone
        returnValue += TimeZone.CurrentTimeZone.GetUtcOffset(returnValue);

        return returnValue;
      }
    }

    #endregion Properties
  }
}

How it works is that it first reads in the old DOS header, at the end of this header is a file offset to the new NT File Header structure. I seek to that position, and read in the NT File Header. From that header I can get the linker time stamp. As this is a general purpose library I also check whether the header is 32 or 64 bit, and read in either the Optional 32 bit Header, or the Optional 64 bit Header, which can then be used however you like. This is how you would get the 32 bit header:

PeHeaderReader reader = new PeHeaderReader("myDllFileLocation");
if (reader.Is32BitHeader)
{
  PeHeaderReader.IMAGE_OPTIONAL_HEADER32 header32 = reader.OptionalHeader32;
}

After completing the STDEV test on Javascript I wondered how fast I could get the C# version to go using some of the new parallel extensions in C# on my dual core non-HT E4500 (so only two threads). Of course the first problem is that the Parallel Extensions CTP has been taken down but fortunately there’s still a way to get PFX via downloading the Reactive Framework. With that done I started out by using Plinq’s simple AsParallel to directly convert the Alglib STDEV implementation to Linq and then add the AsParallel but soon ran into problem with the execution time ballooning out. To eliminate variables I reverted to a simple calculation of the mean. The problem there is that the CPU is so fast that it is difficult to measure so I decided to run the mean over an array of 100M doubles 5 times and take the average. Since I was doing a comparison of various methods anyway I thought I’d thrown in some for vs foreach comparisons and try “hoisting” the length into a local variable. See also: Jon Skeet: For vs Foreach (with hoisting).

The various mean variants I tried are below:

• calc_mean: Basic calculation of the mean with a for loop
• hoist_mean: Take the array.Length out of the loop test into a local variable
• func_mean: What’s the overhead of making a function call, is C# able to inline it?
• f_hoist_mean: Use both a functional call and hoist
• linq_mean: Use Linq
• plinq_mean: Use Parallel Linq
• foreach_mean: Same as calc_mean but with foreach rather than for
• pfor_mean: Use Parallel.For
• pforeach_mean: Use Parallel.ForEach
• part_mean: Use Parallel.invoke with custom problem partitioning (split in half)

Once the mean was done I tried some stdev variants:

• calc_stddev: Implementation for Alglib
• part_stddev: Parallel invoke with custom partition of the data
• part2_stddev: Parallel invoke with custom partition of the algorithm

Lastly I used two different timing techniques:

• Time(Action): Straightforward timing with the Stopwatch (in ms)
• pTime(Action): Timing the runs in parallel

Here’s the C# code:

using System;
using System.Linq;
using System.Text;
using System.Linq.Parallel;
using System.Diagnostics;
using System.Threading;
using System.Threading.Tasks;

using System.Reflection;
using System.Collections.Generic;

namespace ParallelComparison
{
    class Program
    {
        static int array_size = (int)1e8; // 800MiB
        static int loop_count = 5;

        #region Mean

        public static double calc_mean(double[] x)
        {
            double mean = 0.0;

            for (int i = 0; i < x.Length; i++)
            {
                mean += x[i];
            }

            return mean / x.Length;
        }

        public static double hoist_mean(double[] x)
        {
            double mean = 0.0;
            int len = x.Length;

            for (int i = 0; i < len; i++)
            {
                mean += x[i];
            }

            return mean / x.Length;
        }

        private static double add(double a, double b) { return a + b; }

        public static double func_mean(double[] x)
        {
            double mean = 0.0;

            for (int i = 0; i < x.Length; i++)
            {
                mean = add(mean, x[i]);
            }

            return mean / x.Length;
        }

        public static double f_hoist_mean(double[] x)
        {
            double mean = 0.0;
            int len = x.Length;

            for (int i = 0; i < len; i++)
            {
                mean = add(mean, x[i]);
            }

            return mean / x.Length;
        }

        public static double linq_mean(double[] x)
        {
            return x.Average<double>(a => a);
        }

        public static double plinq_mean(double[] x)
        {
            return x.AsParallel().Average<double>(a => a);
        }

        public static double foreach_mean(double[] x)
        {
            double mean = 0.0;

            foreach (double a in x)
            {
                mean += a;
            }
            return mean / x.Length;
        }

        public static double pfor_mean(double[] x)
        {
            double mean = 0.0;

            object syncObj = new object();

            Parallel.For<double>(0,x.Length,() => 0.0, (i, ls, sum) =>
            {
                sum += x[i];
                return sum;
            }, sum => { lock (syncObj) mean += sum; });

            return mean / x.Length;
        }

        public static double pforeach_mean(double[] x)
        {
            double mean = 0.0;

            object syncObj = new object();

            Parallel.ForEach<double, double>(x, () => 0.0, (a, state, sum) =>
            {
                sum += a;
                return sum;
            }, sum => { lock (syncObj) mean += sum; });

            return mean / x.Length;
        }

        public static double part_mean2(double[] x, int start, int end)
        {
            double mean = 0.0;

            for (int i = start; i < end; i++)
            {
                mean += x[i];
            }

            return mean;
        }

        public static double part_mean(double[] x)
        {
            double left = 0.0, right = 0.0;

            Parallel.Invoke(() => left = part_mean2(x, 0, x.Length / 2),
                            () => right = part_mean2(x, x.Length / 2, x.Length));

            return (left + right) / x.Length;
        }

        #endregion

        #region Timing
        static double Time(Action action)
        {
            var sw = Stopwatch.StartNew();
            for (int i = 0; i < loop_count; i++)
            {
                action();
            }
            return sw.Elapsed.TotalMilliseconds / loop_count;
        }

        static double pTime(Action action)
        {
            var sw = Stopwatch.StartNew();
            Parallel.For(0, loop_count, (a) => action());
            return sw.Elapsed.TotalMilliseconds / loop_count;
        }
        #endregion

        #region Stdev
        public static void v2_part(double[] x, double mean, out double v1, out double v2, int start, int end)
        {
            v1 = 0.0; v2 = 0.0;
            for(int i = start; i < end; i++) {
                v2 = v2 + (x[i] - mean);
                v1 = v1 + (x[i] - mean) * (x[i] - mean);
            }
        }

        public static double fv1(double[] x, double mean)
        {
            double v1 = 0.0;
            foreach(double a in x) {
                v1 += (a - mean) * (a - mean);
            }

            return v1;
        }

        public static double fv2(double[] x, double mean)
        {
            double v2 = 0.0;
            foreach(double a in x) {
                v2 += (a - mean);
            }
            return v2;
        }

        public static void calc_stddev(double[] x, int n, out double stddev, out double mean, bool fmean)
        {

            double v1 = 0;
            double v2 = 0;
            double variance;
            mean = 0.0; stddev = 0.0;

            if (fmean)
            {
                for (int i = 0; i < x.Length; i++)
                {
                    mean += x[i];
                }
                mean /= x.Length;
            }

            //
            // Variance (using corrected two-pass algorithm)
            //
            if (n != 1)
            {
                v1 = 0; v2 = 0;

                for (int i = 0; i < n; i++)
                {
                    v2 += (x[i] - mean);
                    v1 += (x[i] - mean) * (x[i] - mean);
                }

                v2 = v2 * v2 / n;
                variance = (v1 - v2) / (n - 1);
                if ((double)(variance) < (double)(0))
                {
                    variance = 0;
                }
                stddev = Math.Sqrt(variance);
            }
        }

        public static void part_stddev(double[] x, int n, out double stddev, double mean, bool fmean)
        {
            mean = 0.0; stddev = 0.0;

            if (fmean)
            {
                mean = part_mean(x);
            }

            //
            // Variance (using corrected two-pass algorithm)
            //
            if (n != 1)
            {
                double v1 = 0;
                double v2 = 0;
                double variance;

                double v1_left = 0.0, v1_right = 0.0;
                double v2_left = 0.0, v2_right = 0.0;

                Parallel.Invoke(() => v2_part(x, mean, out v1_left, out v2_left, 0, x.Length / 2),
                                () => v2_part(x, mean, out v1_right, out v2_right, x.Length / 2, x.Length));

                v1 = v1_left + v1_right;
                v2 = v2_left + v2_right;

                v2 = v2 * v2 / n;
                variance = (v1 - v2) / (n - 1);
                if ((double)(variance) < (double)(0))
                {
                    variance = 0;
                }
                stddev = Math.Sqrt(variance);
            }
        }

        public static void part2_stddev(double[] x, int n, out double stddev, double mean, bool fmean)
        {
            mean = 0.0; stddev = 0.0;

            if (fmean)
            {
                mean = part_mean(x);
            }

            //
            // Variance (using corrected two-pass algorithm)
            //
            if (n != 1)
            {
                double v1 = 0;
                double v2 = 0;
                double variance;

                Parallel.Invoke(() => v1 = fv1(x, mean),
                                () => v2 = fv2(x, mean));

                v2 = v2 * v2 / n;
                variance = (v1 - v2) / (n - 1);
                if ((double)(variance) < (double)(0))
                {
                    variance = 0;
                }
                stddev = Math.Sqrt(variance);
            }
        }
        #endregion

        public struct data { public string name; public double ms; };

        static void Main()
        {
#if(DEBUG)
                    Console.WriteLine("Warning: benchmarking in debug mode will not produce meaningful results!");
#endif

            double[] array = new double[array_size];
            for (int i = 0; i < array.Length; ++i)
            {
                array[i] = i%10 + 1000000000;
            }

            int [] nums = {0,1,2,3};

            List<data>[] all = new List<data>[nums.Length];

            foreach (int i in nums)
            {
                Console.WriteLine(loop_count + " loops over " + sizeof(double) * array_size / 1e6 + " MiB, RAM->CPU predicted memcopy time: " + (sizeof(double) * array_size) / 3.2e6 + " ms");

                var methods = ((MethodInfo[])typeof(Program).GetMethods(BindingFlags.Public | BindingFlags.Static)).Where((a) => (a.Name.EndsWith("_mean")));

                List<data> trial = new List<data>(methods.Count());

                foreach (var m in methods)
                {
                    trial.Add(new data { name = m.Name, ms = Math.Round(Time(() => m.Invoke(null, new object[] { array })),i) });
                }

                trial.Add(new data { name = "pTime(foreach)", ms = Math.Round(pTime(() => foreach_mean(array)),i) });

                all[i] = trial;

                double stddev = 0.0, mean = 0.0;

                Console.Write("part_stdev:\t" + Time(() => part_stddev(array, array.Length, out stddev, mean, true)));
                Console.WriteLine("\t " + stddev);
                Console.Write("part2_stdev:\t" + Time(() => part2_stddev(array, array.Length, out stddev, mean, true)));
                Console.WriteLine("\t " + stddev);
                Console.Write("calc_stdev:\t" + Time(() => calc_stddev(array, array.Length, out stddev, out mean, true)));
                Console.WriteLine("\t " + stddev);

                array_size /= 10;
                loop_count *= 10;

                Array.Resize<double>(ref array, array_size);
            }

            var q = from a in all[0]
                    join b in all[1] on a.name equals b.name
                    join c in all[2] on a.name equals c.name
                    join d in all[3] on a.name equals d.name
                    orderby a.ms descending
                    select new { a.name, m800 = a.ms, m80 = b.ms, m8 = c.ms, m08 = d.ms };

            Console.WriteLine(); Console.WriteLine("\t\t800MiB\t80MiB\t8MiB\t0.8MiB");
            foreach (var r in q) {
                Console.WriteLine(r.name + "\t" + r.m800 + "\t" + r.m80 + "\t" + r.m8 + "\t" + r.m08);
            }
        }
    }
}

The first few times I ran this I noticed a problem – no matter what changes I made, many of the functions all ran at around the same speed. With 100M doubles basically what is happening is that 800MiB is being shuffled across the front side bus from the memory into the CPU and it turns out the FSB speed is the limiting factor. Using CPU-Z we can see that my machine, while DDR2, has only a single DIMM and is thus only single channel, so while the memory is PC2-6400 with a theoretical max bandwidth of 6.4GB/s my machine can only achieve half that.

One way to test the actual memory bandwidth is with the STREAM Benchmark which gives the following results on my machine:

C:\Stream>stream5.8_omp-32.exe
Unknown Minor OS Version (6.1)
OS Build number 7600
Service Pack Level:
-------------------------------------------------------------
STREAM version $Revision: 5.8 $
-------------------------------------------------------------
This system uses 8 bytes per DOUBLE PRECISION word.
-------------------------------------------------------------
Array size = 10000000, Offset = 0
Total memory required = 228.9 MB.
Each test is run 40 times, but only
the *best* time for each is used.
-------------------------------------------------------------
Number of Threads requested = 2
-------------------------------------------------------------
Printing one line per active thread....
Printing one line per active thread....
-------------------------------------------------------------
Your clock granularity/precision appears to be 1 microseconds.
Each test below will take on the order of 67617 microseconds.
   (= 67617 clock ticks)
Increase the size of the arrays if this shows that
you are not getting at least 20 clock ticks per test.
-------------------------------------------------------------
WARNING -- The above is only a rough guideline.
For best results, please be sure you know the
precision of your system timer.
-------------------------------------------------------------
Function      Rate (MB/s)   Avg time     Min time     Max time
Copy:        3013.0788       0.0569       0.0531       0.1064
Scale:       3008.6747       0.0570       0.0532       0.0849
Add:         2904.1052       0.0908       0.0826       0.2312
Triad:       2906.7409       0.0886       0.0826       0.1591
-------------------------------------------------------------
Solution Validates
-------------------------------------------------------------

I decided to run the program multiple times with different sized input arrays from 800MiB down to 0.8MiB, noting that the L2 cache on this machine is 2MB so hopefully the input array will be cached with the smallest array test. I ran through 5,50,500 and 5000 loops as the array size decreased. Here are the results (all times in ms):

And for the STDEV:

We can clearly see the speed limitation imposed by the memory bus bandwidth with many implementations. For the 800,80,8 and 0.8 data sizes my machine’s memory bandwidth of 3.2GB/s predicts minimum memory copy times of 250, 25, 2.5 and 0.25 ms. One way to reduce this would be to copy less data, possibly using ints or even shorts. I didn’t want to re-implement using ints so I made a brief investigation of the available Generic Numeric programming techniques available. Of the two I found, MiscUtil and Generic Numerics only the second implemented operators which would allow me to use generics without changing the code, thus allowing it to be turned on and off for comparison. Conveniently the downloadable project includes a performance comparison between the straight code and generics, and here is the partial result:

Non-generic standard deviation of 1000000 int elements
Method TestListSigma took 00:00:00.0079411
Generic standard deviation of 1000000 int elements
Method TestListSigmaG took 00:00:00.0050658
Generic standard deviation of 1000000 int elements using wrapper structs
Method TestListSigmaG2 took 00:00:00.0147550
Generic standard deviation of 1000000 int elements using overflow checks
Method TestListSigmaG3 took 00:00:00.0059253

Of these methods it is SigmaG2 which offers the use of standard operators … and is also the slowest, taking nearly three times longer than the fastest method. I didn’t want to make code change just for generics and a factor of 3x slowdown would invalidate any performance differences between data types so I abandoned that approach and stuck with doubles.

Conclusion

Mean:

A quick scan of the mean results set shows that VS2010 was generally as fast as VS2008 on everything. Although the larger sizes do benefit a little from caching (able to beat their predicted times by small amounts) the test is primarily memory bandwidth limited. It is only with the smaller tests, where the dataset is cached, that we start to measure the performance of the actual code. However at no input array size do the parallel versions of the code match the straightforward implementations, except where runs are timed in parallel, rather than the actual function itself being parallel and then the boost is significant. For large data sets where there isn’t much work to be done there is no benefit to making the code parallel. Of interest is the poor performance of the Linq and Plinq variants which are vastly slower than any other method. MS has a paper describing when to use Plinq and when to use Parallel.For. In my testing Parallel.For and Parallel.ForEach also didn’t fare well with the various overheads, data partitioning, locking and extra function calls ruining their performance. While these do offer user defined partitioning (Partitioner.Create) Parallel.For and Parallel.Foreach were so slow that I didn’t bother trying anything other than default partitioning. However this is an interesting area and the VisualizePartioning example from the PFX team blog is worth a look. Here’s a screenshot showing how the Parallel.For is load balanced where there are 2 cores:

Stdev:

The story is different here where there is just enough work to be done for us to start to see a benefit to parallel programming although only with VS 2008. With VS 2010 the best implementation was single threaded and it took nearly 4x longer to run the part2_stdev function than VS 2008 did. VS 2010 did have the best performing single threaded results.

Overall

Be careful when performance testing to make sure you test what you planned to, for most array sizes this ended up being a test of how fast the Front Side Bus could push data from memery into the CPU rather than the code’s speed. The various for vs foreach vs function calls vs hoisting variants – little difference on this dataset, although foreach without hoisting or using function calls is consistently near the top. Don’t bother with parallel implementations unless there is a lot of varied work going on in each loop. Making a parallel version of calc_mean was pointless. Don’t convert code to use Linq if you are interested in performance as Linq and Plinq turned out to be the slowest variants by far, with the Plinq mean taking 20x as long as the single threaded foreach on 0.8MiB test. Lastly I’m not sure what happened with the 80MiB calc_stdev test where it took nearly 250ms, looks like an anomaly but really, what dataset would be complete without at least one unexplained anomaly?