MODULE PZFilter EXPORTS Main; (* Filter a curve with "a-posteriori" parametrization. *) (* Last edited on 2000-04-29 10:23:06 by stolfi *) (* This program takes an arbitrary curve "c" and outputs several versions "d[0],d[1],..." of it, geometrically filtered with Gaussian filters of scale "lambda", "2*lambda", "4*lambda", etc. In signal-processing contexts, a Gaussian filter with scale "sigma" is a filter whose kernel is a Gaussian distribution having mean zero and standard deviation "sigma", i.e. | G_sigma(t) = sigma/sqrt(2Pi) exp(-t^2/(2*sigma^2)) `Geometric filtering' a curve "c" with a filter "F" means applying "F" to "c", after reparametrizing "c" in such a way that the filtered curve has unit velocity. Actually, each output curve "d[i]" after "d[0]" is obtained by filtering the preceding curve "d[i-1]" with a Gaussian filter whose scale is computed so that the result will be roughly equivalent to filtering "c" to scale "2^i*lambda". The purpose of this optimization is to reduce the amount of work from "O(n log n)" to "O(n)", where "n" is the number of samples in curve "c[0]". *) IMPORT ParseParams, Process, FileRd, Wr, FileWr, Math; IMPORT OSError, Thread, Stdio, Fmt; IMPORT LR3; IMPORT PZCurve, PZLR3Chain, PZLRChain; FROM Stdio IMPORT stderr; FROM PZTypes IMPORT LONG, LONGS, NAT, BOOL; <* FATAL Wr.Failure, Thread.Alerted, OSError.E *> TYPE Options = RECORD input: TEXT; (* Input float chain file (without ".flc"). *) output: TEXT; (* Output float/curv chain file name (without ".flc"). *) lambdaMin: LONG; (* Finest filtering scale length. *) lambdaRatio: LONG; (* Ratio between successive scales. *) nScales: NAT; (* Number of scales (default +oo). *) tolLength: LONG; (* Max change in length (rel. lambda) per iteration. *) tolStep: LONG; (* Max change in interval size (rel lambda) per iter. *) maxIter: NAT; (* Maximum number of iterations (0 for parametric). *) nRepStages: NAT; (* Number of reparametrization stages. *) writeRaw: BOOL; (* TRUE to write raw approximation to input contour *) writeStep: NAT; (* Write the curve every this many filtering scales. *) maxStep: LONG; (* Maximum sampling step for output (default +oo). *) verbose: BOOL; (* TRUE to make noises while working *) END; PROCEDURE Main() = VAR cOld, cNew: PZCurve.T; lambdaOld, lambdaNew: LONG; totIter: NAT; iScale: NAT := 0; CONST MaxRelativeLambda = 0.45d0; (* Stop when lambda/L is grater than this. *) BEGIN WITH o = GetOptions(), inputFileName = o.input & ".flc", rp = ReadCurve(inputFileName, o.verbose), c = rp.curve, cCmt = rp.cmt & "PZFilter:\n" & " input file = " & inputFileName & "\n" DO IF o.writeRaw THEN WriteRawCurve(o, cCmt, c); END; cOld := c; lambdaOld := 0.0d0; lambdaNew := o.lambdaMin; iScale := 0; (* Filter for successive scales: *) LOOP IF lambdaNew > MaxRelativeLambda * cOld.tPeriod OR iScale >= o.nScales THEN EXIT END; cNew := GeometricFilter( cOld, lambdaOld, lambdaNew, o.nRepStages, o.tolLength, o.tolStep, o.maxIter, totIter, o.verbose ); IF cNew.tPeriod <= lambdaNew THEN EXIT END; IF iScale = o.nScales OR iScale MOD o.writeStep = 0 THEN WriteSampledCurve(o, cCmt, cNew, lambdaNew, totIter) END; cOld := cNew; lambdaOld := lambdaNew; lambdaNew := lambdaNew * o.lambdaRatio; INC(iScale); END; Wr.PutText(stderr, "maximum lambda = " & FLR(lambdaOld,8,5) & "\n" ); END END Main; PROCEDURE GeometricFilter( VAR c: PZCurve.T; lambdaOld: LONG; lambdaNew: LONG; nRepStages: NAT; tolLength: LONG; tolStep: LONG; maxIter: NAT; VAR totIter: NAT; (* Total iterations performed for all stages *) verbose: BOOL; ): PZCurve.T = (* Takes a curve "c" that has already been geometrically filtered to scale "lambdaOld" and returns a copy "d" of "c" processed with a geometric Gaussian filter of scale "lambdaNew > lambdaOld", sampled at "n" uniformly spaced points, where "n" is chosen according to the Nyquist criterion. Upon entry, assumes "c" is parametrized by arclength of the "lambdaOld"-filtered version. Returns the "lambdaNew"-filtered curve "d" with the natural parametrization, and reparametrizes "c" with the same parametrization. *) CONST MaxIntermStageIter = 3; (* Max iterations in intermediate stages. *) VAR d: PZCurve.T; change: LONG; stageIter: NAT; (* Iteration count per stage. *) maxStageIter: NAT; (* Max iterations in current stage. *) maxLength: LONG; (* Max length of "lambdaNew"-filtered curve. *) ra, ta, tb: REF LONGS := NIL; (* Work vectors for FilterStep. *) BEGIN maxLength := c.tPeriod; totIter := 0; FOR k := 1 TO nRepStages DO stageIter := 0; IF k < nRepStages THEN maxStageIter := MaxIntermStageIter ELSE maxStageIter := LAST(NAT) END; WITH exp = FLOAT(nRepStages - k, LONG)/FLOAT(MAX(1, nRepStages - 1), LONG), lambdaStg = lambdaNew * Math.pow(2.0d0, exp), minSamples = CEILING(4.0d0 * maxLength/lambdaStg) (* exp = FLOAT(k, LONG)/FLOAT(nRepStages, LONG), *) (* lambdaStg = lambdaOld * Math.pow(lambdaNew/lambdaOld, exp) *) DO LOOP INC(totIter); INC(stageIter); IF verbose THEN Wr.PutText(stderr, "\n") END; WITH (* Compute the parameter "lambdaInc" of a convolution that applied to "c" has the same effect as applying convolution "lambdaStg" to the original curve "p": *) lambdaInc = Math.sqrt(lambdaStg*lambdaStg - lambdaOld*lambdaOld), TOld = c.tPeriod + 0.0d0 DO IF verbose THEN Wr.PutText(stderr, "iteration " & Fmt.Int(totIter) & " "); Wr.PutText(stderr, "lambdaStg = " & FLR(lambdaStg,8,5) & " "); Wr.PutText(stderr, "lambdaInc = " & FLR(lambdaInc,8,5) & "\n"); END; change := FilterStep(c, lambdaInc, d, minSamples, verbose, ra,ta,tb); <* ASSERT d.tPeriod < lambdaInc OR c.tPeriod = d.tPeriod *> IF d.tPeriod <= lambdaInc OR d.m < 2 THEN Wr.PutText(stderr, "! curve collapsed, stage abandoned !\n"); EXIT END; IF ABS(change) <= tolStep*lambdaStg AND ABS(TOld - d.tPeriod) <= tolLength*lambdaStg THEN Wr.PutText(stderr, "! stage converged !\n"); EXIT END; IF totIter >= maxIter THEN IF maxIter > 1 THEN Wr.PutText(stderr, "! too many iterations, stopped !\n") END; EXIT END; IF stageIter >= maxStageIter THEN Wr.PutText(stderr, "! too many iterations, interm stage stopped !\n"); EXIT END; END; END; END; END; RETURN d END GeometricFilter; PROCEDURE FilterStep( VAR cOrg: PZCurve.T; (* Original curve, possibly reparametrized. *) lambda: LONG; (* Standard deviation of convolution kernel. *) VAR cFil: PZCurve.T; (* Filtered curve (preallocated). *) minSamples: NAT; (* Number of samples by Nyquist's limit *) verbose: BOOL; (* TRUE to mumble while working. *) (* WORK*) ra, ta, tb: REF LONGS; (* Working vector *) ): LONG = (* Convolves the curve "cOrg" with a Gaussian of unit area and characteristic duration "lambda", i.e. "exp(-(t/lambda)^2/2)/A" where "A" is the unit-area normalization term. Samples the filtered curve at about "minSamples" equally spaced times, and stores the result in "cFil". Then reparametrizes "cOrg" so that the time lapse between two samples of "cOrg" is the length of "cFil" between the points corresponding to those samples. (This step is omitted if the filtered curve "cFil" has length less than lambda.) Reallocates "cFil" if it hasn't been initialized, or has too few or too many points. Returns the largest (in absolute value) signed change in the spacing between consecutive sampling times of "cOrg". *) VAR change: LONG; BEGIN IF cFil = NIL OR cFil.m < minSamples OR cFil.m > (12*minSamples) DIV 10 THEN cFil := PZCurve.New(minSamples); IF verbose THEN Wr.PutText(stderr, "resizing to " & Fmt.Int(cFil.m) & " samples\n") END END; IF verbose THEN Wr.PutText(stderr, "period = " & FLR(cOrg.tPeriod,8,5) & " "); Wr.PutText(stderr, "filtering...\n"); END; cFil.setToFiltered(cOrg, lambda); <* ASSERT cFil.tPeriod + 0.0d0 = cFil.tPeriod *> (* To catch NaN *) change := Reparametrize(cOrg, cFil, lambda, verbose, ra, ta, tb); <* ASSERT cFil.tPeriod < lambda OR cOrg.tPeriod = cFil.tPeriod *> RETURN change END FilterStep; <*UNUSED*> PROCEDURE Debug(msg: TEXT; READONLY x: LONGS) = <* FATAL Wr.Failure, Thread.Alerted *> BEGIN Wr.PutText(stderr, msg); FOR i := 0 TO LAST(x) DO IF i > 0 AND i MOD 5 = 0 THEN Wr.PutText(stderr, "\n") ELSE Wr.PutText(stderr, " ") END; Wr.PutText(stderr, Fmt.Pad(Fmt.LongReal(x[i], prec := 8, style := Fmt.Style.Fix), 12) ) END; Wr.PutText(stderr, "\n"); END Debug; PROCEDURE Reparametrize( VAR cOrg: PZCurve.T; VAR cFil: PZCurve.T; minLength: LONG; verbose: BOOL; (* WORK *) VAR ra, ta, tb: REF LONGS; ): LONG = (* Given a curve "cOrg", and a parametrically filtered and resampled version "cFil" of "cOrg", this procedure reparametrizes both by the arc length of "cFil". More precisely, it sets the time "cFil.t[j]" to be the length of the curve "cFil(t)" from sample "cFil.p[0]" to sample "cFil.p[j]". Then it sets the time "cOrg.t[i]" associated with each sample "cOrg.p[i]" to the value "t" such that "cFil(t)" is the point of "cFil" that corresponds to "cOrg.p[i]". The routine then computes the time intervals between the consecutive samples of "cOrg", and compares those intervals with the corresponding intervals of old times. The returned value is the maximum absolute change in those intervals. As a special case, if the length of "cFil" is less than "minLength", the curve "cOrg" is not reparametrized, and the procedure returns "minLength". *) VAR maxChange, minStep, maxStep: LONG; tao, tan: LONG; maxChangeIndex: NAT; BEGIN <* ASSERT cFil.tPeriod = cOrg.tPeriod *> <* ASSERT cFil.tPeriod + 0.0d0 = cFil.tPeriod *> (* To catch NaNs*) <* ASSERT cOrg.tPeriod + 0.0d0 = cOrg.tPeriod *> (* To catch NaNs*) EnsureMinSize(ra, cOrg.m); EnsureMinSize(ta, cFil.m); EnsureMinSize(tb, cOrg.m); WITH tOrg = cOrg.t^, tFil = cFil.t^, TOld = cOrg.tPeriod + 0.0d0, tFilOld = SUBARRAY(ta^, 0, cFil.m), (* Original times of filtered samples. *) tOrgNew = SUBARRAY(tb^, 0, cOrg.m), (* New times of original samples. *) arg = SUBARRAY(ra^, 0, cOrg.m) (* Work array for reparametrization. *) DO (* Save the original times of filtered samples: *) tFilOld := tFil; (* Set the times of filtered samples to be length since "cFil.p[0]": *) cFil.uniformize(); WITH TNew = cFil.tPeriod DO IF verbose THEN Wr.PutText(stderr, "filtered curve length = " & FLR(TNew,8,5) & " "); Wr.PutText(stderr, "change = " & FLR(TNew-TOld,8,5) & "\n"); END; (* Apply the reparametrization "tFilOld -> tFil" to "tOrg" too: *) IF cOrg.m < 2 OR TNew < minLength THEN IF verbose THEN Wr.PutText(stderr, "filtered curve collapsed, "); Wr.PutText(stderr, "length = " & FLR(TNew,8,5) & " "); Wr.PutText(stderr, "- original times not changed \n"); END; RETURN minLength ELSE (* Compute "arg[i]" = fractional index of "tOrgOld[i]" in "tFil": *) PZLRChain.ArgsFromValues(tFilOld, TOld, tOrg, arg); <* ASSERT arg[0] + 0.0d0 = arg[0] *> PZLRChain.ValuesFromArgs(tFil, TNew, arg, tOrgNew); END; cOrg.tPeriod := TNew; (* Compute variation in time intervals from "tOrg" to "tOrgNew": *) maxChange := 0.0d0; maxChangeIndex := 0; minStep := LAST(LONG); maxStep := FIRST(LONG); tao := tOrg[cOrg.m-1] - TOld; tan := tOrgNew[cOrg.m-1] - TNew; FOR i := 0 TO cOrg.m-1 DO WITH tbo = tOrg[i], tbn = tOrgNew[i], stepSize = tbn - tan, stepChange = (tbn-tan)-(tbo-tao) DO <* ASSERT tan <= tbn *> <* ASSERT tao <= tbo *> IF stepSize > maxStep THEN maxStep := stepSize END; IF stepSize < minStep THEN minStep := stepSize END; IF ABS(stepChange) > ABS(maxChange) THEN maxChange := stepChange; maxChangeIndex := i END; tao := tbo; tan := tbn END END; IF verbose THEN Wr.PutText(stderr, "new time steps: min = " & FLR(minStep,8,5) & " "); Wr.PutText(stderr, "max = " & FLR(maxStep,8,5) & "\n"); Wr.PutText(stderr, "max change = " & FLR(maxChange,8,5) & " "); Wr.PutText(stderr, "at sample " & Fmt.Int(maxChangeIndex) & "\n"); END; (* Now set the times in "cOrg" appropriately: *) cOrg.setTimes(tOrgNew, TNew); <* ASSERT cOrg.tPeriod + 0.0d0 = cOrg.tPeriod *> (* To catch NaN *) END; RETURN maxChange END END Reparametrize; PROCEDURE EnsureMinSize(VAR r: REF LONGS; n: NAT) = BEGIN IF r = NIL OR NUMBER(r^) < n THEN r := NEW(REF LONGS, (110 * n) DIV 100 + 10) END END EnsureMinSize; TYPE CurveData = RECORD cmt: TEXT; curve: PZCurve.T END; PROCEDURE ReadCurve(fileName: TEXT; verbose: BOOL): CurveData = VAR L: LONG; BEGIN Wr.PutText(stderr, "reading file " & fileName & " ...\n"); WITH chd = PZLR3Chain.Read(FileRd.Open(fileName), headerOnly := FALSE), p = chd.c^, m = NUMBER(p), c = PZCurve.New(m), s = NEW(REF LONGS, m)^ DO (* Translate the curve's barycenter to the origin: *) WITH b = PZLR3Chain.Barycenter(p) DO FOR i := 0 TO LAST(p) DO p[i] := LR3.Sub(p[i], b) END END; c.setSamples(p); (* Initial times are lengths along the polygonal: *) L := PZLR3Chain.LinearLengths(p, s); IF verbose THEN Wr.PutText(stderr, "polygonal length = " & FLR(L,8,5) & "\n"); END; c.setTimes(s, L); (* Recompute lengths along interpolated curve: *) c.uniformize(); IF verbose THEN Wr.PutText(stderr,"adjusted curve length = " & FLR(c.tPeriod,8,5) & "\n"); END; (* Use the natural parameter values of the original curve as labels: *) c.setLabels(c.t^, c.tPeriod); RETURN CurveData{cmt := chd.cmt, curve := c} END END ReadCurve; PROCEDURE WriteRawCurve( READONLY o: Options; cmt: TEXT; READONLY c: PZCurve.T; ) = BEGIN WriteCurve(o.output & "000", cmt, c); END WriteRawCurve; PROCEDURE WriteSampledCurve( READONLY o: Options; cmt: TEXT; READONLY c: PZCurve.T; lambda: LONG; totIter: NAT; ) = PROCEDURE FmtMaxStep(maxStep: LONG): TEXT = BEGIN IF maxStep # LAST(LONG) THEN RETURN " maxStep = " & FLR(maxStep,8,5) ELSE RETURN "" END END FmtMaxStep; BEGIN WITH name = o.output & Fmt.Pad(Fmt.Int(ROUND(lambda)), 3, '0'), L = c.tPeriod, step = MIN(o.maxStep, lambda/4.0d0), n = CEILING(L/step) DO WriteCurve( name, cmt & " filtered to scale lambda = " & FLR(lambda,8,5) & " (pixels)\n" & " nRepStages = " & Fmt.Int(o.nRepStages) & " tolLength = " & FLR(o.tolLength,8,5) & " tolStep = " & FLR(o.tolStep,8,5) & " maxIter = " & Fmt.Int(o.maxIter) & "\n" & " sampled at " & Fmt.Int(n) & " equally spaced points" & FmtMaxStep(o.maxStep) & "\n" & " totIter = " & Fmt.Int(totIter) & ")" & "\n", SampleCurve(c, n) ); END END WriteSampledCurve; PROCEDURE SampleCurve(c: PZCurve.T; n: NAT): PZCurve.T = (* Resamples a curve at "n" points equally spaced in time. *) BEGIN WITH tau = NEW(REF PZLRChain.T, n)^, p = NEW(REF PZLR3Chain.T, n)^, v = NEW(REF PZLR3Chain.T, n)^, r = NEW(REF PZLRChain.T, n)^, T = c.tPeriod, R = c.rPeriod, d = PZCurve.New(n) DO FOR i := 0 TO n-1 DO tau[i] := FLOAT(i,LONG)/FLOAT(n,LONG) * T END; c.sample(tau, p, v, r); d.setSamples(p); d.setTimes(tau, T); d.setLabels(r, R); RETURN d END END SampleCurve; PROCEDURE WriteCurve( name: TEXT; cmt: TEXT; READONLY c: PZCurve.T; ) = BEGIN WITH fileName = name & ".lbl" DO Wr.PutText(stderr, "writing file " & fileName & " ...\n"); PZLRChain.Write( FileWr.Open(fileName), cmt, c.r^, c.rPeriod, unit := 0.001d0 ) END; WITH fileName = name & ".flc" DO Wr.PutText(stderr, "writing file " & fileName & " ...\n"); PZLR3Chain.Write( FileWr.Open(fileName), cmt, c.p^, unit := 0.001d0 ) END; END WriteCurve; PROCEDURE GetOptions(): Options = VAR o: Options; BEGIN WITH pp = NEW(ParseParams.T).init(stderr) DO TRY pp.getKeyword("-input"); o.input := pp.getNext(); pp.getKeyword("-output"); o.output := pp.getNext(); pp.getKeyword("-lambdaMin"); o.lambdaMin := pp.getNextLongReal(); IF pp.keywordPresent("-nScales") THEN o.nScales := pp.getNextInt(0, LAST(NAT)) ELSE o.nScales := LAST(NAT) END; IF pp.keywordPresent("-parametric") THEN o.maxIter := 0; o.tolLength := 0.0d0; o.tolStep := 0.0d0 ELSE pp.getKeyword("-tolLength"); o.tolLength := pp.getNextLongReal(0.0d0); pp.getKeyword("-tolStep"); o.tolStep := pp.getNextLongReal(0.0d0); IF pp.keywordPresent("-maxIter") THEN o.maxIter := pp.getNextInt(0, LAST(NAT)) ELSE o.maxIter := LAST(NAT) END END; IF pp.keywordPresent("-nRepStages") THEN o.nRepStages := pp.getNextInt(1, LAST(NAT)) ELSE o.nRepStages := 1 END; IF pp.keywordPresent("-maxStep") THEN o.maxStep := pp.getNextLongReal(1.0d-6, LAST(LONG)) ELSE o.maxStep := LAST(LONG) END; IF pp.keywordPresent("-lambdaRatio") THEN o.lambdaRatio := pp.getNextLongReal(); ELSE o.lambdaRatio := 2.0d0; END; o.writeRaw := pp.keywordPresent("-writeRaw"); IF pp.keywordPresent("-writeStep") THEN o.writeStep := pp.getNextInt(1, LAST(NAT)) ELSE o.writeStep := 1 END; o.verbose := pp.keywordPresent("-verbose"); pp.finish(); EXCEPT | ParseParams.Error => Wr.PutText(stderr, "Usage: PZFilter \\\n"); Wr.PutText(stderr, " -input NAME \\\n"); Wr.PutText(stderr, " -output NAME \\\n"); Wr.PutText(stderr, " -lambdaMin NUM \\\n"); Wr.PutText(stderr, " [ -nScales NUM ] \\\n"); Wr.PutText(stderr, " [ -nRepStages NUM ] \\\n"); Wr.PutText(stderr, " [ -lambdaRatio NUM] \\\n"); Wr.PutText(stderr, " { -parametric | \\\n"); Wr.PutText(stderr, " -tolLength NUM -tolStep NUM [ -maxIter NUM ] \\\n"); Wr.PutText(stderr, " } \\\n"); Wr.PutText(stderr, " [ -maxStep NUM ] \\\n"); Wr.PutText(stderr, " [ -verbose ] \\\n"); Wr.PutText(stderr, " [ -writeRaw ] \\\n"); Wr.PutText(stderr, " [ -writeStep NUM ] \n"); Process.Exit(1); END; END; RETURN o END GetOptions; PROCEDURE FLR(x: LONG; width, prec: NAT): TEXT = BEGIN RETURN Fmt.Pad(Fmt.LongReal(x, prec := prec, style := Fmt.Style.Fix), width) END FLR; BEGIN Main() END PZFilter. (* Copyright © 2001 Universidade Estadual de Campinas (UNICAMP). Authors: Helena C. G. Leitão and Jorge Stolfi. This file can be freely distributed, used, and modified, provided that this copyright and authorship notice is preserved, and that any modified versions are clearly marked as such. This software has NO WARRANTY of correctness or applicability for any purpose. Neither the authors nor their employers chall be held responsible for any losses or damages that may result from its use. *)