**Reading Tips:**

**《****Delphi Image Processing****The series focuses on efficiency. The general code is Pascal, and the core code is BaSm.**

**《****C ++ Image Processing****The series focuses on code clarity and readability, and all uses C ++ code.**

**Make sure that the two items are consistent and can be compared with each other.**

**The code in this article must include the imagedata. Pas unit in "Delphi Image Processing-data type and public process.**

**Note: The image Gaussian fuzzy processing code is modified most often. Although this modification does not change the algorithm, the processing process has been modified. This alone can increase the processing speed by more than 40% on the original basis. At the same time, the SSE floating-point operation is used to replace the fixed-point operation of the original assembly. In order to facilitate the comparison, and do not want to destroy the original code, the original code is slightly modified and retained. The modified code is appended to the article.**

In my article "Delphi Image Processing-image convolution", I once introduced the use of common image convolution processes to perform Gaussian blur on the image. The processing effect is not bad, it feels okay to process small images, but it is still too slow to process large images, on my P4 2.8 GB and 1 GB memory machines, perform Gaussian fuzzy q = 3, r = 5 on 10 million pixel images, excluding image loading and pre-data conversion, the processing speed has never been significantly improved after several modifications. The main reason is the general convolution process: the convolution template obtained with r = 5 is 11x11 pixels, and each pixel has four components (32-bit argb ), each pixel must be subject to 11*11*4 = 484 multiplication, 484 addition, and 4 Division. Finally, it must be subject to determining whether four components are super-bounded, it's hard to think about it! If you do not use BaSm to Process Code with fixed points, the processing speed is even harder to imagine.

I searched for image Gaussian Blur optimization algorithms multiple times on the Internet. Many algorithms and processing methods, including code optimization, were not as good as my Gaussian fuzzy processing process. I was very disappointed. When searching for other materials the day before yesterday, a foreign website found that the Gaussian Blur Processing Method for images seemed to be different from conventional algorithms, but there was no detailed information (because they do not understand foreign languages, I seldom access websites outside China, but I still read some formulas and pseudo-code.) After repeated efforts, I can summarize the processing process as follows:

1. Calculate a size + 1 long Gaussian distribution weight data using the given Q and Length Size:

// Calculate the initial data for (I =-radius; I <= radius; I ++) {x = I/Q; weights [I + radius] = exp (-x * X/2)} // sum = 0for (I =-radius; I <= radius; I ++) {sum + = weights [I + radius]} // Data normalization, that is, the sum of the data after normalization is equal to 1for (I =-radius; I <= radius; I ++) {weights [I + radius]/= sum}

2. Use weights to perform vertical fuzzy operations on the original image, that is, taking the pixel (x, y) as the center, right (X, Y-radius) and (X, Y + radius) after the pixels are multiplied by the values of weights, the new pixels are obtained and written to the corresponding points on a temporary image (because the data is normalized, product and Division operations are not required );

3. Use weights to perform horizontal Fuzzy Operations on temporary images, that is, taking the pixel (x, y) as the center, to (x-radius, Y) and (x + radius, Y) after the pixels are multiplied by weights, the new pixels are obtained and written to the corresponding points on the target image.

The process ends.

Since the above processing process only performs a "Ten" operation on each pixel of the image, the operation on each pixel point is greatly reduced, and the greater the fuzzy length, the more reduced. As mentioned above, the Q = 3 and r = 5 Fuzzy Operations only require 11*2*4 = 88 multiplication and 88 addition operations.

I still use BaSm to calculate the number of points according to the above process. The improved Gaussian fuzzy Process Code is as follows:

procedure CrossBlur(var Dest: TImageData; const Source: TImageData; Weights: Pointer; Size: Integer);var height, srcStride: Integer; _weights: Pointer; dstOffset, srcOffset: Integer; reds, greens, blues: Integer;asm push esi push edi push ebx mov _Weights, ecx mov ecx, [edx].TImageData.Stride mov srcStride, ecx call _SetCopyRegs mov height, edx mov dstOffset, ebx push esi push edi push edx push ecx push eax // blur col add ecx, Size // width = Source.Width dec ecx mov edi, _weights // edi = weights@@cyLoop: push ecx@@cxLoop: push ecx push esi push edi xor ebx, ebx mov reds, ebx mov greens, ebx mov blues, ebx mov ecx, Size@@cblurLoop: movzx eax, [esi].TARGBQuad.Blue movzx edx, [esi].TARGBQuad.Green imul eax, [edi] imul edx, [edi] add blues, eax add greens, edx movzx eax, [esi].TARGBQuad.Red movzx edx, [esi].TARGBQuad.Alpha imul eax, [edi] imul edx, [edi] add reds, eax add ebx, edx add edi, 4 add esi, srcStride loop @@cblurLoop pop edi pop esi mov eax, blues mov edx, greens mov ecx, reds shr eax, 16 shr edx, 16 shr ecx, 16 shr ebx, 16 mov [esi].TARGBQuad.Blue, al mov [esi].TARGBQuad.Green, dl mov [esi].TARGBQuad.Red, cl mov [esi].TARGBQuad.Alpha, bl add esi, 4 pop ecx loop @@cxLoop pop ecx dec height jnz @@cyLoop pop srcOffset pop ecx pop height pop edi pop esi // blur row@@ryLoop: push ecx@@rxLoop: push ecx push esi push edi xor ebx, ebx mov reds, ebx mov greens, ebx mov blues, ebx mov ecx, Size mov edi, _weights@@rblurLoop: movzx eax, [esi].TARGBQuad.Blue movzx edx, [esi].TARGBQuad.Green imul eax, [edi] imul edx, [edi] add blues, eax add greens, edx movzx eax, [esi].TARGBQuad.Red movzx edx, [esi].TARGBQuad.Alpha imul eax, [edi] imul edx, [edi] add reds, eax add ebx, edx add edi, 4 add esi, 4 loop @@rblurLoop pop edi pop esi mov eax, blues mov edx, greens mov ecx, reds shr eax, 16 shr edx, 16 shr ecx, 16 shr ebx, 16 mov [edi].TARGBQuad.Blue, al mov [edi].TARGBQuad.Green, dl mov [edi].TARGBQuad.Red, cl mov [edi].TARGBQuad.Alpha, bl add esi, 4 add edi, 4 pop ecx loop @@rxLoop add esi, srcOffset add edi, dstOffset pop ecx dec height jnz @@ryLoop pop ebx pop edi pop esiend;procedure ImageGaussiabBlur(var Data: TImageData; Q: double; Radius: Integer);var src: TImageData; fweights: array of Single; weights: array of Integer; i, size: Integer; fx: Double;begin if Radius <= 0 then begin if Abs(Q) < 1.0 then Radius := 1 else Radius := Round(Abs(Q)) + 2; end; size := Radius shl 1 + 1; SetLength(fweights, size); for i := 1 to Radius do begin fx := i / Q; fweights[Radius + i] := exp(-fx * fx / 2); fweights[Radius - i] := fweights[Radius + i]; end; fweights[Radius] := 1.0; fx := 0.0; for i := 0 to size - 1 do fx := fx + fweights[i]; SetLength(weights, size); for i := 0 to size - 1 do weights[i] := Round(fweights[i] / fx * 65536.0); SetLength(fweights, 0); src := _GetExpandData(Data, Radius); CrossBlur(Data, src, weights, size); FreeImageData(src);end;

The improved Gaussian fuzzy processing process is used on my machine to perform Gaussian fuzzy q = 3, r = 5 on 10 million pixel images, excluding image loading and pre-data conversion, the processing speed was significantly improved with a time consumption of 1390 Ms. I process the image pixels in 32-bit argb color. If I change the color to 24-bit RGB, the time consumption can be reduced. However, RGB colors cannot process 32-bit pixel images such as PNG.

What is the result of Gaussian fuzzy processing using ten words instead of template convolution? Please refer to the following simple example code and processing:

Sample Code:

procedure TForm1.Button3Click(Sender: TObject);var bmp: TGpBitmap; g: TGpGraphics; data: TImageData;begin bmp := TGpBitmap.Create('..\media\56-3.jpg'); g := TGpGraphics.Create(Canvas.Handle); g.DrawImage(bmp, 0, 0); data := LockGpBitmap(bmp); ImageGaussiabBlur(Data, 3, 6); UnlockGpBitmap(bmp, data); g.DrawImage(bmp, data.Width, 0); g.Free; bmp.Free;end;

Process the source image:

Comparison between processing effect and Photoshop Gaussian fuzzy processing:

On the top left is the Photoshop radius 3.0 Gaussian blur, and on the top right is the process q = 3.0, r = 6 Gaussian blur.

The lower left is the 5.0 Gaussian blur of Photoshop radius, and the lower right is the process q = 5.0, r = 9 Gaussian blur.

How about it? The effect is good!

Unfortunately, I failed to find the method for automatically calculating the Blur radius according to Q, so two parameters Q and radius are provided in the processing process.

**The following is the SSE code after this modification. The principle and algorithm are the same as above, but the processing method is somewhat different: Because the upper and lower sides of the Gaussian fuzzy matrix are symmetric, the radius point is the center, add the upper and lower symmetric rows (during column processing) or the left and right symmetric columns (during row processing) and multiply them with the Gaussian distribution weight data. In this way, except for the central row (column, only the previous 50% processing is required.**

procedure CrossBlur(var Dest: TImageData; const Source: TImageData; Weights: Pointer; Radius: Integer);var height, srcStride: Integer; dstOffset, srcOffset: Integer;asm push esi push edi push ebx push ecx mov ecx, [edx].TImageData.Stride mov srcStride, ecx call _SetCopyRegs mov height, edx mov srcOffset, eax mov dstOffset, ebx pop ebx pxor xmm7, xmm7 push esi // pst = Source.Scan0 push edi push edx push ecx // blur col mov eax, srcStride mov edx, eax shr edx, 2 // width = Source.Width mov edi, Radius shl edi, 1 imul edi, eax add edi, esi // psb = pst + Radius * 2 * Source.Stride@@cyLoop: push edx@@cxLoop: push esi push edi push ebx mov ecx, Radius pxor xmm0, xmm0 // sum = 0@@cblurLoop: movd xmm1, [esi] // for (i = 0; i < Radius; i ++) movd xmm2, [edi] // { punpcklbw xmm1, xmm7 punpcklbw xmm2, xmm7 paddw xmm1, xmm2 // ps = pst + psb punpcklwd xmm1, xmm7 cvtdq2ps xmm1, xmm1 // pfs (flaot * 4) = ps (int * 4) mulps xmm1, [ebx] // pfs *= Weights[i] addps xmm0, xmm1 // sum += pfs add ebx, 16 add esi, eax // pst += Source.Stride sub edi, eax // psb -= Source.Stride loop @@cblurLoop // } movd xmm1, [esi] punpcklbw xmm1, xmm7 punpcklwd xmm1, xmm7 cvtdq2ps xmm1, xmm1 // pfs (flaot * 4) = pst (int * 4) mulps xmm1, [ebx] // pfs *= Weights[Radius] addps xmm0, xmm1 // sum += pfs pop ebx pop edi pop esi cvtps2dq xmm0, xmm0 // ps (int * 4) = sum (flaot * 4) packssdw xmm0, xmm7 packuswb xmm0, xmm7 movd [esi], xmm0 // pst (byte * 4) = ps (int * 4) pask add esi, 4 add edi, 4 dec edx jnz @@cxLoop pop edx dec height jnz @@cyLoop pop edx pop height pop edi // pd = Dest.Scan0 pop esi // psl = pst mov eax, Radius shl eax, 1+2 add eax, esi // psr = psl + Radius * 2 // blur row@@ryLoop: push edx // width = Dest.Width@@rxLoop: push esi push ebx push eax mov ecx, Radius pxor xmm0, xmm0 // sum = 0@@rblurLoop: movd xmm1, [esi] // for (i = 0; i < Radius; i ++) movd xmm2, [eax] // { punpcklbw xmm1, xmm7 punpcklbw xmm2, xmm7 paddw xmm1, xmm2 // ps = psl + psr punpcklwd xmm1, xmm7 cvtdq2ps xmm1, xmm1 // pfs (flaot * 4) = ps (int * 4) mulps xmm1, [ebx] // pfs *= Weights[i] addps xmm0, xmm1 // sum += pfs add ebx, 16 add esi, 4 // psl ++ sub eax, 4 // psr -- loop @@rblurLoop // } movd xmm1, [esi] punpcklbw xmm1, xmm7 punpcklwd xmm1, xmm7 cvtdq2ps xmm1, xmm1 // pfs (flaot * 4) = psl (int * 4) mulps xmm1, [ebx] // pfs *= Weights[Radius] addps xmm0, xmm1 // sum += pfs cvtps2dq xmm0, xmm0 // ps (int * 4) = sum (flaot * 4) packssdw xmm0, xmm7 packuswb xmm0, xmm7 movd [edi], xmm0 // pd (byte * 4) = ps (int * 4) pask pop eax pop ebx pop esi add eax, 4 add esi, 4 add edi, 4 dec edx jnz @@rxLoop add eax, srcOffset add esi, srcOffset add edi, dstOffset pop edx dec height jnz @@ryLoop pop ebx pop edi pop esiend;// --> st x// <-- st e**x = 2**(x*log2(e))function _Expon: Extended;asm fldl2e // y = x*log2e fmul fld st(0) // i = round(y) frndint fsub st(1), st // f = y - i fxch st(1) // z = 2**f f2xm1 fld1 fadd fscale // result = z * 2**i fstp st(1)end;function GetWeights(var Buffer, Weights: Pointer; Q: Single; Radius: Integer): Integer;const _fcd1: Single = 0.1; _fc1: Single = 1.0; _fc2: Single = 2.0; _fc250: Single = 250.0; _fc255: Single = 255.0;var R: Integer; v, QQ2: double;asm mov R, ecx mov ecx, eax fld Q fabs fcom _fcd1 fstsw ax sahf jae @@1 fld _fcd1 fstp st(1) // if (Q < 0.1) Q = 0.1 jmp @@2@@1: fcom _fc250 fstsw ax sahf jbe @@2 fld _fc250 fstp st(1) // if (Q > 250) Q = 250@@2: fst Q fmul Q fmul _fc2 fstp QQ2 // QQ2 = 2 * Q * Q fwait mov eax, R test eax, eax jg @@10 push eax // if (radius <= 0) fld1 // { fadd Q // radius = Abs(Q) + 1 fistp [esp].Integer fwait pop eax@@testRadius: // while (TRUE) mov R, eax // { fldz // sum = 0@@testLoop: // for (R = radius; R > 0; R ++) fild R // { fld st(0) fmulp st(1), st fdiv QQ2 fchs call _Expon // tmp = Exp(-(R * R) / (2.0 * Q * Q)); cmp R, eax jne @@3 fst v // if (R == radius) v = tmp@@3: faddp st(1), st(0) // sum += tmp dec R jnz @@testLoop // } fmul _fc2 // sum *= 2 fadd _fc1 // sum += 1 fdivr v fmul _fc255 fistp R cmp R, 0 je @@4 // if ((INT)(v / sum * 255 + 0.5) = 0) break inc eax // radius ++ jmp @@testRadius // }@@4: dec eax jnz @@5 inc eax@@5: mov R, eax // }@@10: inc eax shl eax, 4 add eax, 12 push edx push ecx mov edx, eax mov eax, GHND call GlobalAllocPtr pop ecx pop edx test eax, eax jz @@Exit mov [ecx], eax // buffer = GlobalAllocPtr(GHND, (Radius + 1) * 16 + 12) add eax, 12 and eax, -16 mov [edx], eax // weights = ((char* )buffer + 12) & 0xfffffff0 mov ecx, R // ecx = radius mov edx, eax // edx = weights fldz // for (i = radius, sum = 0; i > 0; i --)@@clacLoop: // { fild R fld st(0) fmulp st(1), st fdiv QQ2 fchs call _Expon fstp [edx].Double // weights[i] = Expon(-(i * i) / (2 * Q * Q)) fadd [edx].Double // sum += weights[i] add edx, 16 dec R jnz @@clacLoop // } fmul _fc2 // sum *= 2 fld1 fstp [edx].Double // weights[radius] = 1 fadd [edx].Double // sum += weights[radius] push ecx inc ecx@@divLoop: // for (i = 0; i <= Radius; i ++) fld st(0) // weights[i] = Round(weights[i] / sum) fdivr [eax].Double fst [eax].Single fst [eax+4].Single fst [eax+8].Single fstp [eax+12].Single add eax, 16 loop @@divLoop ffree st(0) fwait pop eax // return Radius@@Exit:end;procedure ImageGaussiabBlur(var Data: TImageData; Q: Single; Radius: Integer);var Buffer, Weights: Pointer; src: TImageData;begin Radius := GetWeights(Buffer, Weights, Q, Radius); if Radius = 0 then Exit; if Data.AlphaFlag then ArgbConvertPArgb(Data); src := _GetExpandData(Data, Radius); CrossBlur(Data, src, Weights, Radius); FreeImageData(src); GlobalFreePtr(Buffer); if Data.AlphaFlag then PArgbConvertArgb(Data);end;

**For details about the use of GDI + units and descriptions in the "Delphi image processing" series, see the article 《****GDI + for VCL basics-GDI + and VCL****.**

**Due to limited levels, errors are inevitable. Correction and guidance are welcome. Email Address:****Maozefa@hotmail.com**

**Here, you can access "Delphi Image Processing-Article Index".**