Retrieving global forces/moments/cop using ForcePlatformWrenchFilter

[Ben van Basten]

Hello everyone, 

I'm trying to retrieve the global forces, moments and CoP from a C3D (attached, TYPE-2 force plate from Nexus). 

Basically, the code is as follows: applying a ForcePlatformExtractor, feeding its results to a ForcePlatformWrenchFilter and retrieving those wrenches.

As far as I see, the resulting wrenches should contain the global force, moment and CoP. However, I can't seem to get the right CoP (possibly due to wrong moments). Global forces appear to be fine. Code can be seen below.

First issue:

The wrench position is always the same, it seems to be the position of the force plate itself. Why is this?

I tried recomputing the CoP using the global moments from the wrenches, which, in the Vicon coordinate system would be:

CoPX = -MomentY / ForceZ;

CoPY = MomentX / ForceZ;

However, the CoPs are incorrect so I guess the moments are not in lab space (but in a sort-off force plate coordinate system). 

I then decided to try a GroundReactionWrenchFilter.

Second issue

I replaced my ForcePlatformWrenchFilter with a GroundReactionWrenchFilter.  Here, the wrench position is changing over time and the resulting wrench position is very close to the CoP as reported by Vicon (and Visual3D), but not exactly right. Could this be the difference between point-of-application and center-of-pressure? 

If this is the case, is there a way to use the GroundReactionWrenchFilter such that it reports the CoP instead?

Vicon:

cop:294.078, 308.603, 0

Visual3D

cop: 0.294078, 0.308603,0.000000

BTK/Mokka

cop: 293.646, 308.658, 0

Our software 

cop: -0.294078, 0, 0.308603

  

Third issue

The reported moments differ (not even close) to the global moments as reported by other software (or compared to the equation above). I think I can derive the the global moments myself (using the free moment), but it is unclear what the WrenchFilters are reporting. Are these the moments in force plate space? 

We really would like to use BTK in our projects, as our implementation is limited to Vicon C3D's, so any help would be greatly appreciated.

Regards,

Ben van Basten

// RETRIEVING THE WRENCH COLLECTION

if (!m_P->GetBtkC3DFileIO()->CanReadFile(filename))

       return false;

 

m_P->GetBtkC3DFileIO()->Read(filename, m_P->GetBtkAcquisition());

// Retrieve force platforms

auto forceExtractor = btk::ForcePlatformsExtractor::New();

forceExtractor->SetInput(m_P->GetBtkAcquisition());

forceExtractor->Update();

btk::ForcePlatformCollection::Pointer forcePlatforms = forceExtractor->GetOutput();

// Precompute global forces and moments

auto wrenchFilter = btk::ForcePlatformWrenchFilter::New();

wrenchFilter->SetTransformToGlobalFrame(true);

wrenchFilter->SetInput(forcePlatforms);

wrenchFilter->Update();

btk::WrenchCollection::Pointer pBtkwrenches = wrenchFilter->GetOutput();

When I need the global force and moment from a particular frame I query the wrench collection as follows

// QUERYING THE WRENCH COLLECTION

bool C3DFile::GetForceChannelValues(size_t plate, size_t frame, std::vector<std::vector<double> >& results) const

{

if ((int)plate >= m_P->m_pBtkwrenches->GetItemNumber())

return false;

      results.resize(GetNbAnalogSamplesPerFrame(), std::vector<double>(9, 0.0));

for ( size_t s = 0; s < GetNbAnalogSamplesPerFrame(); s++)

{

size_t offset = frame*GetNbAnalogSamplesPerFrame();

if ( offset + s >= size_t(m_P->m_pBtkwrenches->GetItem(plate)->GetForce()->GetValues().size()) )

return false;

if ( offset + s >= size_t(m_P->m_pBtkwrenches->GetItem(plate)->GetMoment()->GetValues().size()) )

return false;

auto force = m_P->m_pBtkwrenches->GetItem(plate)->GetForce()->GetValues().row(offset + s);

results[s][0] = force[0];

results[s][1] = force[1];

results[s][2] = force[2];

auto moment = m_P->m_pBtkwrenches->GetItem(plate)->GetMoment()->GetValues().row(offset + s);

results[s][3] = moment[0];

results[s][4] = moment[1];

results[s][5] = moment[2];

                auto position = m_P->m_pBtkwrenches->GetItem(plate)->GetPosition()->GetValues().row(offset + s);
                results[s][6] = position[0];
                results[s][7] = position[1];
                results[s][8] = position[2];

}

return true;

}

[Arnaud Barré]

Dear Ben,

The COP is making assumptions that the PWA does not do. The COP’s formula assumes that horizontal forces are null which might be true (or close to) for a gait task but not for other tasks (e.g. running, jump, etc.). I decided then to use the PWA.

Regarding your third issue, you need to use also the result computed by the class GroundReactionWrenchFilter. The moments computed by the class ForcePlatformWrenchFilter are expressed to the centre surface origin of the force plate. In case of the class GroundReactionWrenchFilter, the result is computed at the PWA.

Regards,

Arnaud

> --
> You received this message because you are subscribed to the Google Groups "BTK Users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to btk-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <Walk01_2.c3d>

[Ben van Basten]

Dear Arnaud,

Thanks for your swift reply. I'll see if I can subclass GroundReactionWrenchFilter such that it calculates the COP (I think it would be a nice enhancement if the user can select which wrench position (COP or PWA) is being calculated.)

Regards, Ben

[Arnaud BARRE]

Dear Ben,

I have just introduced a new library called OpenMA which might simplify your work. It is in heavy development, but you can already compute the force and moment at the COP and not the PWA (in fact you have the choice).

// C++11 code
// Load the content of the C3D file

auto root = ma::io::read("yourfile.c3d");

// Extract all the force plates found

auto fps = root.findChildren<ma::instrument::ForcePlate*>();

// For each force plates, compute the wrench at the centre of pressure expressed in the glbal frame.

std::vector<ma::TimeSequence*> wrenches;

for (const auto& fp : fps)

  wrenches.push_back(fp->wrench(ma::instrument::Location::CentreOfPressure));

// ...

// IMPORTANT: Delete the variable root as this was allocated on the heap.

delete root;

This example will replace all the BTK code you use (the result of a ma::TimeSequence for a wrench is an array of 10 by n, where n is the number of samples. The 10th column is for the validity of the sample.).

I hope this will help you in your work

Regards,

Arnaud

[Ben van Basten]

Hi Arnaud,

that looks great! Is the C3D support and force computation stable enough? I'll try it out and see if it is suitable for our products.

I do, however, also want to go the BTK route (for its stability)... Would it be difficult for me to adjust the GroundReactionForceFilter to calculate CoP instead of PWA?

Ben

Op maandag 29 augustus 2016 19:51:44 UTC+2 schreef Arnaud BARRE:

[Arnaud Barré]

Hi Ben,

The code is stable, already with unit tests and results were verified. All the force plates type are not yet implemented only due to a matter of time. What type of force plate do you use in your C3D file? We can compare results between our implementations if you are interested.

Regarding the computation of force and moment at the COP instead of the PWA in BTK, you can inherit from the class btk::GroundReactionWrenchFilter and overload the function FinishGRWComputation. You have to copy its content and modify lines 101-102 [1]. A better solution would be to modify the class btk::GroundReactionWrenchFilter and offer the choice to compute force and moment at the COP or PWA. However, this is out of the scope of the support I planned for BTK.

Arnaud

[1] https://github.com/Biomechanical-ToolKit/BTKCore/blob/master/Code/BasicFilters/btkGroundReactionWrenchFilter.h#L101

[Ben van Basten]

Hi Arnaud, 

I am currently working on type 2 (Vicon) and type 4 (MotionAnalysis and Qualisys). 

I'll modify GroundReactionWrenchFilter accordingly (add a COP/PWA option as well) to see if I can get it to work. We can then compare to resulting CoP in BTK and openMA.

Regards,

Ben

Op dinsdag 30 augustus 2016 15:16:58 UTC+2 schreef Arnaud BARRE: