Scale Ground Contact Model using Matlab

[Evan Dooley]

Hello,

I couldn't find the post, but there was some discussion about how to scale a ground contact model using the Matlab interface, as opposed to changing sphere locations in the OSIM file. I have been using this method that implements the contact model discussed in the Miller2021 paper (as closely as I could from the paper description: https://peerj.com/articles/11960/), and scales it to the size of each subject.

There are hard-coded lines that other users would need to update w.r.t. where the center of mass for the foot segments in the specific base model they are using. Also you need to have your model loaded into matlab, and set to a variable called "model" (but perhaps that was self explanatory).

Hope this can be useful to some of you. Happy to answer any questions.

Code: Select all

% contact parameters =
stiffness = 3067776; % stiffness from Umberger's paper
dissipation = 2; % miller model uses this stiffness & dissipation
    % Get a reference to the model's ground body
    ground = model.getGround();
    
    %% Pull existing foot bodies
    % Get right foot
    R_calc = model.getBodySet().get('calcn_r');
    R_toes = model.getBodySet().get('toes_r');
    
    % Get left foot
    L_calc = model.getBodySet().get('calcn_l');
    L_toes = model.getBodySet().get('toes_l');
    
    %% Get Foot Scale Factors
    CoM = osimVec3ToArray(L_toes.getMassCenter()); % either toes would work, calcn misses Z factor because CoM aligns with calc origin
    scaleX = CoM(1)/0.0346; % hard coding in location from base model = HumanPiG_v11.osim; base model is same for all subjects
    scaleY = CoM(2)/0.006;
    scaleZ = CoM(3)/0.0175;
    
    %% ---------------------------------------
    %% Construct Contact GEOMETRY
    %% ------------------------------------------
    %% FLOOR Contact
    % Set top of floor as a "contact half space", where a half space is
    % everything to one side of an infinite plane. In its local coordinate
    % system, all points for which x>0 are considered to be inside the
    % geometry. Its location and orientation properties can be used to move and
    % rotate it to represent other half spaces
    groundContactLocation = Vec3(0, 0, 0); %raise it to the top of the boxes representing the floor, sitting at global y = 0
    groundContactOrientation = Vec3(0,0,-1.57); % want the x>0 space to be everything in y>0 in the global frame, so rotate the surface about Z axis
    groundContactSpace = ContactHalfSpace(groundContactLocation,...
                                           groundContactOrientation,...
                                           ground); %should be an infinite surface, so I think I only have to attach it to one body
    groundContactSpace.setName('FloorContact'); 
    model.addContactGeometry(groundContactSpace);
    
    %% Spheres
    r = 0.01 * max([scaleX, scaleY, scaleZ]); % adjust radius based on the largest of the scale factors
    scaleCMx = 0.01 * scaleX;
    scaleCMy = 0.01 * scaleY;
    scaleCMz = 0.01 * scaleZ; % adjust the measurement = 1 cm to the scaled values

    % RIGHT
    % Posterior Heel
    thisSphere = ContactSphere();
    thisSphere.setRadius(r); % set the radius of the sphere, in meters
    thisSphere.setFrame(R_calc); % this is the body object to attach to
    thisSphere.setLocation(Vec3(scaleCMx*1.5,0,0)); % this is the location of the center of the sphere relative to the coord frame of that body
    thisSphere.setName('R_PostHeel_Contact'); % the name that'll appear in the 'ContactGeometrySet' list
    model.addContactGeometry(thisSphere); % add it to the model object
    sphereList.R_PostHeel_Contact = thisSphere; % add the "ContactSphere" object to a list, to ease creating the contact force objects
    % Medial Heel
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(R_calc);
    thisSphere.setLocation(Vec3(scaleCMx*3,-scaleCMy,-scaleCMz*2.5));
    thisSphere.setName('R_MedHeel_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.R_MedHeel_Contact = thisSphere;
    % Lateral Heel
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(R_calc);
    thisSphere.setLocation(Vec3(scaleCMx*3,-scaleCMy,scaleCMz*2.5));
    thisSphere.setName('R_LatHeel_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.R_LatHeel_Contact = thisSphere;
    % Proximal Met1
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(R_toes);
    thisSphere.setLocation(Vec3(-scaleCMx*7.5,0,-scaleCMz*2.5));
    thisSphere.setName('R_ProxMet1_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.R_ProxMet1_Contact = thisSphere;
    % Proximal Met3
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(R_toes);
    thisSphere.setLocation(Vec3(-scaleCMx*8.5,-scaleCMy/2,0));
    thisSphere.setName('R_ProxMet3_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.R_ProxMet3_Contact = thisSphere;
    % Proximal Met5
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(R_toes);
    thisSphere.setLocation(Vec3(-scaleCMx*10,-scaleCMy,scaleCMz*2.5));
    thisSphere.setName('R_ProxMet5_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.R_ProxMet5_Contact = thisSphere;
    % Distal Met1
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(R_toes);
    thisSphere.setLocation(Vec3(0,-scaleCMy,-scaleCMz*2.5));
    thisSphere.setName('R_DistMet1_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.R_DistMet1_Contact = thisSphere;
    % Distal Met3
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(R_toes);
    thisSphere.setLocation(Vec3(-scaleCMx*2,-scaleCMy,scaleCMz));
    thisSphere.setName('R_DistMet3_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.R_DistMet3_Contact = thisSphere;
    % Distal Met5
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(R_toes);
    thisSphere.setLocation(Vec3(-scaleCMx*4.5,-scaleCMy,scaleCMz*3.5));
    thisSphere.setName('R_DistMet5_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.R_DistMet5_Contact = thisSphere;
    % Hallux
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(R_toes);
    thisSphere.setLocation(Vec3(scaleCMx*2,-scaleCMy/2,-scaleCMz*2));
    thisSphere.setName('R_Hallux_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.R_Hallux_Contact = thisSphere;
    % Other Toes
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(R_toes);
    thisSphere.setLocation(Vec3(scaleCMx*1,0,scaleCMz*2));
    thisSphere.setName('R_Toes_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.R_Toes_Contact = thisSphere;

    % LEFT
    % Posterior Heel
    thisSphere = ContactSphere();
    thisSphere.setRadius(r); % set the radius of the sphere, in meters
    thisSphere.setFrame(L_calc); % this is the body object to attach to
    thisSphere.setLocation(Vec3(scaleCMx*1.5,0,0)); % this is the location of the center of the sphere relative to the coord frame of that body
    thisSphere.setName('L_PostHeel_Contact'); % the name that'll appear in the 'ContactGeometrySet' list
    model.addContactGeometry(thisSphere); % add it to the model object
    sphereList.L_PostHeel_Contact = thisSphere; % add the "ContactSphere" object to a list, to ease creating the contact force objects
    % Medial Heel
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(L_calc);
    thisSphere.setLocation(Vec3(scaleCMx*3,-scaleCMy,scaleCMz*2.5));
    thisSphere.setName('L_MedHeel_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.L_MedHeel_Contact = thisSphere;
    % Lateral Heel
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(L_calc);
    thisSphere.setLocation(Vec3(scaleCMx*3,-scaleCMy,-scaleCMz*2.5));
    thisSphere.setName('L_LatHeel_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.L_LatHeel_Contact = thisSphere;
    % Proximal Met1
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(L_toes);
    thisSphere.setLocation(Vec3(-scaleCMx*7.5,0,scaleCMz*2.5));
    thisSphere.setName('L_ProxMet1_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.L_ProxMet1_Contact = thisSphere;
    % Proximal Met3
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(L_toes);
    thisSphere.setLocation(Vec3(-scaleCMx*8.5,-scaleCMy/2,0));
    thisSphere.setName('L_ProxMet3_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.L_ProxMet3_Contact = thisSphere;
    % Proximal Met5
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(L_toes);
    thisSphere.setLocation(Vec3(-scaleCMx*10,-scaleCMy,-scaleCMz*2.5));
    thisSphere.setName('L_ProxMet5_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.L_ProxMet5_Contact = thisSphere;
    % Distal Met1
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(L_toes);
    thisSphere.setLocation(Vec3(0,-scaleCMy,scaleCMz*2.5));
    thisSphere.setName('L_DistMet1_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.L_DistMet1_Contact = thisSphere;
    % Distal Met3
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(L_toes);
    thisSphere.setLocation(Vec3(-scaleCMx*2,-scaleCMy,-scaleCMz));
    thisSphere.setName('L_DistMet3_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.L_DistMet3_Contact = thisSphere;
    % Distal Met5
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(L_toes);
    thisSphere.setLocation(Vec3(-scaleCMx*4.5,-scaleCMy,-scaleCMz*3.5));
    thisSphere.setName('L_DistMet5_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.L_DistMet5_Contact = thisSphere;
    % Hallux
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(L_toes);
    thisSphere.setLocation(Vec3(scaleCMx*2,-scaleCMy/2,scaleCMz*2));
    thisSphere.setName('L_Hallux_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.L_Hallux_Contact = thisSphere;
    % Other Toes
    thisSphere = ContactSphere();
    thisSphere.setRadius(r);
    thisSphere.setFrame(L_toes);
    thisSphere.setLocation(Vec3(scaleCMx*1,0,-scaleCMz*2));
    thisSphere.setName('L_Toes_Contact');
    model.addContactGeometry(thisSphere);
    sphereList.L_Toes_Contact = thisSphere;

    %% Add Contact Force --> SmoothSphereHalfSpaceForce
    % Similar to the Hunt Crossley Model, but makes smooth transition from zero
    % to some force, and enforces non-null derivatives. This is necessary for
    % solving with the Moco optimization solver
    % --- Varying parameters ----------
    staticFriction      = 0.8;
    dynamicFriction     = 0.8;
    viscousFriction     = 0.5;
    transitionVelocity  = 0.2;
    % smoothing parameters --> defaults
    hertzSmoothing = 300; %The parameter that determines the smoothness of the transition of the tanh used to smooth the Hertz force. The larger the steeper the transition but also the more discontinuous-like, default is 300.
    huntCrossleySmoothing = 50; % The parameter that determines the smoothness of the transition of the tanh used to smooth the Hunt-Crossley force. The larger the steeper the transition but also the more discontinuous-like, default is 50.
    
    % call = 
    % matlabVariabletypeSSHSF = SmoothSphereHalfSpaceForce('name',contact sphere object, contact half space object);
    
    % pull all of the contact geometries
    allContactGeoObj = model.getContactGeometrySet();
    % make force for each
    for s = 1:allContactGeoObj.getSize()-1 % opensim starts counting at 0, but the force at the first index (i = 0) is the ground half space
        % get contact object
        currentSphereGeo = allContactGeoObj.get(s);
        sphereName = char(currentSphereGeo.getName()); % object type = "char"
        currentSphere = sphereList.(sphereName); % object type = "ContactSphere"
        forceName = strrep(sphereName,'Contact','Force');
        % create contact force
        thisContactForce = SmoothSphereHalfSpaceForce(forceName,currentSphere,groundContactSpace);% set all of the parameters for the force object
        thisContactForce.set_stiffness(stiffness);
        thisContactForce.set_dissipation(dissipation);
        thisContactForce.set_static_friction(staticFriction);
        thisContactForce.set_dynamic_friction(dynamicFriction);
        thisContactForce.set_viscous_friction(viscousFriction);
        thisContactForce.set_transition_velocity(transitionVelocity);
        thisContactForce.set_hertz_smoothing(hertzSmoothing);
        thisContactForce.set_hunt_crossley_smoothing(huntCrossleySmoothing);
        % add to model
        model.addForce(thisContactForce);
        % clean up... just in case
        clear thisContactForce currentSphereGeo currentSphere sphereName forceName
    end % next contact sphere

[Ross Miller]

Looks nice, thanks Evan.

The supplementary material/data for that paper should have the actual .osim file included but I can send it or post it here if it's not there, if you want to see some other detail of the contact model. They look like the code below.

The stiffness and damping we used were just whatever the values in the original 2D_gait.osim model were. For walking and running at least, I don't find those parameters have much effect, unless I do something silly like damping = A BILLION.

Ton once suggested changing all the friction coefficients to the same value, to avoid high-frequency oscillations in the shear forces. I still have different values here, 0.8 for static and dynamic and 0.5 for viscous. It's a fairly common problem in predictive simulations although I don't always see it. Have never seen it in a tracking simulation in Moco.

Ross

Code: Select all

<SmoothSphereHalfSpaceForce name="contactHeel_r">
	<!--Path to a Component that satisfies the Socket 'sphere' of type ContactSphere (description: The sphere participating in this contact.).-->
	<socket_sphere>/contactgeometryset/heel_r</socket_sphere>
	<!--Path to a Component that satisfies the Socket 'half_space' of type ContactHalfSpace (description: The half-space participating in this contact.).-->
	<socket_half_space>/contactgeometryset/floor</socket_half_space>
	<!--The stiffness constant (i.e., plain strain modulus), default is 1 (N/m^2)-->
	<stiffness>3067776</stiffness>
	<!--The dissipation coefficient, default is 0 (s/m).-->
	<dissipation>2</dissipation>
	<!--The coefficient of static friction, default is 0.-->
	<static_friction>0.80000000000000004</static_friction>
	<!--The coefficient of dynamic friction, default is 0.-->
	<dynamic_friction>0.80000000000000004</dynamic_friction>
	<!--The coefficient of viscous friction, default is 0.-->
	<viscous_friction>0.5</viscous_friction>
	<!--The transition velocity, default is 0.01 (m/s).-->
	<transition_velocity>0.20000000000000001</transition_velocity>
	<!--The parameter that determines the smoothness of the transition of the tanh used to smooth the Hertz force. The larger the steeper the transition but the worse for optimization, default is 300.-->
	<hertz_smoothing>300</hertz_smoothing>
	<!--The parameter that determines the smoothness of the transition of the tanh used to smooth the Hunt-Crossley force. The larger the steeper the transition but the worse for optimization, default is 50.-->
	<hunt_crossley_smoothing>50</hunt_crossley_smoothing>
</SmoothSphereHalfSpaceForce>

Code: Select all

<ContactSphere name="heel_r">
	<!--Path to a Component that satisfies the Socket 'frame' of type PhysicalFrame (description: The frame to which this geometry is attached.).-->
	<socket_frame>/bodyset/calcn_r</socket_frame>
	<!--Location of geometry center in the PhysicalFrame.-->
	<location>0.0013075275819317961 -0.0080999999999999996 0</location>
	<!--Radius of the sphere (default: 0).-->
	<radius>0.014999999999999999</radius>
</ContactSphere>

[Nicholas Bianco]

Thanks Evan!

[Evan Dooley]

Glad to hear others have had some discussion about how the friction coefficients relate to the high-frequency oscillations. I do see them in my simulations using MocoTrack and have started trying to reduce them. My hypothesis is that I get them because I have no goal related to matching recorded GRX forces. I want to be able to see how GRX changes, so I wanted to be able to validate that my simulation method could predict GRX forces to match recorded forces with no knowledge of the recorded forces.

The model that I found in the supplementary materials for the Miller2021 paper only has 5 contact objects on each foot (https://simtk.org/projects/umocod/):

Code: Select all

<ContactGeometrySet name="contactgeometryset">
			<objects>
                <ContactSphere name="heel_r">
                    <socket_frame>/bodyset/calcn_r</socket_frame>
                    <radius>0.035000000000000003</radius>
                    <location>0.031307527581931796 0.010435842527310599 0</location>
                </ContactSphere>
                <ContactSphere name="mh1_r">
                    <socket_frame>/bodyset/calcn_r</socket_frame>
                    <radius>0.014999999999999999</radius>
                    <location>0.1820 -0.0081 -0.0230</location>
                </ContactSphere>
                <ContactSphere name="mh3_r">
                    <socket_frame>/bodyset/calcn_r</socket_frame>
                    <radius>0.014999999999999999</radius>
                    <location>0.1658 -0.0081 0.0116</location>
                </ContactSphere>
                <ContactSphere name="mh5_r">
                    <socket_frame>/bodyset/calcn_r</socket_frame>
                    <radius>0.014999999999999999</radius>
                    <location>0.1430 -0.0081 0.0370</location>
                </ContactSphere>
                <ContactSphere name="hallux_r">
                    <socket_frame>/bodyset/toes_r</socket_frame>
                    <radius>0.014999999999999999</radius>
                    <location>0.0562 -0.0102 -0.0181</location>
                </ContactSphere>
                <ContactSphere name="othertoes_r">
                    <socket_frame>/bodyset/toes_r</socket_frame>
                    <radius>0.014999999999999999</radius>
                    <location>0.0441 -0.006 0.0242</location>
                </ContactSphere>
                
                <ContactSphere name="heel_l">
                    <socket_frame>/bodyset/calcn_l</socket_frame>
                    <radius>0.035000000000000003</radius>
                    <location>0.031307527581931796 0.010435842527310599 0</location>
                </ContactSphere>
                <ContactSphere name="mh1_l">
                    <socket_frame>/bodyset/calcn_l</socket_frame>
                    <radius>0.014999999999999999</radius>
                    <location>0.1820 -0.0081 0.0230</location>
                </ContactSphere>
                <ContactSphere name="mh3_l">
                    <socket_frame>/bodyset/calcn_l</socket_frame>
                    <radius>0.014999999999999999</radius>
                    <location>0.1658 -0.0081 -0.0116</location>
                </ContactSphere>
                <ContactSphere name="mh5_l">
                    <socket_frame>/bodyset/calcn_l</socket_frame>
                    <radius>0.014999999999999999</radius>
                    <location>0.1430 -0.0081 -0.0370</location>
                </ContactSphere>
                <ContactSphere name="hallux_l">
                    <socket_frame>/bodyset/toes_l</socket_frame>
                    <radius>0.014999999999999999</radius>
                    <location>0.0562 -0.0102 0.0181</location>
                </ContactSphere>
                <ContactSphere name="othertoes_l">
                    <socket_frame>/bodyset/toes_l</socket_frame>
                    <radius>0.014999999999999999</radius>
                    <location>0.0441 -0.006 -0.0242</location>
                </ContactSphere>
                
                <ContactHalfSpace name="floor">
                    <socket_frame>/ground</socket_frame>
                    <orientation>0 0 -1.5707963267948966</orientation>
                </ContactHalfSpace>
            </objects>
			<groups />
		</ContactGeometrySet>

Maybe I was looking in the wrong place? It would be great to be able to see the contact locations from the model with 11 spheres under each foot.

Thanks,
Evan

[Ross Miller]

These are the ones I'm using recently:

Code: Select all

<ContactSphere name="heel_r">
                    <socket_frame>/bodyset/calcn_r</socket_frame>
                    <radius>0.015</radius>
                    <location>0.001307527581931796 -0.0081 0</location>
                </ContactSphere>
                <ContactSphere name="latheel_r">
                    <socket_frame>/bodyset/calcn_r</socket_frame>
                    <radius>0.015</radius>
                    <location>0.021307527581931796 -0.0081 0.025</location>
                </ContactSphere>
                <ContactSphere name="medheel_r">
                    <socket_frame>/bodyset/calcn_r</socket_frame>
                    <radius>0.015</radius>
                    <location>0.021307527581931796 -0.0081 -0.025</location>
                </ContactSphere>
                <ContactSphere name="mid_r">
                    <socket_frame>/bodyset/calcn_r</socket_frame>
                    <radius>0.015</radius>
                    <location>0.08355376 -0.0081 0.0056</location>
                </ContactSphere>
                <ContactSphere name="latmid_r">
                    <socket_frame>/bodyset/calcn_r</socket_frame>
                    <radius>0.015</radius>
                    <location>0.07215376 -0.0081 0.031</location>
                </ContactSphere>
                <ContactSphere name="medmid_r">
                    <socket_frame>/bodyset/calcn_r</socket_frame>
                    <radius>0.015</radius>
                    <location>0.09165376 -0.0081 -0.023</location>
                </ContactSphere>
                <ContactSphere name="mh1_r">
                    <socket_frame>/bodyset/calcn_r</socket_frame>
                    <radius>0.014999999999999999</radius>
                    <location>0.1820 -0.0081 -0.0230</location>
                </ContactSphere>
                <ContactSphere name="mh3_r">
                    <socket_frame>/bodyset/calcn_r</socket_frame>
                    <radius>0.014999999999999999</radius>
                    <location>0.1658 -0.0081 0.0116</location>
                </ContactSphere>
                <ContactSphere name="mh5_r">
                    <socket_frame>/bodyset/calcn_r</socket_frame>
                    <radius>0.014999999999999999</radius>
                    <location>0.1430 -0.0081 0.0370</location>
                </ContactSphere>
                <ContactSphere name="hallux_r">
                    <socket_frame>/bodyset/toes_r</socket_frame>
                    <radius>0.014999999999999999</radius>
                    <location>0.0562 -0.0102 -0.0181</location>
                </ContactSphere>
                <ContactSphere name="othertoes_r">
                    <socket_frame>/bodyset/toes_r</socket_frame>
                    <radius>0.014999999999999999</radius>
                    <location>0.0441 -0.006 0.0242</location>
                </ContactSphere>