jmaths.pas

定点数函数集

	adx := ax1 - ax0;
	ady := ay1 - ay0;
	bdx := bx1 - bx0;
	bdy := by1 - by0;
	TWO := 2 shl 16;

	if (adx = 0) and (bdx = 0) then begin
  // both vertical lines
	 dist := Abs(JFixDiv((ax0 + ax1) - (bx0 + bx1),TWO));
	 result := dist = 0;
  exit;
	end else
 if adx = 0 then begin
  // A  vertical
  xa := JFixDiv((ax0 + ax1),TWO);
  xmb := JFixDiv(bdy,bdx); // slope segment B
  xbb := by0 - JFixMul(bx0,xmb); // y intercept of segment B
  JVxIntersect := xa;
  JVyIntersect := (JFixMul(xmb,JVxIntersect)) + xbb;
	end else
 if bdx = 0 then begin // B vertical
  xb := JFixDiv((bx0 + bx1),TWO);
  xma := JFixDiv(ady,adx);           // slope segment A
  xba := ay0 - (JFixMul(ax0,xma)); // y intercept of segment A
  JVxIntersect := xb;
  JVyIntersect := (JFixMul(xma,JVxIntersect)) + xba;
	end else begin
  xma := JFixDiv(ady,adx); // slope segment A
  xba := ay0 - (JFixMul(ax0,xma)); // y intercept of segment A
  xmb := JFixDiv(bdy,bdx);           // slope segment B
  xbb := by0 - (JFixMul(bx0,xmb)); // y intercept of segment B
  // parallel lines?
  if xma = xmb then begin
		 // Need trig functions
		 dist := Abs(JFixMul((xba-xbb),(JFixCos(JFixArcTan(JFixDiv((xma + xmb),TWO))))));
		 if dist < (1 shl 16) then begin
		  result := true;
    exit;
		 end else begin
		  result := false;
    exit;
		 end;
	 end else begin
		 // Calculate points of intersection
		 // At the intersection of line segment A and B, XA=XB=XINT and YA=YB=YINT
		 if (xma - xmb) = 0 then begin
		  result := false;
    exit;
		 end;
		 JVxIntersect := JFixDiv((xbb-xba),(xma-xmb));
		 JVyIntersect := (JFixMul(xma,JVxIntersect)) + xba;
	 end;
	end;
	// After the point or points of intersection are calculated, each
	// solution must be checked to ensure that the point of intersection lies
	// on line segment A and B.
	minxa := Min(ax0, ax1);
	maxxa := Max(ax0, ax1);
	minya := Min(ay0, ay1);
	maxya := Max(ay0, ay1);
	minxb := Min(bx0, bx1);
	maxxb := Max(bx0, bx1);
	minyb := Min(by0, by1);
	maxyb := Max(by0, by1);
	result := ((JVxIntersect >= minxa) and (JVxIntersect <= maxxa) and
            (JVyIntersect >= minya) and (JVyIntersect <= maxya) and
           	(JVxIntersect >= minxb) and (JVxIntersect <= maxxb) and
            (JVyIntersect >= minyb) and (JVyIntersect <= maxyb));
end;


{ **************************************************************************** }



const
FIX_SHIFT : Longint = 16;

function Float2_Fix(x: double): FIXED;
begin
 Longint(result) :=  Trunc(x * FIX_SCALEF);
end;

function Fix_2Float(Value: FIXED): double;
begin
 result := Longint(Value)* FIX_SCALEF_;
end;

function FIX_FRAC(Value : FIXED): Longint;
begin
 result :=  Longint(Value) and (FIX_SCALE - 1);
end;

function Fix_Scale: Longint;
begin
 result := 1 shl FIX_SHIFT;
end;

function FIX_SCALEF: double;
begin
 result := FIX_SCALE;
end;

function FIX_SCALEF_: double;
begin
 result := 1.0 / FIX_SCALE;
end;

function Fix_Mult(a,b: FIXED): FIXED;
begin
 result := FIXED((Longint(a) * Longint(b)) shl FIX_SHIFT);
end;

function Fix_Div_(a,b: FIXED): FIXED;
begin
 Longint(result) := Trunc((Longint(a) shl FIX_SHIFT) / Longint(b));
end;

function Fix_Sqrt(Num : FIXED): FIXED;  register;// uses Newton’s method to find  num
var
s : Fixed;
i : integer;
begin
 Longint(s) := (Longint(num) + 65536) shr 1;  //divide by 2
 for i := 0 to 5 do  //converge six times
  Longint(s) := (Longint(s) + (Longint(num) div Longint(s))) shr 1;
 result := s;
end;

// ------------------------ geometry function -----------------------------

function Calc2PPixelLength(Pt1,Pt2 : TPoint): Longint;
var
XLen,YLen : Longint;
begin
 XLen := Abs(Pt2.X - Pt1.X);
 YLen := Abs(Pt2.Y - Pt1.Y);
 if XLen >= YLen then result := XLen
                 else result := YLen;
end;

function Calc2PLength(Pt1,Pt2: TFPoint): double;
begin
 result := sqrt(sqr(Pt2.X - Pt1.X) + sqr(Pt2.Y - Pt1.Y));
end;

function Calc2PLength(Pt1,Pt2: TPoint): double;
var
fPt1,fPt2 : TFPoint;
begin
 fPt1.X := Pt1.X; fPt1.Y := Pt1.Y;
 fPt2.X := Pt2.X; fPt2.Y := Pt2.Y;
 result := sqrt(sqr(fPt2.X - fPt1.X) + sqr(fPt2.Y - fPt1.Y));
end;

function Calc2PAngleEx(Pt1,Pt2: TFPoint): double;
begin
 result := Arc2Tangent(Pt2.X - Pt1.X, Pt2.Y - Pt1.Y);
end;

function Calc2PRadian(Pt1,Pt2: TPoint): double;
var
FPt1,FPt2 : TFPoint;
begin
 FPt1.X := Pt1.X; FPt1.Y := Pt1.Y;
 FPt2.X := Pt2.X; FPt2.Y := Pt2.Y;
 result := Arc2Tangent(FPt2.X - FPt1.X, FPt2.Y - FPt1.Y);
end;

function Calc2PAngleEx(Pt1,Pt2: TPoint): double;
var
FPt1,FPt2 : TFPoint;
begin
 FPt1.X := Pt1.X; FPt1.Y := Pt1.Y;
 FPt2.X := Pt2.X; FPt2.Y := Pt2.Y;
 result := Arc2Tangent(FPt2.X - FPt1.X, FPt2.Y - FPt1.Y);
end;

function CalcVerticalPoint(Pt1,Pt2,Pt3: TFPoint): TFPoint;
var
aAng,bAng,P2Ang,P3Ang : double;
hLen,aLen : double;
begin
 hLen := Calc2PLength(Pt1,Pt3);
 P2Ang := Calc2PAngleEx(Pt1,Pt2);
 P3Ang := Calc2PAngleEx(Pt1,Pt3);
 bAng := P2Ang - P3Ang;
 aAng := bAng - PiDivide2;
 aLen := Round(Sin(aAng) * hLen);
 Result.x := Pt1.X - aLen * cos(P2Ang);
 Result.y := Pt1.Y + aLen * sin(P2Ang);
end;

// ---------------------- Contour 等值线（等高线）------------------------------

procedure CalcContour(Memo: TMemo;
                      D: TContourMatrix ; // 2D - Data field
                      ilb,iub, // west - east ilb lower bound
                      jlb,jub : Integer; // iub upper bound north -
                                         // south jlb lower bound jub upper bound
                      x : TContourVector; // coord. vector west - east
                      y : TContourVector; // coord. vector north - south
                      nc: Integer; // nc number of cut levels
                      z : TContourVector); // values of cut levels
const
 im : array [0..3] of Integer = (0,1,1,0);   // coord. cast array west - east
 jm : array [0..3] of Integer = (0,0,1,1);   // coord. cast array north - south
var
m1,m2,m3,deside : Integer;
dmin,dmax,x1,x2,y1,y2 : double;
lcnt,i,j,k,m : integer;
casttab : TContourCastArray;
h : TContourVectorL4D;
sh : TContourVectorL4I;
xh,yh : TContourVectorL4D;
temp1,temp2 : double ;
r : Byte;
 // ------- service xsec west east lin. interpol -------------------------------
 function xsec(p1,p2:Integer):Double;
 begin
  result:=(h[p2]*xh[p1]-h[p1]*xh[p2])/(h[p2]-h[p1]);
 end;
 //------- service ysec north south lin interpol -------------------------------
 function ysec(p1,p2:Integer):Double;
 begin
  result := (h[p2]*yh[p1]-h[p1]*yh[p2])/(h[p2]-h[p1]);
 end;
begin
 if not Assigned(Memo) then exit;
 Memo.Clear;
 // set casting array
 casttab[0,0,0]:=0;casttab[0,0,1]:=0;casttab[0,0,2]:=8;
 casttab[0,1,0]:=0;casttab[0,1,1]:=2;casttab[0,1,2]:=5;
 casttab[0,2,0]:=7;casttab[0,2,1]:=6;casttab[0,2,2]:=9;
 casttab[1,0,0]:=0;casttab[1,0,1]:=3;casttab[1,0,2]:=4;
 casttab[1,1,0]:=1;casttab[1,1,1]:=3;casttab[1,1,2]:=1;
 casttab[1,2,0]:=4;casttab[1,2,1]:=3;casttab[1,2,2]:=0;
 casttab[2,0,0]:=9;casttab[2,0,1]:=6;casttab[2,0,2]:=7;
 casttab[2,1,0]:=5;casttab[2,1,1]:=2;casttab[2,1,2]:=0;
 casttab[2,2,0]:=8;casttab[2,2,1]:=0;casttab[2,2,2]:=0;
 // set line counter
 lcnt:=0;
 for j:=jub-1 downto jlb do begin // over all north - south and + for j
  for i:=ilb To iub-1 do begin // east - west coordinates of datafield + for i
   // set casting bounds from array
   temp1 := min(D[i,j],D[i,j + 1]);
   temp2 := min(D[i + 1,j],D[i + 1,j + 1]);
   dmin  := min(temp1,temp2);
   temp1 := max(D[i,j],D[i,j + 1]);
   temp2 := max(D[i + 1,j],D[i + 1,j + 1]);
   dmax  := max(temp1,temp2);
   If not ((dmax >= z[0]) and (dmin <= z[nc - 1])) then continue;
   // ask horzintal cut avail.    +If dmin && dmax in z[0] .. z[nc-1]
   for k:=0 To nc-1 do begin
    // over all possible cuts ---- +for k
    If (z[k] > dmin) and (z[k] <= dmax) then continue;

    // aks for cut intervall ----- +If z[k] in dmin .. dmax
    for m := 4 downto 0 do begin  // deteriening the cut casts and set the + for m
     If m > 0 then begin // height and coordinate vectors
      h[m] := D[i + im[m - 1],j + jm[m - 1]] - z[k];
      xh[m] := x[i + im[m - 1]];
      yh[m] := y[j + jm[m - 1]];
     end else begin
      h[0] := (h[1] + h[2] + h[3] + h[4]) / 4;
      xh[0] := (x[i] + x[i + 1]) / 2;
      yh[0] := (y[j] + y[j + 1]) / 2;
     end; // If m>0 then else
     If h[m] > 0 then sh[m] := 1 else
     If h[m] < 0 then sh[m] := -1 else sh[m]:= 0;
    end;

    for m := 1 to 4 do begin // set directional casttable
     // Note: at this stage the relative heights of the corners and the
     // centre are in the h array, and the corresponding coordinates are
     // in the xh and yh arrays. The centre of the box is indexed by 0
     // and the 4 corners by 1 to 4 as shown below.
     // Each triangle is then indexed by the parameter m, and the 3
     // vertices of each triangle are indexed by parameters m1,m2,and
     // m3.
     // It is assumed that the centre of the box is always vertex 2
     // though this isimportant only when all 3 vertices lie exactly on
     // the same contour level, in which case only the side of the box
     // is drawn.
     //      AS ANY BODY NOWS IST FROM THE ORIGINAL
     //      vertex 4 +-------------------+ vertex 3
     //               | \               / |
     //               |   \    m-3    /   |
     //               |     \       /     |
     //               |       \   /       |
     //               |  m=2    X   m=2   |       the centre is vertex 0
     //               |       /   \       |
     //               |     /       \     |
     //               |   /    m=1    \   |
     //               | /               \ |
     //      vertex 1 +-------------------+ vertex 2
     //               Scan each triangle in the box
     m1 := m; m2 := 0;
     If not (m = 4) then m3 := m + 1 else m3 := 1;
     deside := casttab[sh[m1] + 1 ,sh[m2] + 1, sh[m3] + 1];
     If not (deside = 0) then begin // ask is there a desition available -------- +If If not(deside=0)
      Case deside of // ------- determin the by desided cast cuts ------------ +Case deside;
        1: begin
            x1 := xh[m1];
            y1 := yh[m1];
            x2 := xh[m2];
            y2 := yh[m2];
           end;
        2: begin
            x1 := xh[m2];
            y1 := yh[m2];
            x2 := xh[m3];
            y2 := yh[m3];
           end;
        3: begin
            x1 := xh[m3];
            y1 := yh[m3];
            x2 := xh[m1];
            y2 := yh[m1];
           end;
        4: begin
            x1 := xh[m1];
            y1 := yh[m1];
            x2 := xsec(m2,m3);
            y2 := ysec(m2,m3);
           end;
        5: begin
            x1 := xh[m2];
            y1 := yh[m2];
            x2 := xsec(m3,m1);
            y2 := ysec(m3,m1);
           end;
        6: begin
            x1 := xh[m3];
            y1 := yh[m3];
            x2 := xsec(m1,m2);
            y2 := ysec(m1,m2);
           end;
        7: begin
            x1 := xsec(m1,m2);
            y1 := ysec(m1,m2);
            x2 := xsec(m2,m3);
            y2 := ysec(m2,m3);
           end;
        8: begin
            x1 := xsec(m2,m3);
            y1 := ysec(m2,m3);
            x2 := xsec(m3,m1);
            y2 := ysec(m3,m1);
           end;
        9: begin
            x1 := xsec(m3,m1);
            y1 := ysec(m3,m1);
            x2 := xsec(m1,m2);
            y2 := ysec(m1,m2);
           end;
      end;
      // ----------do someting with the results ----------------------------
      Memo.Lines.Add(Format('%2.2f %2.2f %2.2f %2.2f %2.2f',[z[k],x1,y1,x2,y2]));
     end;
    end;
   end;
  end;
 end;
end;

procedure CalcContourTest(Memo : TMemo);
const
dimx = 100;  // dimension west - east
dimy = 100;  // dimenstion north west
dimh = 10;   // dimension for contour levels
var
Mat : TContourMatrix;  // 2D - Datafield
scx : TContourVector;  // scaling vector west - east
scy : TContourVector;  // scaling vector north - west
hgt : TContourVector;  // vector for the countur levels
i,j : Integer;  // adress indexes
x,y : Double;   // coord. values
mi,ma : Double; // for minimum & maximum
begin
 if not Assigned(Memo) then exit;
 Memo.Clear;
 setlength(scx,dimx); // create dynamicly the vectors and datafield
 setlength(scy,dimy);
 setlength(hgt,dimh);
 setlength(mat,dimx);
 for i:= 0 to dimx - 1 do Setlength(mat[i],dimy);
 for i:= 0 to dimx - 1 do scx[i] := i * 10; // set scaling vector west - east
 for i:= 0 to dimy - 1 do scy[i] := i * 10; // set scaling vector north - south

 for i := 0 to dimx - 1 do  // set 2d data field
  for j := 0 to dimy - 1 do begin
   x := i - dimx / 2;
   y := j - dimy / 2;
   mat[i,j] := (sin(x / dimx * 4 * pi) * cos(y / dimy * 4 * pi)) +
               (sin(x / dimx * 2 * pi) * cos(y / dimy * 2 * pi)) +
               (sin(x / dimx * 1 * pi) * cos(y / dimy * 1 * pi)) +
               (sin(x / dimx * 0.5 * pi) * cos(y / dimy * 0.5 * pi)) +
               (sin(x / dimx * 0.25 * pi) * cos(y / dimy * 0.25 * pi));
  end;
 mi := 1e16; // Set the minimunm and maximum fof the data field
 ma := -1e16;
 for i := 0 to dimx - 1 do
  for j := 0 to dimy - 1 do begin
   if mat[i,j] < mi then mi := mat[i,j];
   if mat[i,j] > ma then ma := mat[i,j];
  end;
 for i := 0 to dimh - 1 do hgt[i] := mi + i * (ma - mi); // (dimh - 1); // create cut levels
 CalcContour(Memo,mat,0,dimx - 1,0,dimy - 1,scx,scy,dimh,hgt); //call the contour algorithm
end;



end.