使用汇编和 SSE2/SSE3 指令进行绘图优化
4.92/5 (24投票s)
介绍用于加速预绘图过程的算法和技术。
引言
本文的目的是证明最新的并不意味着更快或更好，并且汇编对于时间关键的代码仍然非常有用。因为用户不喜欢等待 - 绘图以及所有与数据预绘图处理相关的操作都应该以超音速完成。但是，不可能对绘图本身做任何事情，无论您使用 GDI、GDI+、Direct X 还是 OpenGL 进行绘图，这都无关紧要。但优化预绘图数据处理是可能的。我说预绘图处理是什么意思？预绘图处理包括以下步骤：
- 仿射变换
- 裁剪
- 绘制点简化
因此，让我们看看标准 .Net 库中可以用于上述操作的内容。
- 对于仿射变换，最佳候选是
System::Drawing::Drawing2d::Matrix。但不幸的是，它的变换速度太慢了。 - 对于剪裁操作，我没有找到任何东西。
System::Drawing::Graphics具有ClipBounds属性，但如果您尝试绘制一个坐标超出可见区域的图元，它就会崩溃。 - 对于简化，标准库也没有提供任何东西。
结果是所有预绘图操作都让开发者头疼，在这篇文章中，我将尝试展示如何将这种头疼减到最少。
背景
让我们来看下一个任务，输入是一个由大量点表示的曲线，比如说 1,000,000 个点。我们需要交互式地绘制这条曲线。所以所有负责准备原始曲线以供绘制的操作都应该非常快速地完成。这些预绘图操作包括引言段中提到的所有三个操作。为了解决这个任务，可以假设几种方法。
创建您自己的 Matrix 类来高速执行仿射变换。创建功能或类来执行剪裁。创建功能或类来执行点简化。当然，这种方法可以符合 OOP 思想并且是通用的，并且将为每个步骤解决速度问题。但是，让我们看看这种实现的缺点。每个步骤都需要为存储结果而进行内存分配，当然当前步骤的结果将被用作下一个步骤的输入数据。因此，总的来说，在最坏的情况下，这种实现至少需要三次内存分配操作，并且在最坏的情况下需要遍历所有点三次，所以总时间复杂度将是 O(3n)。
我想介绍一种可能不那么通用和面向对象的方法，但它没有上述方法的缺点。主要思想是在一个循环中执行所有预绘图操作，因此每个循环迭代都会执行变换、剪裁和简化操作。关键点是尽可能少地执行内存重新分配，并尽可能少地访问内存。当然，动态内存分配是好的，但另一方面，它会导致额外的操作来寻找合适大小的内存块进行分配。因此，例程中的内存（重新）分配越少，性能就越好。
那么，让我们详细看看这种方法。我想介绍的第一个类是托管的 Transformer
第二个要介绍的是一组用汇编语言实现的非托管函数，用于执行仿射变换、剪裁、结果点简化和输出。实现这些函数要点：
- 不执行任何动态内存分配。
- 在一个循环中执行变换、剪裁和简化。因此，实现的算法时间复杂度为 O(n)。
- 减少内存访问次数（汇编中的远地址寻址）。这是执行速度的一个相当重要的因素，因为寄存器访问比内存访问快很多倍。因此，在最坏的情况下，每次迭代内存最多只访问 4 次。
- 将当前点加载到 XMM 寄存器。由于 SSE2 指令，这只需要一次内存访问。
- 写入结果点。这得益于 MMX 指令。
- 增加结果点的计数器，因为该值位于堆栈上。
- 设置不连续标记。（仅当有必要时）。
- 最小化内存访问的另一个方面是，变换矩阵值和剪裁边界在主循环之前只加载到 xmm 寄存器一次。这对于执行时间非常重要，因为编译器无法进行这种优化，并且可以减少循环迭代期间的内存访问。如果您用 C、C++ 或 C# 或任何其他高级编程语言编写代码，这种优化是不可能的。
- 如上所述，MMX/SSE2/SSE3 指令的使用允许缓存每个迭代使用的数据到 mmx/xmm 寄存器中，并提高计算速度和计算精度。
现在是时候仔细看看算法的主要参与者了。在下面的代码示例中，我将描述实现的关键点以及为什么这样做的原因。因此，下面的代码是 TransfromClipReducePoints 函数的一部分，负责初始化。
;load transformation matrix to xmm0-xmm2 registers in a way:
; register Low | High
;-------------------------
; xmm0 m00 | m10
; xmm1 m01 | m11
; xmm2 m20 | m21
movsd xmm0,_m00
movhpd xmm0,_m10
movsd xmm1, _m01
movhpd xmm1, _m11
movsd xmm2, _m20
movhpd xmm2, _m21
;loading clipping bounds to xmm6, xmm7 registers in a way
; register Low | High
;-------------------------
; xmm6 xMin | yMin
; xmm7 xMax | yMax
movsd xmm6, _xMin
movhpd xmm6, _yMin
movsd xmm7, _xMax
movhpd xmm7, _yMax
;discontinuity marker loading to mm7 register
mov eax, 80000000h
push eax
push eax
movq mm7, [esp]
add esp,8h
xor eax, eax
push eax ;stack allocation for output points count.
xor ebx, ebx ;initialization to store discontinuity flag.
要提及的关键点：
- 变换矩阵加载方式。
原始矩阵。m00 m01 m10 m11 m20 m21 xmm0 低 xmm0 高 xmm1 低 xmm1 高 xmm2 低 xmm2 高 m00 m10 m01 m11 m20 m21 元素如何进入 xmm 寄存器。
矩阵加载方式是为了减少每个点加载时重复的操作次数。因为输入数组中的数据是 x0,y0,x1,y1,...xn,yn，所以可以通过一次内存访问将一个点的坐标 x,y 加载到 xmm 寄存器中。并且因为矩阵以上述方式加载到
xmm寄存器中，所以可以按原样加载数据，而无需额外的重组来进行计算。 不连续标记是输出点坐标的定义值，用于指定曲线超出剪裁区域。值
80000000h对应于托管代码的int::MinValue。标记值存储在mm7寄存器中，以便一次内存访问输出。
下一部分代码对每个点执行仿射变换。
;---------------------------------------------------------------------
;Function performs affine transformation of point located by ptr[esi]
;increments [esi] to point no next point
;calculated point is stored in xmm3 register, in a way Low - X, High-Y
;---------------------------------------------------------------------
TransformPoint PROC PRIVATE
movupd xmm3, [esi] ;xmm3 = x,y
movapd xmm4,xmm3 ;xmm4 = x,y
;Calculation
mulpd xmm3, xmm0 ;xmm3 = M00*X | M10*Y
mulpd xmm4, xmm1 ;xmm4 = M01*X | M11*Y
haddpd xmm3,xmm4 ;xmm3 = M00*X + M10*Y | M01X+M11Y
addpd xmm3, xmm2 ;xmm3 = M00*X + M10*Y+ M20 | M01*X + M11*Y + M20
add esi, 10h ;increment esi till next call of this function.
ret 0
TransformPoint ENDP要提及的关键点：
该函数使用
SSE2/SSE3指令计算两个坐标。此函数的实现方式对矩阵元素的加载方式至关重要。
在整个代码中，只有此函数执行输入数据读取，这就是为什么它会递增
esi寄存器（算法中的第一次内存访问）。此函数的输出存储在xmm3寄存器中，供将来使用。该函数仅在输入数据类型为
double且大小为 8 字节时才能正常工作。movupd指令用于从内存加载数据，因为不保证数据与 16 字节对齐。对于 .Net C++ 项目，我找到了这个构建选项，但对于 C# 项目，我没有找到。
接下来的 FilterPoint 和 FilterPointDiscont 函数执行范围点过滤。区别仅在于 FilterPoint 函数重置不连续标志，而 FilterPointDiscont 函数设置此标志。这样做的原因是，通过算法逻辑可以清楚地知道何时需要不连续标志，何时不需要，因此为了避免不必要的检查并减少执行时间，功能几乎相同的功能被分成了两个函数。
在上图的蓝色和红色方块代表要输出的结果点。可以看到，对于一个结果点，有几个原始点，因此有必要正确处理何时需要不连续标记以及何时不需要。最有趣的是最上面一行，因为即使曲线在一个情况下（左上方的正方形）超出绘图区域，也不需要不连续标记，但对于另一种情况（红色方块）则需要。因此，结果将包含两条要绘制的独立曲线。让我们看看它是如何处理的。
;------------------------------------------------------------------------
;Function performs points simplification by range. If rounded to int point
;is same as previous, point is ignored. Point is outputted to ptr [edi]
;and increments edi register till next use.
;Function increments the counter of outputted points.
;Input:
;Point to examine located in xmm4 register
;Function uses MMX registers mm1, and mm0 to perform rounding and to store
;previously added point values.
;Function set discontinuity flag.
;------------------------------------------------------------------------
FilterPointDiscont PROC PRIVATE
shufpd xmm4,xmm4,1h ;xmm4= Xt,Yt
CVTPD2PI mm1,xmm4
cmp dword ptr[esp+4h], 0h ;check for first point
je __Add
;if not first point perform filtering
movq mm2,mm0
pcmpeqd mm2,mm1 ;if equal all bits of mm2 = 1b
movq mm3,mm2
PSRLQ mm3,20h
pand mm2,mm3
;if both coordinates are equal value will be FFFFFFFF,
;if one of them are not equal value will be 0
movd eax,mm2
cmp eax,0h
jne __Skip
__Add:
call OutputPoint
__Skip:
mov ebx, 01b
ret 0
FilterPointDiscont ENDP
;--------------------------------------------------------------------
;Function outputs result points.
;Input : mm1 contains point to be outputted.
ebx contains discontinuity flag 1 set discontinuity marker 0 do not set
edi register points to output buffer
[esp+8] counter of outputted points
;--------------------------------------------------------------------
OutputPoint PROC PRIVATE
test ebx,01b
jz AddPoint
movq [edi],mm7
add edi,8h
inc dword ptr[esp+8h]
AddPoint:
movq [edi],mm1
add edi,8h
movq mm0,mm1
inc dword ptr[esp+8h]
ret 0
OutputPoint ENDP
要提及的点：
函数
FilterPoint重置不连续标志，上面列出的FilterPointDiscont函数在点过滤完成后设置不连续标志，该标志的值将影响下一个点的输出。函数
OutputPoint检查用作不连续标志的ebx寄存器，如果标志为 1，则在输出当前点之前输出不连续标记。
最后是实现算法的顶层函数的代码。
;---------------------------------------------------------------------
;Main Function to perform transfromations
;Features: Performs in one loop Affine transformation, curve clipping,
;and drawing points simplification
;Iteraction with used functions is performed throght registers.
;Parameters:
; _values -pointer to arrray of double values to be transformed, assume that values are located in maner
; x0,y0,x1,y1,...,xn-1,yn-1,xn,yn
; _count - unsigned integer values of points in _values array (!) count of Points= count of pairs X,Y
; but not the count of elements in _values array.
;_res - pointer to array of integers used as output buffer.
;_m00~_m21 - double numbers, representing elements of transformation matrix.
; Number - index of element(row,column)
;_xMin,xMax,yMin,yMax - double numbers representing clipping area.
;returns - count of points(!) after transformation.
;Notes:
;Because of curve clipping there possible cases then curve goes out of clipping area
;and than comes back and so on. To handle this cases,
;points of discontinuity are marked as 80000000h -32bit integer min value.
;----------------------------------------------------------------------
TransfromClipReducePoints PROC PUBLIC _values :DWORD,
_count :DWORD,
_res :DWORD,
_m00 :QWORD,
_m01 :QWORD,
_m10 :QWORD,
_m11 :QWORD,
_m20 :QWORD,
_m21 :QWORD,
_xMin :QWORD,
_yMin :QWORD,
_xMax :QWORD,
_yMax :QWORD
push ebx
push esi
push edi
;load parameters to corrensponding registers
mov esi, dword ptr[_values] ;input data
mov ecx, dword ptr[_count] ;input points count
mov edi, dword ptr[_res] ;destination pointer
;load transformation matrix to xmm0-xmm2 registers in a way:
; register Low | High
;-------------------------
; xmm0 m00 | m10
; xmm1 m01 | m11
; xmm2 m20 | m21
movsd xmm0,_m00
movhpd xmm0,_m10
movsd xmm1, _m01
movhpd xmm1, _m11
movsd xmm2, _m20
movhpd xmm2, _m21
;loading clipping bounds to xmm6, xmm7 registers in a way
; register Low | High
;-------------------------
; xmm6 xMin | yMin
; xmm7 xMax | yMax
movsd xmm6, _xMin
movhpd xmm6, _yMin
movsd xmm7, _xMax
movhpd xmm7, _yMax
;discontinuity marker loading to mm7 register
mov eax, 80000000h
push eax
push eax
movq mm7, [esp]
add esp,8h
xor eax, eax
push eax ;stack allocation for ouput points count.
xor ebx, ebx ;initialization to store discontinuity flag.
call TransformPoint
;calculating clipping bit code for calculated point and return value in [al] register.
call CalculateBitCode
;in this function during calculations xmm3 contains current point, and xmm5 previos point.
;Store calculated values as previous point
movapd xmm5, xmm3
;registers dl,dh are used to store clip bit codes
; dh - for point located in xmm5 (previos point)
; dl - for point located in xmm3 (current point)
mov dh,al
dec ecx ;because firs point is already processed
StartMainLoop:
call TransformPoint
call CalculateBitCode
mov dl,al ;dl - set dl to current point clip code
;Clipping. Check for two trivial cases
;1st trivial case: both points are outside and almoust on same side
;(same bit is on for both points)
test dl,dh ;test operation performs bitwise AND operation,
;so if same bit in both operators is set result will be not 0(zero)
jnz Continue
;put previous point to xmm4, as is Y,X
movapd xmm4,xmm5 ;because AddPoint use as input xmm4
;2nd trivial case, both points are inside the clipping area
mov al,dh
or al,dl
jnz ClipPrevPoint
;mov al,0h
call FilterPoint
jmp Continue
ClipPrevPoint:
cmp dh,0h
jnz ClipPrev
;mov al,0h
call FilterPoint
jmp ClipCurrentPoint
ClipPrev:
mov al,dh
call CalcClip ;uses eax, xmm3,xmm5 result stored in xmm4
;mov al,0h
call FilterPoint
ClipCurrentPoint:
movapd xmm4,xmm3
cmp dl,0h
jz Continue
mov al,dl
call CalcClip
;mov al,1h
call FilterPointDiscont
Continue:
movapd xmm5,xmm3 ;Store current point to previous
mov dh,dl ;Store current point clip code to previous
loop StartMainLoop
pop eax ;pop return value
emms ;clear MMX state.
pop edi
pop esi
pop ebx
ret 5ch
TransfromClipReducePoints ENDP
如您所见，实现相当线性，并且如您所测，提供的代码以真正超音速执行所需的操作。但是，当然，并非一切都那么美好。提供的代码有一些限制：
- 为了正确工作，有必要实现一个表示二维点的类。例如提供的
PointD类。重要的是，假定该类仅包含两个类型为 double 的成员变量。变量的顺序是 X,Y。 - 因为即使在示例代码中
PointD类也是一个托管类，所以有必要提供获取实例有效地址的功能。 - 根据您的需要，您可以继承提供的
PointD类来扩展其功能，或者创建自己的类，但不要添加成员变量，因为这会导致实例大小发生变化，从而导致数据损坏。
Using the Code
即使主代码是用汇编语言编写的，但该功能可以从托管代码中轻松访问。不仅可以由 Windows Forms 应用程序使用，也可以由 ASP.Net 应用程序使用。下面的代码演示了 Tranformer::TransformDblToInt 函数的实现。当然，您可以创建自己的函数，但需要遵循以下要求才能正常工作：
- 非托管代码接收一个类型为 double 的数组作为输入，并假定数组中的值代表二维点坐标，例如：x0,y0,x1,y1，依此类推。
- 有必要提供足够的输出缓冲区大小。
- 非托管代码的结果是
int类型的数组，值存储方式为：x0,y0,x1,y1，依此类推。 - 由于剪裁和可能的不连续性，结果包含不连续标记，值为 xi 和 yi 被设置为
int::MinValue(80000000h)。
有关更多详细信息，请参阅源代码。
//Function performs points transformation, clipping and simplification
//Function returns the list of System::Drawing::Point arrays.
List<array<Point>^>^ Transformer::TransformDblToInt(
array<PointD>^ dPoints, //input points
double m00,//matrix elements
double m10,
double m20,
double m01,
double m11,
double m21,
double dXMin,//clipping area.
double dYMin,
double dXMax,
double dYMax)
{
//Compare is there enough space in transformation buffer
// multiply 4 because in worst case count of clipped point can be twice more
//than original points, also each point contains 2(two) coordinates(x,y)
if (piConvertBuffer->Length < dPoints->Length*4)
{
piConvertBuffer = gcnew array<int>(dPoints->Length*4);
}
//getting pointers to points array to send to unmanaged code
// result array of integers
pin_ptr<int> piRes = πConvertBuffer[0];
// input array of PointD.
pin_ptr<double> piInput = dPoints[0].GetEffectiveAddress();
//call of unmanaged function.
int iTotal= ::TransfromClipReducePoints(piInput,dPoints->Length, piRes,
m00,m01,m10,m11,m20,m21,dXMin,dYMin,dXMax,dYMax);
//Packing result to send back to call function.
List<array<Point>^>^ oRes =gcnew List<array<Point>^>();
List<Point> oTemp = gcnew List<Point>();
//loop throught all returned points
for (int i=0;i< iTotal; i++)
{
//check for discontinuity marker.
if((int::MinValue == piConvertBuffer[i+i])&&(int::MinValue == piConvertBuffer[i+i+1]))
{
if ( 0 != oTemp.Count )
{
oRes->Add(oTemp.ToArray());
oTemp.Clear();
}
continue;
}
oTemp.Add(Point(piConvertBuffer[i+i],piConvertBuffer[i+i+1]));
}
//check and add last points to result array
if ( 0 != oTemp.Count )
{
oRes->Add(oTemp.ToArray());
oTemp.Clear();
}
oTemp.Clear();
if (0 == oRes->Count)
{
oRes = nullptr;
return nullptr;
}
return oRes;
}关注点
本文的主要观点是证明标准库仍然不完美，并且开发人员应该始终考虑优化。我还试图表明，如果您想充分利用 CPU 的强大功能，即使现在也只有汇编可以帮助您。简而言之，提供的代码涵盖了以下方面：
- 使用
SSE2/SSE3和MMX指令。 - 托管和非托管代码之间的交互。
- 使用 VS.Net 编译汇编文件（不是内联汇编，而是纯 *.asm 文件）。
注释
在演示程序中，仅计算预绘图处理花费的时间。由于不是本文的主题，因此在 System::Drawing::Drawing2d::Matrix::TransformPoints 方法完成变换后，不进行绘图、剪裁和点缩减。仅为 ::TransfromClipReducePoints 函数的结果进行绘图。
演示程序通过在 y 轴上添加随机噪声来生成 y=cos(x) 曲线。
最好在 ::TransfromClipReducePoints 函数返回的结果上执行 Douglas Pecker 简化算法，但这取决于您。我的代码中没有使用此算法，因为该算法不是线性的，并且当曲线点存在很大噪声时可能会导致相当深的递归。
演示程序不执行 CPUID 检查 SSE3 指令支持，因此请确保您的 CPU 支持 SSE3 指令。简单来说，您的 CPU 必须是 P4 或更高版本。对于 AMD CPU - 与 P4 同代的 CPU 不支持 SSE3 指令。
使用的资源：
历史
版本 1.0，发布日期 2009 年 7 月 28 日