Commit 731f5a89 authored by Xavier de Blas's avatar Xavier de Blas

Added Joan Charmant Kinovea code for Butterworth filtering

parent 7b26123a
d=read.csv2("sample-signal-output.csv")
xini = 6545
xend = 7263
#trajectory raw and filtered
plot(d$Raw.coordinate..mm., type="l", xlim=c(xini,xend), axes=F, ylab="")
lines(d$Filtered.coordinate..mm., type="l", xlim=c(xini,xend), col="yellow")
#speed
speed = d$Velocity..mm.ms.
speedMax = max(abs(speed[xini:xend]), na.rm=T)
par(new=T)
plot(speed, type="l", xlim=c(xini,xend), col="green", ylim=c(-speedMax * 1.1, speedMax * 1.1), axes=F, ylab="speed")
axis(2)
#accel
accel = d$Acceleration..mm.ms.. * 1000
accelMax = max(abs(accel[xini:xend]), na.rm=T)
par(new=T)
plot(accel, type="l", xlim=c(xini,xend), col="red", ylim=c(-accelMax * 1.1, accelMax * 1.1), axes=F, ylab="")
mtext(side=4, "accel")
axis(4)
#end
box()
This diff is collapsed.
This diff is collapsed.
#region License
/*
Copyright © 2017 Joan Charmant
The MIT license.
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/
#endregion
using System;
using System.Collections.Generic;
namespace Kinovea.Filtering
{
/// <summary>
/// Butterworth low-pass filter. Fourth order, zero-lag low pass filter. Two 2nd order passes (forward/backward) are used to reset the phase shift.
/// Ref: "Biomechanics and Motor control of human movement", David A. Winter, 4th ed. 2009.
///
/// To initialize the filter the trajectory is extrapolated for 10 data points each side using reflected values around the end points.
/// The extrapolated points are then removed from the filtered results.
/// Ref: "Padding point extrapolation techniques for the butterworth digital filter", Gerald Smith, J. Biomechonics Vol. 22, No. s/9, pp. 967-971, 1989.
///
/// The filter is passed at various cutoff frequencies between 0.5Hz and the Nyquist frequency, and all results are kept.
///
/// The best-guess cutoff frequency is found by estimating autocorrelation of residuals and taking the frequency yielding the least autocorrelated residuals.
/// Ref: "A procedure for the automatic determination of filter cutoff frequency for the processing of biomechanical data", John Challis, JAB, 1999, 15, 303-317.
///
/// The autocorrelation of residuals is estimated using the Durbin-Watson statistic.
/// </summary>
public class ButterworthFilter
{
private double a0, a1, a2;
private double b1, b2;
private double x0, x1, x2, y0, y1, y2;
private double correctionFactor;
/// <summary>
/// Filter a list of samples and return a list of lists of filtered values at various test cutoff frequencies.
/// The number of samples must be >= 10.
/// </summary>
/// <param name="samples">The raw values to be filtered</param>
/// <param name="fs">The sampling frequency</param>
/// <param name="fcTests">The number of cutoff frequencies to test. These will span between 0.5Hz and the Nyquist frequency.</param>
/// <param name="bestCutoff">The best guess cutoff frequency found at minimal autocorrelation of residuals.</param>
public List<FilteringResult> FilterSamples(double[] samples, double fs, int fcTests, out int bestCutoffIndex)
{
if (samples.Length <= 10)
throw new ArgumentException("Number of samples must be superior to 10");
List<FilteringResult> results = new List<FilteringResult>();
double nyquist = fs / 2;
UpdateCorrectionFactor(2);
int padding = 10;
double[] padded = AddPadding(samples, padding);
// Compute filtered result for a range of fc.
// We keep the whole array of results so the user can manually change cutoff frequency without a
// complete recomputation, and it is also necessary to compute the best-guess value.
double min = 0.5;
double max = nyquist;
double step = (max - min) / fcTests;
double bestScore = 1;
int bestIndex = -1;
int index = 0;
for (double fc = min; fc < max; fc += step)
{
double[] filteredPadded = FilterSamples(padded, fs, fc);
double[] filtered = RemovePadding(filteredPadded, padding);
double[] residuals = samples.Subtract(filtered);
double dw = StatsHelper.DurbinWatson(residuals);
if (double.IsNaN(dw))
continue;
double dwNormalized = Math.Abs(2 - dw) / 2;
FilteringResult result = new FilteringResult(fc, filtered, dwNormalized);
results.Add(result);
if (dwNormalized < bestScore)
{
bestScore = dwNormalized;
bestIndex = index;
}
index++;
}
bestCutoffIndex = bestIndex;
return results;
}
private double[] AddPadding(double[] samples, int padding)
{
// Extrapolation of trajectory using reflection of values around the end points.
// Ref: "Padding point extrapolation techniques for the butterworth digital filter". Smith 1989.
double[] padded = new double[samples.Length + 2 * padding];
for (int i = 0; i < padding; i++)
padded[i] = samples[0] + samples[0] - samples[padding - i];
for (int i = 0; i < samples.Length; i++)
padded[padding + i] = samples[i];
for (int i = 0; i < padding; i++)
padded[padding + samples.Length + i] = samples[samples.Length - 1] + samples[samples.Length - 1] - samples[samples.Length - 1 - i - 1];
return padded;
}
private double[] RemovePadding(double[] samples, int padding)
{
double[] result = new double[samples.Length - 2 * padding];
for (int i = padding; i < samples.Length - padding; i++)
result[i - padding] = samples[i];
return result;
}
private double[] FilterSamples(double[] samples, double fs, double fc)
{
UpdateCoefficients(fs, fc);
double[] forward = ForwardPass(samples);
double[] backward = BackwardPass(forward);
return backward;
}
private void UpdateCorrectionFactor(int passes)
{
// Ref: Chapt. 3.4.4.2 of "Biomechanics and motor control of human movement".
correctionFactor = Math.Pow((Math.Pow(2, 1.0 / passes) - 1), 0.25);
}
private void UpdateCoefficients(double fs, double fc)
{
// Ref: Chapt. 2.2.4.4 of "Biomechanics and motor control of human movement".
double o = Math.Tan(Math.PI * fc / fs) / correctionFactor;
double k1 = MathHelper.SQRT2 * o;
double k2 = o * o;
a0 = k2 / (1 + k1 + k2);
a1 = 2 * a0;
a2 = a0;
double k3 = 2 * a0 / k2;
b1 = -2 * a0 + k3;
b2 = 1 - a0 - a1 - a2 - b1;
}
private void ResetValues()
{
x0 = x1 = x2 = y0 = y1 = y2 = 0;
}
private double[] ForwardPass(double[] raw)
{
ResetValues();
y1 = x1 = raw[0];
y0 = x0 = raw[1];
double[] filtered = new double[raw.Length];
filtered[0] = y1;
filtered[1] = y0;
for (int i = 2; i < raw.Length; i++)
{
double f = FilterSample(raw[i]);
filtered[i] = f;
}
return filtered;
}
private double[] BackwardPass(double[] forward)
{
ResetValues();
y1 = x1 = forward[forward.Length - 1];
y0 = x0 = forward[forward.Length - 2];
double[] filtered = new double[forward.Length];
filtered[0] = y1;
filtered[1] = y0;
for (int i = forward.Length - 3; i >= 0; i--)
{
double f = FilterSample(forward[i]);
filtered[forward.Length - 1 - i] = f;
}
Array.Reverse(filtered);
return filtered;
}
private double FilterSample(double sample)
{
x2 = x1;
x1 = x0;
x0 = sample;
y2 = y1;
y1 = y0;
y0 = a0 * x0 + a1 * x1 + a2 * x2 + b1 * y1 + b2 * y2;
return y0;
}
}
}
#region License
/*
Copyright © 2017 Joan Charmant
The MIT license.
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/
#endregion
using System;
using System.Collections.Generic;
namespace Kinovea.Filtering
{
public static class Extensions
{
public static List<double> Subtract(this List<double> a, List<double> b)
{
if (a.Count != b.Count)
throw new ArgumentException("Lists must have the same number of elements");
List<double> result = new List<double>();
for (int i = 0; i < a.Count; i++)
result.Add(a[i] - b[i]);
return result;
}
public static double[] Subtract(this double[] a, double[] b)
{
if (a.Length != b.Length)
throw new ArgumentException("Lists must have the same number of elements");
double[] result = new double[a.Length];
for (int i = 0; i < a.Length; i++)
result[i] = a[i] - b[i];
return result;
}
}
}
\ No newline at end of file
#region License
/*
Copyright © 2017 Joan Charmant
The MIT license.
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/
#endregion
using System;
using System.Collections.Generic;
using System.Drawing;
namespace Kinovea.Filtering
{
/// <summary>
/// Contains filtered and unfiltered data for one time series of 2D points and the corresponding cutoff frequencies.
/// </summary>
public class FilteredTrajectory
{
/// <summary>
/// Number of samples.
/// </summary>
public int Length { get; private set; }
/// <summary>
/// Whether the trajectory can be filtered, depends on the number of samples.
/// </summary>
public bool CanFilter { get; private set; }
/// <summary>
/// Time coordinates.
/// </summary>
public long[] Times { get; private set; }
/// <summary>
/// Raw coordinates.
/// </summary>
public double[] RawXs { get; private set; }
/// <summary>
/// Raw coordinates.
/// </summary>
public double[] RawYs { get; private set; }
/// <summary>
/// Filtered coordinates at the currently selected cutoff frequency.
/// </summary>
public double[] Xs
{
get
{
return XCutoffIndex < 0 ? RawXs : FilterResultXs[XCutoffIndex].Data;
}
}
/// <summary>
/// Filtered coordinates at the currently selected cutoff frequency.
/// </summary>
public double[] Ys
{
get
{
return YCutoffIndex < 0 ? RawYs : FilterResultYs[YCutoffIndex].Data;
}
}
/// <summary>
/// Filtered X coordinates time series at various cutoff frequencies.
/// </summary>
public List<FilteringResult> FilterResultXs { get; set; }
/// <summary>
/// Best-guess cutoff frequency index for Xs series.
/// </summary>
public int XCutoffIndex { get; set; }
/// <summary>
/// Filtered time series at various cutoff frequencies.
/// </summary>
public List<FilteringResult> FilterResultYs { get; set; }
/// <summary>
/// Best-guess cutoff frequency index for Ys series.
/// </summary>
public int YCutoffIndex { get; set; }
//public Circle BestFitCircle { get; private set; }
/// <summary>
/// Initialize the data and filter it if possible.
/// </summary>
public void Initialize(List<TimedPoint> samples, double captureFramesPerSecond)
{
this.Length = samples.Count;
Times = new long[samples.Count];
RawXs = new double[samples.Count];
RawYs = new double[samples.Count];
XCutoffIndex = -1;
YCutoffIndex = -1;
// Raw coordinates.
for (int i = 0; i < samples.Count; i++)
{
//PointF point = calibrationHelper.GetPointAtTime(samples[i].Point, samples[i].T);
PointF point = samples[i].Point;
RawXs[i] = point.X;
RawYs[i] = point.Y;
Times[i] = samples[i].T;
}
this.CanFilter = samples.Count > 10;
if (this.CanFilter)
{
//double framerate = calibrationHelper.CaptureFramesPerSecond;
double framerate = captureFramesPerSecond;
ButterworthFilter filter = new ButterworthFilter();
// Filter the results a hundred times at various cutoff frequency and store all data along with the best cutoff frequency.
int tests = 100;
int bestCutoffIndexX;
FilterResultXs = filter.FilterSamples(RawXs, framerate, tests, out bestCutoffIndexX);
XCutoffIndex = bestCutoffIndexX;
int bestCutoffIndexY;
FilterResultYs = filter.FilterSamples(RawYs, framerate, tests, out bestCutoffIndexY);
YCutoffIndex = bestCutoffIndexY;
}
//BestFitCircle = CircleFitter.Fit(this);
}
public PointF RawCoordinates(int index)
{
return new PointF((float)RawXs[index], (float)RawYs[index]);
}
public PointF Coordinates(int index)
{
float x = XCutoffIndex < 0 ? (float)RawXs[index] : (float)FilterResultXs[XCutoffIndex].Data[index];
float y = YCutoffIndex < 0 ? (float)RawYs[index] : (float)FilterResultYs[YCutoffIndex].Data[index];
return new PointF(x, y);
}
}
}
#region License
/*
Copyright © 2017 Joan Charmant
The MIT license.
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/
#endregion
using System;
namespace Kinovea.Filtering
{
public class FilteringResult
{
public double CutoffFrequency { get; private set; }
public double[] Data { get; private set; }
public double DurbinWatson { get; private set; }
public FilteringResult(double fc, double[] data, double dw)
{
this.CutoffFrequency = fc;
this.Data = data;
this.DurbinWatson = dw;
}
}
}
#region License
/*
Copyright © 2017 Joan Charmant
The MIT license.
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/
#endregion
using System;
namespace Kinovea.Filtering
{
/// <summary>
/// A list of possible kinematics values. These are used as keys for the time series collections.
/// </summary>
public enum Kinematics
{
X,
Y,
XRaw,
YRaw,
LinearDistance, // Accumulated distance since the starting point.
LinearHorizontalDisplacement, // Straight line displacement from the current point to the starting point, horizontal component.
LinearVerticalDisplacement, // Straight line displacement from the current point to the starting point, vertical component.
LinearSpeed,
LinearHorizontalVelocity,
LinearVerticalVelocity,
LinearAcceleration, // Instantaneous acceleration in the direction of the velocity vector.
LinearHorizontalAcceleration,
LinearVerticalAcceleration,
AngularPosition, // Absolute or relative value of the angle at this time.
AngularDisplacement, // Displacement with regards to previous point of measure.
TotalAngularDisplacement, // Total displacement with regards to first point of measure.
AngularVelocity,
TangentialVelocity,
AngularAcceleration,
TangentialAcceleration,
CentripetalAcceleration,
ResultantLinearAcceleration
}
}
<?xml version="1.0" encoding="utf-8"?>
<Project ToolsVersion="14.0" DefaultTargets="Build" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
<Import Project="$(MSBuildExtensionsPath)\$(MSBuildToolsVersion)\Microsoft.Common.props" Condition="Exists('$(MSBuildExtensionsPath)\$(MSBuildToolsVersion)\Microsoft.Common.props')" />
<PropertyGroup>
<Configuration Condition=" '$(Configuration)' == '' ">Debug</Configuration>
<Platform Condition=" '$(Platform)' == '' ">AnyCPU</Platform>
<ProjectGuid>{BF61942A-D05A-4299-BA24-C2861C8A0E36}</ProjectGuid>
<OutputType>Library</OutputType>
<AppDesignerFolder>Properties</AppDesignerFolder>
<RootNamespace>Filterer</RootNamespace>
<AssemblyName>Filterer</AssemblyName>
<TargetFrameworkVersion>v4.5</TargetFrameworkVersion>
<FileAlignment>512</FileAlignment>
</PropertyGroup>
<PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Debug|AnyCPU' ">
<DebugSymbols>true</DebugSymbols>
<DebugType>full</DebugType>
<Optimize>false</Optimize>
<OutputPath>bin\Debug\</OutputPath>
<DefineConstants>DEBUG;TRACE</DefineConstants>
<ErrorReport>prompt</ErrorReport>
<WarningLevel>4</WarningLevel>
</PropertyGroup>
<PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release|AnyCPU' ">
<DebugType>pdbonly</DebugType>
<Optimize>true</Optimize>
<OutputPath>bin\Release\</OutputPath>
<DefineConstants>TRACE</DefineConstants>
<ErrorReport>prompt</ErrorReport>
<WarningLevel>4</WarningLevel>
</PropertyGroup>
<ItemGroup>
<Reference Include="System" />
<Reference Include="System.Core" />
<Reference Include="System.Drawing" />
<Reference Include="System.Xml.Linq" />
<Reference Include="System.Data.DataSetExtensions" />
<Reference Include="Microsoft.CSharp" />
<Reference Include="System.Data" />
<Reference Include="System.Net.Http" />
<Reference Include="System.Xml" />
</ItemGroup>
<ItemGroup>
<Compile Include="Butterworth.cs" />
<Compile Include="Extensions.cs" />
<Compile Include="FilteredTrajectory.cs" />
<Compile Include="FilteringResult.cs" />
<Compile Include="Kinematics.cs" />
<Compile Include="LinearKinematics.cs" />
<Compile Include="MathHelper.cs" />
<Compile Include="MovingAverage.cs" />
<Compile Include="Properties\AssemblyInfo.cs" />
<Compile Include="StatsHelper.cs" />
<Compile Include="TimedPoint.cs" />
<Compile Include="TimeSeriesCollection.cs" />
<Compile Include="TimeSeriesPadder.cs" />
</ItemGroup>
<Import Project="$(MSBuildToolsPath)\Microsoft.CSharp.targets" />
<!-- To modify your build process, add your task inside one of the targets below and uncomment it.
Other similar extension points exist, see Microsoft.Common.targets.
<Target Name="BeforeBuild">
</Target>
<Target Name="AfterBuild">
</Target>
-->
</Project>
\ No newline at end of file
#region License
/*
Copyright © 2017 Joan Charmant
The MIT license.
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/
#endregion
using System;
using System.Collections.Generic;
using System.Drawing;
namespace Kinovea.Filtering
{
/// <summary>
/// Build time series data for various linear kinematics values.
/// Note: This class was simplified compared to Kinovea code, to focus on 1D and only compute speed and acceleration.
/// </summary>
public class LinearKinematics
{
public TimeSeriesCollection BuildKinematics(FilteredTrajectory traj, float intervalMilliseconds, bool enableSecondPassFiltering)
{
TimeSeriesCollection tsc = new TimeSeriesCollection(traj.Length);
if (traj.Length == 0)
return tsc;
tsc.AddTimes