-- DynFitImpl.mesa -- Last edited by Michael Plass August 23, 1982 9:18 am -- Michael Plass 30-Sep-81 DIRECTORY JaMFnsDefs USING [GetJaMBreak], Vector, Complex, DynFit, Seq; DynFitImpl: PROGRAM IMPORTS JaMFnsDefs, Vector, Complex EXPORTS DynFit = BEGIN OPEN Seq; PointNumber: TYPE = DynFit.PointNumber; -- These global variables hold the samples NearestZ: PUBLIC PROCEDURE [z: ComplexSequence, v:Complex.Vec] RETURNS [PointNumber] = BEGIN closest:PointNumber _ 0; best:REAL _ 10e+20; n:NAT _ z.length; IF z#NIL THEN FOR i:NAT IN [0..n) DO dd:REAL = Complex.SqrAbs[Complex.Sub[v,z[i]]]; IF dd < best THEN {best _ dd; closest _ i} ENDLOOP; RETURN[closest]; END; LeastSquares: PROCEDURE [numberOfPoints: REAL, sum, sumSqr: Complex.Vec, sumSqrAbs: REAL] RETURNS [REAL] = BEGIN -- finds the sum of the squares of the distances of a set of points to a line, where the line is chosen such that this sum is minimized. The points are specified as complex numbers, and this routine requires their (complex) sum and the sum of their (complex) squares, as well as the sum of the squares of their absolute values. For the best numerical stability, the origin should be chosen as close to the center of mass of the points as is practical. OPEN Complex; RETURN[(sumSqrAbs - SqrAbs[sum]/numberOfPoints - Abs[Sub[sumSqr,Vector.Div[Sqr[sum],numberOfPoints]]])/2] END; BestLine: PUBLIC PROCEDURE [z: ComplexSequence, firstPoint, lastPoint: PointNumber] RETURNS [r: DynFit.Line] = BEGIN -- finds the line that minimizes the sum of the squares of the distances to the points. The equation of the resulting line is given by centerOfMass + t*direction, where t is a real-valued parameter. OPEN Complex; n:NAT _ z.length; f:NAT = IF firstPoint<=lastPoint THEN firstPoint ELSE lastPoint; l:NAT = lastPoint + firstPoint - f; m:NAT _ l - f + 1; a:Complex.Vec _ [0,0]; v:Complex.Vec _ [0,0]; IF m > n THEN m _ n; FOR i:NAT IN [f..l] DO a _ Add[a,z[i]]; ENDLOOP; a _ Vector.Div[a,m]; FOR i:NAT IN [f..l] DO v _ Add[v,Sqr[Sub[z[i],a]]]; ENDLOOP; v _ Vector.Div[v,m]; r.centerOfMass _ a; r.direction _ SqRt[v]; END; tryAll: BOOLEAN _ FALSE; FitSegments: PUBLIC PROCEDURE [z: ComplexSequence, tolerance: REAL] RETURNS [segments: DynFit.Segments, totalBadness: REAL] = BEGIN n:NAT _ z.length; bestPrev:NatSequence = NEW[NatSequenceRec[n]]; segCount:NatSequence = NEW[NatSequenceRec[n]]; totSqrD:RealSequence = NEW[RealSequenceRec[n]]; IF n<2 THEN RETURN [NIL,0]; bestPrev[0] _ 0; segCount[0] _ 0; totSqrD[0] _ 0.0; FOR m:NAT IN [1..n) DO sum,sumSqr: Complex.Vec _ [0,0]; sumSqrAbs: REAL _ 0; b:REAL _ 0; lowbad:REAL _ 1.0E+20; bestl:NAT; s:REAL _ 0; sOpt:REAL; sBound:REAL _ 0.001; FOR l:NAT DECREASING IN [0..m) WHILE tryAll OR s<=sBound DO deltaZ:Complex.Vec = Complex.Sub[z[l],z[m]]; sum _ Complex.Add[sum,deltaZ]; sumSqr _ Complex.Add[sumSqr,Complex.Sqr[deltaZ]]; sumSqrAbs _ sumSqrAbs + Complex.SqrAbs[deltaZ]; s _ LeastSquares[m-l+1,sum,sumSqr,sumSqrAbs]; b _ s + totSqrD[l] + tolerance * (segCount[l] + 1); sBound _ sBound + tolerance; IF b GOTO Quit; epsilon: REAL _ 1.0; Sample: PROCEDURE [v: Complex.Vec] = BEGIN newSamples[len] _ v; len _ len + 1; END; newSamples _ NEW [ComplexSequenceRec[segments.length*2]]; len _ 0; FOR m:NAT DECREASING IN (0..segments.length) DO OPEN Complex; -- uses Vec,Add,Sub,Mul,Abs first:NAT _ segments[m-1]; last:NAT _ segments[m]; line:DynFit.Line _ BestLine[z,first,last]; towardsLast:Vec _ Sub[z[last],line.centerOfMass]; t:Vec = -- tangent vector directed roughly toward last point Vector.Mul[line.direction,Vector.Dot[line.direction,towardsLast]]; s:REAL = Abs[t]; unit:Vec = IF s>0 THEN Vector.Div[t,s] ELSE t; normal:Vec = Mul[unit,[0,1]]; max,min:REAL _ 0; -- extreme components along tangent oldNormalComp,oldTangentComp:REAL _ 0; FOR i:NAT IN [first..last] DO zi:Vec = Sub[z[i],line.centerOfMass]; tangentComp:REAL = Vector.Dot[zi,unit]; normalComp:REAL = Vector.Dot[zi,normal]; candidate:REAL _ 0; IF normalComp < epsilon THEN candidate _ tangentComp ELSE IF oldNormalComp*normalComp < 0 THEN -- find crossing candidate _ (normalComp*oldTangentComp - oldNormalComp*tangentComp) / (normalComp - oldNormalComp); IF candidate < min THEN min _ candidate; IF candidate > max THEN max _ candidate; oldNormalComp _ normalComp; oldTangentComp _ tangentComp; ENDLOOP; Sample[Add[line.centerOfMass,Vector.Mul[unit,max]]]; Sample[Add[line.centerOfMass,Vector.Mul[unit,min]]]; ENDLOOP; EXITS Quit => {} END; END. Michael Plass August 23, 1982 9:18 am. Added check for jam break. Removed CountSamples and SamplesToZ. CountSamples: PUBLIC PROCEDURE [s: CurveDefs.SampleHandle] RETURNS [NAT] = BEGIN i:NAT _ 0; FOR sa:CurveDefs.SampleHandle _ s, sa.next UNTIL sa=NIL DO i_i+1; ENDLOOP; RETURN[i]; END; SamplesToZ: PUBLIC PROCEDURE [s: CurveDefs.SampleHandle] RETURNS [z: ComplexSequence] = BEGIN sa: CurveDefs.SampleHandle; i:PointNumber; n:NAT _ CountSamples[s]; z _ NEW[ComplexSequenceRec[n]]; i _ 0; FOR sa_s, sa.next UNTIL sa=NIL DO z[i] _ sa.xy; i_i+1; ENDLOOP; END; Κz– "Mesa" style˜Iprocš ΟcfœΟk œžœCžœžœžœ žœžœžœ+œΟnœžœž œ%žœžœ!žœžœžœžœžœžœžœžœžœžœ+žœ žœžœžœ žœŸ œž œžœ5žœžœžœžœΖœžœ žœhžœŸœžœž œ=žœžœΘœžœ žœžœžœžœ žœžœ"žœ?žœžœ žœžœžœžœžœžœžœžœžœ$žœFžœ žœžœŸ œžœž œ!žœžœ+žœžœžœ%žœ-žœ-žœžœžœžœžœ=žœžœžœžœ1žœ žœžœžœžœ žœ žœ žœžœž œžœžœžœ žœΖžœ žœ$žœfžœžœ žœžœCœžœžœžœ#žœ-žœžœž œžœžœ+žœžœžœŸœžœž œ2žœ$žœžœ1œžœžœžœžœ Ÿœž œžœ+žœžœ5žœžœž œžœžœžœ œžœžœ{5œJžœžœžœžœ.žœ$œžœžœžœžœžœ9žœ&žœ'žœ žœžœžœžœžœœ{žœžœžœžœUžœržœžœ žœžœ%œ%Ÿ œŸ œŸ œžœž œžœžœžœžœžœ(žœžœžœ žœžœžœŸ œžœž œ#žœžœ1žœžœ"žœžœžœžœžœžœ˜Χ.—…—ZΪ