所求问题：

给定任意单位轴q(向量),求向量p(x,y,z)(或点p)饶q旋转theta角度的变换后的新向量p'(或点p')：

1.用四元数工具：
-------------------------------------------------------------------------
结论：构造四元数变换p'= q*p*q-1，（p,q是由向量p,q扩展成的四元数）。那么，p'转换至对应的向量(或点)就是变换后的新向量p'(或点p')。

其中，p',q,p,q-1均为四元数。q由向量q扩展,为q=(cos(theta/2),sin(theta/2)*q)，p由向量p扩展,为p=(0,x,y,z),q-1为q的逆，因为q为单位四元数，所以q-1=q*=(cos(theta/2),-sin(theta/2)*q)。
-------------------------------------------------------------------------

（这个结论的证明过程可以在网上找到。这里略去。）
下面看其时间复杂度：

首先有个三角函数的计算时间，这个可以预先计算好，花费时间不计。考虑n个四元数相乘需进行4*4*(n-1)=16*(n-1)次乘法，15*(n-1)次加法，因为加法化费时间较少，这里仅考虑乘法。这里涉及到三个四元数的乘法，设一次乘法的时间为T,故花费16*2=32T

2.旋转矩阵工具：
-------------------------------------------------------------------------
结论：构造旋转矩阵变换Trot,则变换后的新向量p'(或点p')为p'= p*Trot

其中，p'（x',y',z',1）,p(x,y,z,1)为向量p'，p的4D齐次坐标表示，Trot =

|t*x*x + c,       t*x*y + s*z,      t*x*z - s*y,     0|
|t*x*y - s*z,    t*y*y + c,          t*y*z + s*x,    0|
|t*x*z + s*y,   t*y*z - s*x,       t*z*z + c,        0|
|0,                    0,                       0,                    1|   

c=cos(theta), s=sin(theta),t=1-c.
-------------------------------------------------------------------------
（这个结论的证明过程可以在网上找到。这里略去。）

下面看其时间复杂度：

三角函数的计算时间不计。矩阵本身的元素乘法主要是计算t*x和s*x之类，需进行12+3=15次乘法。两个矩阵相乘的需进行n*n*n次乘法,这里n=4，所以花费4*4*4=64次乘法时间,但这里有个注意的地方就是旋转矩阵的第4维无须考虑，即可简化为3X3的矩阵。故花费3*3*3=27次乘法时间，总共的时间为15+27=42次乘法时间。cost=42T.

比较来看，还是使用四元数的效率要高出一些，这在要变换的点的数目越大时，体现的越明显。实际上，有很多3D引擎为了利用四元数运算的高效率性，一般先将矩阵转换成四元数后进行运算再转回为矩阵后，再传给DirectX或OpenGL库函数。

关于四元数和矩阵在向量（方向）之间的进行插值的效率比较，见下篇：
探讨：物体绕任意向量的旋转-四元数法VS.旋转矩阵法的性能比较

a.透视投影一般的视景体为棱台，相机空间的物体会投影到视平面z=d，这里考虑左手坐标系，矩阵使用行优先方式。如图所示，由相似三角形知识可知相机空间中的物体投影到视平面上的坐标为：

xp = xc*(dx/zc)
yp = yc*(dy/zc)

其中,xc,yc,zc为相机空间坐标,xp,yp,zp为视平面坐标，dx，dy为x,y轴向的视距view distance，视平面到camera的距离,
故相机空间投影到视平面上的矩阵Tcp为：

|dx 0  0 0  |

|0  dy 0 0  |

|0  0   1 1  |

|0  0   0 0  |

(验证:Tcp右乘点p(xc,yc,zc,1)得点p'(xc*dx, yc*dy, zc, zc),转换为3D坐标为(xc*dx/zc, yc*dy/zc, 1),正确。)

********************************************************************
注：因为转换过程中点使用的是4D齐次坐标,所以最后需转换为3D坐标。4D齐次坐标(x,y,z,w)转换为3D坐标的方法为除以w分量，即对应3D坐标为(x/w,y/w,z/w)。
********************************************************************

考虑dx/zc和dy/zc项，如果dx != dy,则投影后x,y的比例会发生变化（原因：投影前坐标比例为xc/yc,投影后为xp/yp = xc*(dx/zc)/yc*(dy/zc) = xc*dx/yc*dy），从而投影后的图像的x,y比例会发生变形。

---------------------------------------------
结论1：所以，一般都会令d=dx=dy,即x,y向的视距相同。否则，图像失真。
---------------------------------------------

考虑视角（view angle,或视野filed of view）的问题，视角的大小不会影响到物体投影后的坐标，只会影响可视的范围。

在视距一样的情况下，x,y轴的视角可以不一样。如果一样，那么视平面就是一个正方形的。于是有：

tan(theta_x/2) = width_p/d
tan(theta_y/2) = height_p/d 

其中，theta_x，theta_y为x,y轴向的视角，width_p，height_p为视平面z=d的宽度（x轴）和高度(y轴)。
----------------------------------------------------------------
结论2：视平面的宽高比rp=width_p/height_p = tan(theta_x/2)/tan(theta_y/2)。
----------------------------------------------------------------

2.视平面到屏幕的变换

下面就是视平面到屏幕的变换了，这是一个2D到2D的变换（视平面的坐标需伸缩以填满整个屏幕空间,即在视平面中出现的所有的点要变换到屏幕上去，同时x,y轴向的比例还要维持和变换前一样，以免比例失真，同时视平面的坐标原点和屏幕中心点（x0=(width_s)/2, y0=(height_s)/2）对应），其实，就是一个坐标缩放加平移的过程:

xs = xp*kx + x0
ys = -yp*ky + y0
 
矩阵Tps为：

|kx     0      0 0  |

|0      -ky    0 0  |

|0      0       1 0  |

|x0   y0      0 1  |

(验证:Tps右乘点p(xp,yp,zp,1)得点p'(xp*kx + x0, -yp*ky + y0, zp, 1),转换为3D坐标为(xp*kx + x0, -yp*ky + y0, zp)，正确。)

其中，kx,ky(kx>0,ky>0)为x,y轴向的缩放因子（kx=(width_s)/(width_p)， ky = (height_s)/(height_p），和视距相同，kx,ky的值必须一样，否则图像的x,y比例又会发生变形。这里-yp*ky是因为一般屏幕的y轴是向下的，跟相机空间和视平面坐标系中的y轴方向相反。
------------------------------------------------------------------------------------------------
结论3: 一般令k=kx=ky，即x,y向的缩放因子相同。否则，图像失真。
于是有，width_s/width_p = height_s/height_p，变化等式得，rp = width_p/height_p = width_s/height_s = rs
所以，在x,y轴视距一样的情况下，要想最后变换到屏幕上的图像x,y比例不失真，必须rp=rs,即视平面的宽高比和屏幕的宽高比一样。
-----------------------------------------------------------------------------------------------

********************************************************************
注：若屏幕宽高为W,H，当W != H，即屏幕不为方形的时候，要确保投影到屏幕上的图像x,y比例不失真，则x,y轴视角(或视野FOV)肯定不能相等。
原因： 由结论2,3知，rp=width_p/height_p = tan(theta_x/2)/tan(theta_y/2)=width_s/height_s=rs=W/H。 故由W/H != 1 => theta_x != theta_Y.
********************************************************************

3.综合：

相机空间点p转换到屏幕空间点p'，变换公式为：

xs = xc*(dx/zc)*kx + x0 = xc*(d/zc)*k + x0
ys = -yc*(dy/zc)*ky + y0 = -yc*(d/zc)*k + y0

综合变换矩阵(相机空间到屏幕空间)Tcs为：
 
Tcs = Tcp*Tps =

|d*k    0       0 0  |

|0      -d*k    0 0  |

|x0     y0      1 1  |

|0      0         0  0|

 其中，d为视距，k为屏幕宽高比或视平面宽高比，x0,y0为屏幕中心，注：最后需转换为3D坐标。

(验证:Tcs右乘点p(xc,yc,zc,1)得点p'(xc*d*k + x0*zc, -yc*d*k + y0*zc, zc, zc),转换为3D坐标为(xc*(d/zc)*k + x0, -yc*(d/zc)*k + y0, 1)，正确。)

4.一般情形：
 ************************************
视距为1，x轴视角为90度，屏幕宽高为W,H.
************************************
 
代入d=1，theta_x = PI/2,x0= W/2,y0=H/2,则视平面宽高为width_p = 2。
要确保屏幕上的图像x,y比例不失真，即rs=rp,有
height_p = 2/rp=2/rs=2H/W,
k=kx=ky=width_s/width_p = W/2.

于是,矩阵为:

Tcs1 = 

|W/2    0       0 0  |

|0      -W/2    0 0  |

|W/2    H/2   1 1  |

|0      0         0 0  |

或

   |W/2    0                  0 0  |

|0      -H/2*(W/H)    0 0  |

|W/2    H/2              1 1  |

|0      0                    0 0  |

(可以看到，y轴的缩放因子中乘上了宽高比(aspect ratio))
 这个矩阵较常用。

---------------------
有什么问题,欢迎探讨.
 

// INCLUDES //////////////////////////////////////////////////////////////////////////
#include "stdio.h"
#include "malloc.h"

// DEFINES //////////////////////////////////////////////////////////////////////////
#define N 4

// TYPES //////////////////////////////////////////////////////////////////////////
typedef struct BLOCK_TYP
{
 int x,y,z;
 int area;
 int height;
 
 BLOCK_TYP *pNext;
} BLOCK, *BLOCK_PTR;

typedef struct TOWER_TYP
{
 int height;
 int num_block;

 BLOCK_PTR pTowerTop;
} TOWER, *TOWER_PTR;

// PROTOTYPES //////////////////////////////////////////////////////////////////////////
TOWER_PTR BuildTower(BLOCK_PTR pBlock, int n);

// FUNCTIONS //////////////////////////////////////////////////////////////////////////
int main(int argc, wchar_t* argv[])
{
 BLOCK blocks[N] = {{1,1,1},{2,2,2},{3,3,3},{4,4,4}};

 printf("We have %d kinds of blocks, they are:\n", N);
 int i;
 for(i = 0; i < N; i++)
 {
  printf("(%d,%d,%d)\n", blocks[i].x, blocks[i].y, blocks[i].z);
 }
 printf("Building the tower of babylon...\n");

 TOWER_PTR pT = BuildTower(blocks, N);
 printf("The max height of tower is: %d.\nTotally %d blocks are used, they are:\n", pT->height, pT->num_block );

 BLOCK_PTR pBlock = pT->pTowerTop;
 for(i = 0; i < pT->num_block; i++)
 {
        printf("(%d) block (%d,%d,%d), area = %d, height = %d\n", i+1, pBlock->x, pBlock->y, pBlock->z, pBlock->area, pBlock->height); 
  pBlock = pBlock->pNext;
 }

 getchar();

 return 0;
}

////////////////////////////////////////////////////////////////////////////////
// Algorithm of building the Tower Of Babylon.
// Input Parameters: pBlock: pointer variable that identifies blocks sequence.
// n: int variable that identifies how many kinds of blocks.
// Output Parameters: None.
// Return: pointer variable that identifies the built tower.
////////////////////////////////////////////////////////////////////////////////
TOWER_PTR BuildTower(BLOCK_PTR pBlock, int n)
{
 int index_block = 0;
 TOWER_PTR pTower = new TOWER();
 BLOCK_PTR pTop = new BLOCK(); 
 pTower->pTowerTop = pTop;

 // initialize tower
 pTower->pTowerTop->x = pBlock->x;
 pTower->pTowerTop->y = pBlock->y;
 pTower->pTowerTop->z = pBlock->z;

 pTower->pTowerTop->area = (pBlock->x)*(pBlock->y);
 pTower->pTowerTop->height = pBlock->z;
 pTower->height = pTower->pTowerTop->height;
 pTower->num_block = 1;
   
 for(int i = 1; i < 3*n; i++)
 {
  index_block = i/3;
  if (i%3 == 0) // condition 1, area = x*y, height = z.
  {
   (pBlock+index_block)->area = ((pBlock+index_block)->x)*((pBlock+index_block)->y);
   (pBlock+index_block)->height = (pBlock+index_block)->z;
  }
  else if (i%3 == 1) // condition 2, area = x*z, height = y.
  {
   (pBlock+index_block)->area = ((pBlock+index_block)->x)*((pBlock+index_block)->z);
   (pBlock+index_block)->height = (pBlock+index_block)->y;
  }
  else // condition 3, area = y*z, height = z.
  {
   (pBlock+index_block)->area = ((pBlock+index_block)->y)*((pBlock+index_block)->z);
   (pBlock+index_block)->height = (pBlock+index_block)->x;
  }
 
  bool bNeedToBeAdded = true; 

  BLOCK_PTR pB = pTower->pTowerTop;
  BLOCK_PTR pPrev = pTower->pTowerTop;
  while(pB != NULL)
  {
   if ((pBlock+index_block)->area < (pB->area))
   {
 // insert new block
 BLOCK_PTR pNewBlock = new BLOCK(); 
    *pNewBlock = *(pBlock+index_block);

 if(pB == pPrev)
 {
  pNewBlock->pNext = pB;
  pTower->pTowerTop = pNewBlock;
 }
 else
 {
  pNewBlock->pNext = pPrev->pNext;
  pPrev->pNext = pNewBlock;
 }
 
    // increase number of blocks
    pTower->num_block++;
    // increase height of tower
    pTower->height += (pBlock+index_block)->height;
   
    bNeedToBeAdded = false;
    break;
   }
   else if ((pBlock+index_block)->area == (pB->area))
   {
    if (pB->height < ((pBlock+index_block)->height))
    {
     // increase height of tower
     pTower->height += (pBlock+index_block)->height - pB->height;
     // replace blocks
  BLOCK_PTR pNewBlock = new BLOCK(); 
  *pNewBlock = *(pBlock+index_block);

  if (pB == pPrev)
  {
  pNewBlock->pNext = pB->pNext;
  pTower->pTowerTop = pNewBlock;
  }
  else
  {
   pNewBlock->pNext = pB->pNext;
   pPrev->pNext = pNewBlock;
  }
    }

    bNeedToBeAdded = false;
    break;
   }

   pPrev = pB;
   pB = pB->pNext;
  }
 
  if(bNeedToBeAdded)
  {
   // add new block at the end
   BLOCK_PTR pNewBlock = new BLOCK(); 
   *pNewBlock = *(pBlock+index_block);
  
   pPrev->pNext = pNewBlock;
   pNewBlock->pNext = NULL;

   // increase number of blocks
   pTower->num_block++;
   // increase height of tower
   pTower->height += (pBlock+index_block)->height;
  }
 }

 return pTower;
}

// INCLUDES ///////////////////////////////////////////////
#include "stdio.h"

// DEFINES ////////////////////////////////////////////////
#define N 11

// PROTOTYPES /////////////////////////////////////////////
void QuickSort(int *seq, int lenght);
void InsertSort(int *seq, int lenght);
void InsertSort2(int *seq, int lenght);
void ShellSort(int *seq, int lenght);
void IncrementInsertSort(int *seq, int length, int increment);
void HeapSort(int *seq, int length);
void SelectionSort(int *seq, int length);
void BubbleSort(int *seq, int length);
void BiBubbleSort(int *seq, int length);
void AdjustHeap(int *seq, int startIndex, int endIndex);

// GLOBALS ////////////////////////////////////////////////
int nAssignmentOperation = 0; // assignment counter
int nComparsionOperation = 0; // comparsion counter

// FUNCTIONS //////////////////////////////////////////////

int main(int argc, char* argv[])
{
 printf("0: Original Sequence\n");
 printf("1: Quick Sort\n");
 printf("2: Insert Sort\n");
 printf("3: Shell Sort\n");
 printf("4: Insert Sort2\n");
 printf("5: Heap Sort\n");
 printf("6: Selection Sort\n");
 printf("7: Bubble Sort\n");
 printf("8: Bi-Bubble Sort\n");
 printf("-1: Exit\n"); 

 int cat = 0; 
 while (cat != -1)
 {
  // clear counter
  nAssignmentOperation = 0;
  nComparsionOperation = 0;

  //int array[N] = {600, 64,6,78,246,445,345,445,7,4885,1,600, 64,6,78,246,445,345,445,7,4885,1,600, 64,6,78,246,445,345,445,7,4885,1,600, 64,6,78,246,445,345,445,7,4885,1,600, 64,6,78,246,445,345,445,7,4885,1,
  //600, 64,6,78,246,445,345,445,7,4885,1,600, 64,6,78,246,445,345,445,7,4885,1,600, 64,6,78,246,445,345,445,7,4885,1,600, 64,6,78,246,445,345,445,7,4885,1,600, 64,6,78,246,445,345,445,7,4885,1};
  int array[N] = {600, 64,6,78,246,445,345,445,7745,4885,1};

  printf("\nPlease select sort method:");
  scanf("%d", &cat);
  printf("-------------------------------------------------\n");

  switch(cat)
  {
   case 0:printf("No Sort. The Original Sequence is\n");
    break;
   case 1:printf("Quick Sort\n"); QuickSort(array, N);
    break;
   case 2:printf("Insert Sort\n"); InsertSort(array, N);
    break;
   case 3:printf("Shell Sort\n"); ShellSort(array, N);
    break;
   case 4:printf("Insert Sort2\n"); InsertSort2(array, N);
    break;
   case 5:printf("Heap Sort\n"); HeapSort(array, N);
    break;
   case 6:printf("Selection Sort\n"); SelectionSort(array, N);
    break;
   case 7:printf("Bubble Sort\n"); BubbleSort(array, N);
    break;
   case 8:printf("Bi-Bubble Sort\n"); BiBubbleSort(array, N);
    break;
   default:printf("Exit...\n");
    return 0;
  }

  for(int i = 0; i < N; i++)
  {
   printf("int = %d\n", array[i]);
  }

  printf("Count of comparsion operation  = %d\n", nComparsionOperation);
  printf("Count of assignment operation  = %d\n", nAssignmentOperation);
  printf("-------------------------------------------------\n");
 }

 return 0;
}

inline void Swap(int *pCurrPost, int *pCurrKey)
{
 nAssignmentOperation += 3;
 int temp;

 //swap value
 temp = *pCurrPost;
 *pCurrPost = *pCurrKey;
 *pCurrKey = temp;
}

////////////////////////////////////////////////////////////////////////////////
// Quick Sort algorithm.  (By SoRoMan, 2006.7.31)
// 快速排序的基本思想是基于分治策略的。对于输入的子序列L[p..r]，如果规模足够小则直接进行排序（比如用前述的冒泡、选择、插入排序均可），否则分三步处理：
// 1.分解(Divide)：将待排序列L[p..r]划分为两个非空子序列L[p..q]和L[q+1..r]，使L[p..q]中任一元素的值不大于L[q+1..r]中任一元素的值。具体可通过这样的途径实现：在序列L[p..r]中选择数据元素L[q]，经比较和移动后，L[q]将处于L[p..r]中间的适当位置，使得数据元素L[q]的值小于L[q+1..r]中任一元素的值。
// 2.递归求解(Conquer)：通过递归调用快速排序算法，分别对L[p..q]和L[q+1..r]进行排序。
// 3.合并(Merge)：由于对分解出的两个子序列的排序是就地进行的，所以在L[p..q]和L[q+1..r]都排好序后不需要执行任何计算L[p..r]就已排好序，即自然合并。
// 这个解决流程是符合分治法的基本步骤的。因此，快速排序法是分治法的经典应用实例之一。
// Input Parameters: seq: pointer variable that identifies which sequence will be sorted.
// length: int that identifies the length of the sequence.
// Output Parameters: None.
// Return: None.
////////////////////////////////////////////////////////////////////////////////
void QuickSort(int *seq, int length)
{
 if(length <= 1)
 {
  return;
 }
 else
 {
  int post = *seq; // set post
  int *pCurrPost = seq; // set current position of post
  int *pCurrKey; // current position of key

  for(int j = length - 1, i = 1; i <= j;)
  {
   // handle the right sub-sequence
   while(i <= j)
   {
    pCurrKey = seq + j;

    if(*pCurrKey < post)
    {
     Swap(pCurrPost, pCurrKey);
     pCurrPost = pCurrKey;

     break;
    }
    else
    {
     j--;
    }

    nComparsionOperation ++;
   }

   // handle the left sub-sequence
   while(i <= j)
   {
    pCurrKey = seq + i;

    if(*pCurrKey > post)
    {
     Swap(pCurrPost, pCurrKey);
     pCurrPost = pCurrKey;

     break;
    }
    else
    {
     i++;
    }

    nComparsionOperation ++;
   }
  }
  // handle two sub sequences recursively
  int lleng = pCurrPost - seq; // left sub sequence
  int rleng = length - lleng - 1; // right sub sequence
  QuickSort(seq, lleng);
  QuickSort(seq + lleng + 1, rleng);
 }
}

////////////////////////////////////////////////////////////////////////////////
// Insert Sort algorithm (A special increment insert sort, increment = 1).  (By SoRoMan, 2006.7.31)
//插入排序的基本思想是，经过i-1遍处理后,L[1..i-1]己排好序。第i遍处理仅将L[i]插入L[1..i-1]的适当位置，使得L[1..i]又是排好序的序列。要达到这个目的，我们可以用顺序比较的方法。首先比较L[i]和L[i-1]，如果L[i-1]≤ L[i]，则L[1..i]已排好序，第i遍处理就结束了；否则交换L[i]与L[i-1]的位置，继续比较L[i-1]和L[i-2]，直到找到某一个位置j(1≤j≤i-1)，使得L[j] ≤L[j+1]时为止。
//简言之,插入排序就是每一步都将一个待排数据按其大小插入到已经排序的数据中的适当位置，直到全部插入完毕。插入排序方法分直接插入排序和折半插入排序两种，这里只介绍直接插入排序，折半插入排序留到“查找”内容中进行
// Input Parameters: seq: pointer variable that identifies which sequence will be sorted.
// length: int that identifies the length of the sequence.
// Output Parameters: None.
// Return: None.
////////////////////////////////////////////////////////////////////////////////
void InsertSort(int *seq, int length)
{
 IncrementInsertSort(seq, length, 1);
}

void InsertSort2(int *seq, int length)
{
 for(int i = 1; i < length; i++)
 {
  for (int j = 0; j < i; j++)
  {
   if(*(seq + i) < *(seq + j))
   { 
    int temp = *(seq + i);
    nAssignmentOperation ++;

    // insert, move (i-j) items
    for(int k = i; k > j; k-- )
    {
     *(seq + k) = *(seq + k - 1);
     nAssignmentOperation ++;
    }

    *(seq + k) = temp;
    nAssignmentOperation ++;

    nComparsionOperation ++;

    break;
   }
  }
 }
}

void IncrementInsertSort(int *seq, int length, int increment)
{
 for(int k = 0; k < increment; k++)
 {
  for (int i = increment + k; i < length; i+=increment + k)
  {
   int temp = *(seq + i);
   for (int j = i; j > increment - 1;)
   {
    nComparsionOperation ++;

    if(temp > *(seq + j - increment))
    {    
     *(seq + j) = *(seq + j - increment);
     nAssignmentOperation ++;
     j-=increment;
    }
    else
    {
     break;
    }
   }

   *(seq + j) = temp;
   nAssignmentOperation ++;    
  }
 } 
}

////////////////////////////////////////////////////////////////////////////////
// Shell Sort algorithm (Increment insert sort, increment = increment/3 + 1.).  (By SoRoMan, 2006.7.31)
// Input Parameters: seq: pointer variable that identifies which sequence will be sorted.
// length: int that identifies the length of the sequence.
// Output Parameters: None.
// Return: None.
////////////////////////////////////////////////////////////////////////////////
void ShellSort(int *seq, int length)
{
 for(int increment = length; increment > 1;)
 {
  increment = increment/3 + 1;
  IncrementInsertSort(seq, length, increment);
 }
}

////////////////////////////////////////////////////////////////////////////////
// Heap Sort algorithm.  (SoRoMan, 2006.7.31)
//
//1.堆的概念：
//堆是一个关键字序列{k1,k2,…,kn},它具有如下特性：
//ki ≤k2i, ki ≤ k2i+1
//或ki ≥ k2i,ki ≥ k2i+1(i=1,2,…, └n/2┘)
//堆的特性在完全二叉树中解释为：完全二叉树中任一结点的关键字都小于等于（或大于等于）它的左右孩子的关键字。
//如下列关键字序列均是堆： {5,23,16,68,94,72,71,73} (小顶堆) {96,83,27,38,11,9} （大顶堆） 由堆的定义可知：若关键字序列{k1,k2,…,kn}是堆，则堆顶元素是n个元素序列中的最小（或最大）元素。

//2.基本思想：
//首先将初始待排序记录序列建堆，则堆顶元素必为含最大关键字或最小关键字的记录，输出该记录，然后将剩余记录再调整为堆，如此反复，直至排序结束。
//在堆排序主要需解决两个问题：
//① 如何建堆？ 在完全二叉树中，所有序号i>└n/2┘的结点都是叶子，因此，以这些结点为根的子树均已是堆，这样只需将以序号为└n/2┘, └n/2┘-1, └n/2┘-2,…,1的结点为根、的子树都调整为堆即可。在按此次序调整每个结点时，其左、右子树均已是堆。
//② 若ki的左、右子树已经是堆，如何将以ki为根的完全二叉树也调整为堆？ 因ki的左、右子树已经是堆，所以必须在ki 和它的左、右孩子中选出最小（或最大）的结点放到ki的位置上，不妨设k2I关键字最小，将ki与k2I交换位置，而交换之后有可能导致以k2I为根的子树不再为堆，于是可重复上述过程，将以k2I为根的子树调整为堆，……,如此逐层下去，最多可能一直调整到树叶，此方法称为"筛选法"。

//3.性能分析：
//堆排序的时间主要由建立初始堆和不断重复建堆这两部分的时间开销构成。其中：建立初始堆时，总共需进行的关键字比较次数 ≤ 4n =O(n) ;排序过程中进行n-1次重建堆所进行的总比较次数不超过下式： 2*(└ log2(n-1) ┘+└ log2(n-2) ┘+…+ └ log22 ┘) < 2n └ log2n ┘=O(n log2n)因此堆排序总的时间复杂度是 O(n+ n log2n)= O(n log2n)。
//堆排序在最坏情况下的时间复杂度也是O(nlog2n)，相对于快速排序来说，这是堆排序的最大优点。
//另外，堆排序仅需一个记录大小供交换用的辅助空间。由于建初始堆所需的比较次数较多，所以堆排序不适宜于记录较少的文件，但对n较大的文件还是很有效的。
// Input Parameters: seq: pointer variable that identifies which sequence will be sorted.
// length: int that identifies the length of the sequence.
// Output Parameters: None.
// Return: None.
////////////////////////////////////////////////////////////////////////////////
void HeapSort(int *seq, int length)
{
 for(int n = length; n > 0; n--)
 {
  for(int j = n/2; j > 0; j--)
  {
   AdjustHeap(seq, j, n);
  }

  Swap(seq, seq+n-1);  
 } 
}

void AdjustHeap(int *seq, int startIndex, int endIndex)
{
 while(2*startIndex <= endIndex)
 {
  int i =  startIndex - 1;
  startIndex = startIndex*2;

  nComparsionOperation ++;
  if (startIndex < endIndex && *(seq + startIndex -1) < *(seq + startIndex))
  {
   startIndex++;
  }

  nComparsionOperation ++;
  if(*(seq + startIndex -1) > *(seq + i))
  {
   Swap(seq + startIndex - 1, seq + i);
  }  
  else
  {
   break;
  }
 }
}

////////////////////////////////////////////////////////////////////////////////
// Selection Sort algorithm.  (SoRoMan, 2006.7.31)
// Input Parameters: seq: pointer variable that identifies which sequence will be sorted.
// length: int that identifies the length of the sequence.
// Output Parameters: None.
// Return: None.
////////////////////////////////////////////////////////////////////////////////
void SelectionSort(int *seq, int length)
{
 for(int i = 0; i < length; i++)
 {
  int k = i;
  for(int j = i+1; j < length; j++)
  {
   nComparsionOperation ++;
   if (*(seq + j) < *(seq + k))
   {
    k = j;
   }
  }

  Swap(seq + k, seq + i);
 }
}

////////////////////////////////////////////////////////////////////////////////
// Bubble Sort algorithm.  (SoRoMan, 2006.7.31)
// Input Parameters: seq: pointer variable that identifies which sequence will be sorted.
// length: int that identifies the length of the sequence.
// Output Parameters: None.
// Return: None.
////////////////////////////////////////////////////////////////////////////////
void BubbleSort(int *seq, int length)
{
 for(int i = 0; i < length; i++)
 {
  for(int j = length-1; j > i; j--)
  {
   nComparsionOperation ++;
   if (*(seq + j) < *(seq + j - 1))
   {
    Swap(seq + j, seq + j - 1);
   }
  }
 }
}

////////////////////////////////////////////////////////////////////////////////
// Bi-Directinal Bubble Sort algorithm.  (SoRoMan, 2006.7.31)
//思想：
//基于在一趟排序中，位于最后一次交换关键字的的后面的关键字序列已经是有序的，所以不用再排序。
//相对于一般的冒泡排序，可以减少关键字比较次数，但交换次数不变。
// Input Parameters: seq: pointer variable that identifies which sequence will be sorted.
// length: int that identifies the length of the sequence.
// Output Parameters: None.
// Return: None.
////////////////////////////////////////////////////////////////////////////////
void BiBubbleSort(int *seq, int length)
{
 int k,low,up;
 low = 0;
 up = length-1;
 while (up > low)
 {
  k = low;  

  for(int i = low; i < up; i++)
  {
   nComparsionOperation ++;
   if (*(seq + i) > *(seq + i + 1))
   {
    Swap(seq + i, seq + i + 1);

    k = i;
   }
  }

  up = k;

  for(int j = up; j > low; j--)
  {
   nComparsionOperation ++;
   if (*(seq + j) < *(seq + j - 1))
   {
    Swap(seq + j, seq + j - 1);

    k = j;
   }
  } 

  low = k;
 }

}

基本上Bresenham画线算法的思路如下:

// 假设该线段位于第一象限内且斜率大于0小于1,设起点为(x1,y1),终点为(x2,y2).
// 根据对称性,可推导至全象限内的线段.
1.画起点(x1,y1).
2.准备画下个点。x坐标增1，判断如果达到终点，则完成。否则，由图中可知,下个要画的点要么为当前点的右邻接点,要么是当前点的右上邻接点.
2.1.如果线段ax+by+c=0与x=x1+1的交点的y坐标大于M点的y坐标的话,下个点为U(x1+1,y1+1)
2.2.否则,下个点为B(x1+1,y1+1)
3.画点(U或者B).
4.跳回第2步.
5.结束.

这里需要细化的是怎么判断下个要画的点为当前点的右邻接点还是当前点的右上邻接点.

设线段方程：ax+by+c=0(x1<x<x2,y1<y<y2)
令dx=x2-x1,dy=y2-y1
则：斜率-a/b = dy/dx.

从第一个点开始，我们有F(x,1,y1) = a*x1+b*y1+c=0
下面求线段ax+by+c=0与x=x1+1的交点:
由a*(x1+1)+b*y+c = 0, 求出交点坐标y=(-c-a(x1+1))/b
所以交点与M的y坐标差值Sub1 = (-c-a(x1+1))/b - (y1+0.5) = -a/b-0.5，即Sub1的处始值为-a/b-0.5。

则可得条件当 Sub1 = -a/b-0.5>0时候,即下个点为U.
反之,下个点为B.
代入a/b,则Sub1 = dy/dx-0.5.

因为是个循环中都要判断Sub,所以得求出循环下的Sub表达式,我们可以求出Sub的差值的表达式.下面求x=x1+2时的Sub，即Sub2
1.如果下下个点是下个点的右上邻接点，则
Sub2 = (-c-a(x1+2))/b - (y1+1.5) = -2a/b - 1.5
故Sub差值Dsub = Sub2 - Sub1 = -2a/b - 1.5 - (-a/b-0.5) = -a/b - 1.代入a/b得Dsub = dy/dx -1;
2.如果下下个点是下个点的右邻接点，
Sub2 = (-c-a(x1+2))/b - (y1+0.5) = -2a/b - 0.5
故Sub差值Dsub = Sub2 - Sub1 = -2a/b - 0.5 - (-a/b-0.5) = -a/b. 代入a/b得Dsub = dy/dx;

于是，我们有了Sub的处始值Sub1 = -a/b-0.5 = dy/dx-0.5,又有了Sub的差值的表达式Dsub = dy/dx -1 （当Sub1 > 0）或 dy/dx（当Sub1 < 0）.细化工作完成。

于是pcode可以细化如下：

// Pcode for Bresenham Line
// By SoRoMan
x=x1;
y=y1;
dx = x2-x1;
dy = y2-y1;
Sub = dy/dx-0.5; // 赋初值，下个要画的点与中点的差值

DrawPixel(x, y); // 画起点
while(x<x2)
{
 x++;
 if(Sub > 0) // 下个要画的点为当前点的右上邻接点
 {
  Sub += dy/dx - 1; //下下个要画的点与中点的差值
  y++; // 右上邻接点y需增1
 }
 else// 下个要画的点为当前点的右邻接点
 {
  Sub += dy/dx; 
 }
// 画下个点
DrawPixel(x,y);
}

PS:一般优化：
为避免小数转整数以及除法运算，由于Sub只是用来进行正负判断，所以可以令Sub = 2*dx*Sub = 2dy-dx,则
相应的DSub = 2dy - 2dx或2dy.

思考1:如果Sub = 0时，会产生取两个点都可以的问题。这个问题还没深入。

思考：计算机图形学中的左上像素填充约定（top-left filling convention for filling geometry）

我们知道目前显示器显示的原理，光栅化显示最小的单元是pixel(或者说精度为piexl），一个pixel所占的区域一般是个正方形，显示器上的画面就由一个或多个这样的小方格组成。如果以每个小方格为一个单元来制定屏幕坐标的话，问题就来了，因为方格本身有大小，不是一个严格意义上的"点"。

问题来源：屏幕坐标点的选取以及像素方格本身的大小
如果以方格中心点为基准制定坐标点（即对于任意一整数坐标点p(x,y),p必为某个像素方格的中心），以一个方格为跨度来制定屏幕坐标系的话。如图，屏幕上3X3的方格组成一个正方形的画面，如果按照方格中心点为基准制定坐标点，其区域屏幕坐标为（0，0，2，2）（注意：0，1，2都位于方格中心，其中，左上角坐标为(0,0)，右下角坐标为（2，2）），而按照数学运算，其上下跨距就是2-0 = 2 pixels，左右跨距2-0 = 2 pixels，即总共包含2X2=4个方格.而实际画面中的跨距是3X3，即共包含9个方格（从最左上端到最右下端），问题来了，前后矛盾了。左上像素填充约定可以解决这个问题。

问题解决：左上像素填充约定
还是举刚才的例子，如果给定左上角坐标为(0,0)，右下角坐标为（2，2）的矩形区域，叫显示设备去绘制的话，按照左上像素填充约定，在屏幕上填充的像素点将是:（0，0）,（0，1）,（1，0）,（1，0）这个2X2区域，而非图中所示的3X3区域。反过来说图中所示的3X3区域要用坐标表示的话应该是（0，0，3，3）。简单来说，在左上像素填充约定下，要填充（x1,y1,x2,y2）（x1,y1,x2,y2为整数）的矩形，实际填充的区域会是（x1-0.5,y1-0.5,x2-0.5,y2-0.5）,这个区域的左上角像素坐标为（x1,y1）,右下角像素坐标为(x2-1,y2-1)，即实际的填充区域会向左和向上各偏移0.5个像素。故称左上像素填充约定。

  0 1 2
   0 口口口
   1 口口口
   2 口口口

D3D, GDI, OpenGL等图形库使用左上填充约定的，如果接触过这些图形库，你可能会知道这个约定。还记得画一个全屏幕的矩形吧：DrawRectangle(0, 0, Screen_Width-1, Screen_Height-1)。

想法一：想到别的填充约定，比如右上，左下，右下等，这些不知道有哪些系统使用。

想法二：如果以像素方格的左上角为基准制定坐标点（即对于任意一整数坐标点p(x,y),p必为某个像素方格的左上角），情况怎么样？这时，给定（0，0，2，2），则填充时可以还会遇到一样的问题，填充（0，0）,（0，1）,（1，0）,（1，0），（2，0），（2，1），（0，2），（1，2），（2，2），最后还是包含了3X3= 9个像素。

想法三：
to be continued...