delphi FMX图像的直方图统计增强

//直方图统计增强 3x3
procedure tongji(b : TBitmap);
var
   b_read : TBitmap;
   //增强的参数  c
   k0,k1,k2,e : single ;
   //均值  和  方差  c
   mg,ug,m3,u3 : single ;

   x, y , I: Integer;
   wdata , rdata : TBitmapData ;
   p0,p,p1,pw: PByteArray;
   temp : array [0..8] of Byte;
   //次数   c
   n : array [0..255] of Integer;
   //概率   p
   np : array [0..255] of single;
begin
   b_read := TBitmap.Create;
   b_read.Assign(b);
   //初始化参数   全图的均值 和方差   c
   k0 := 0.4;
   k1 := 0.001;
   k2 := 0.4;
   e := 5.0;
   mg := 0;
   ug := 0;
   //初始化数组 c
   for I := 0 to 255 do
   begin
     n[i] := 0;
   end;

   if  b.Map( TMapAccess.ReadWrite, wdata)and b_read.Map( TMapAccess.Read,rdata)  then
   begin
        //统计出现的次数  t
        for y := 0 to rdata.Height - 1 do
        begin
            p := rdata.GetScanline(y);
            for x := 0 to rdata.Width - 1 do
            begin
               n[p[x*4]] := n[p[x*4]]+1;
            end;
        end;
        //计算概率    p
        for I := 0 to 255 do
        begin
            np[i] := n[i]/b.Width/b.Height;
        end;
        //算出均值    j
        for I := 0 to 255 do
        begin
            mg := mg + i*np[i];
        end;
        //算出方差    f
        for I := 0 to 255 do
        begin
            ug := ug + power((i-mg),2)*np[i];
        end;
        //开始遍历图片  小窗口 3x3
        for y := 1 to rdata.Height - 2 do
        begin
            p0:= rdata.GetScanline(y-1);
            p := rdata.GetScanline(y);
            p1:= rdata.GetScanline(y+1);
            pw := wdata.GetScanline(y);
            //3x3
            for x := 1 to rdata.Width - 2 do
            begin
               //初始化局部均值和方差  c
               m3 := 0;
               u3 := 0;
               temp[0] := p0[x*4-4];
               temp[1] := p0[x*4];
               temp[2] := p0[x*4+4];
               temp[3] := p[x*4-4];
               temp[4] := p[x*4];
               temp[5] := p[x*4+4];
               temp[6] := p1[x*4-4];
               temp[7] := p1[x*4];
               temp[8] := p1[x*4+4];
               //局部均值    j
               for I := 0 to 8 do
               begin
                   m3 := m3+temp[i]/9;
               end;
               //局部方差    f
               for I := 0 to 8 do
               begin
                   u3 := u3+power((temp[i]-m3),2)/9;
               end;
              //判断是否增强  p
               if (m3<=k0*mg)and(k1*ug<=u3)and(u3<=k2*ug) then
               begin
                 pw[x*4] :=  round(p[x*4]*e);
                 pw[x*4+1] := pw[x*4];
                 pw[x*4+2] := pw[x*4];
               end;
            end;
        end;

        b.Unmap(wdata);
        b_read.Unmap(rdata);
   end;

end;