The squashed feet thing is exactly what I am seeing. It is as if the lower bones have been compressed upwards. The bones are to blame...
To test, I assigned ALL the feet bones to lower leg. This let the feet show correctly, but the knee joint was wrong. It also allowed the feet to touch the ground. This confirms that the Y scaling is hitting the target OK. Assign all the verts to 'Thigh' bones, and the legs appear OK...but he walks like a scarecrow

The file has a date on it of 24/5/07 as last modified, so I am guessing this is when you mailed it over.

Code:
# -----------------------------------------------------------------------------------
#    Computes Euler angles from quaternions.                                        |
# -----------------------------------------------------------------------------------
def computeeulers( quatfloats ) :

    nfloats                = len( quatfloats )
    nvecs                  = nfloats / 4              
    eulers                 = array.array( 'f' )

    for ivec in range( nvecs ) :
        idx                = ivec * 4
        q1                 = quatfloats[idx+0]
        q2                 = quatfloats[idx+1]
        q3                 = quatfloats[idx+2]
        q4                 = quatfloats[idx+3]
                           
        Q11                = 1 - 2 * ( q2*q2 + q3*q3 )
        Q12                =     2 * ( q1*q2 - q3*q4 )
        Q13                =     2 * ( q1*q3 + q2*q4 )
        Q21                =     2 * ( q1*q2 + q3*q4 )
        Q22                = 1 - 2 * ( q1*q1 + q3*q3 )
        Q23                =     2 * ( q2*q3 - q1*q4 )
        Q31                =     2 * ( q1*q3 - q2*q4 )
        Q32                =     2 * ( q2*q3 + q1*q4 )
        Q33                = 1 - 2 * ( q1*q1 + q2*q2 )
                           
#        # 231 method.      
#        spsi               = Q21
#        cpsi               = ( 1 - spsi*spsi )**0.5
#        stheta             = -Q31 / cpsi
#        ctheta             = +Q11 / cpsi
#        sphi               = -Q23 / cpsi
#        cphi               = +Q22 / cpsi
#                           
#        phi_rad            = math.atan2( sphi,   cphi )       
#        theta_rad          = math.atan2( stheta, ctheta )     
#        psi_rad            = math.atan2( spsi,   cpsi )       
#                           
#        # Try.             
#        eulers.append( phi_rad )
#        eulers.append( -theta_rad )
#        eulers.append( -psi_rad )
#                           
#        # Hey! this way  works!

        # GrumpyOldMan's code.
        sint               = 2 * ( q2 * q4 - q1 * q3 )
        cost_tmp           = 1.0 - sint*sint
        if abs( cost_tmp ) > 0.0001 :
            cost           = cost_tmp**0.5
        else :
            cost           = 0.0

        if abs( cost ) > 0.01 :
            sinv           = 2 * ( q2 * q3 + q1 * q4 )     / cost
            cosv           = ( 1 - 2 * ( q1*q1 + q2*q2 ) ) / cost
            sinf           = 2 * ( q1 * q2 + q3 * q4 )     / cost
            cosf           = ( 1 - 2 * ( q2*q2 + q3*q3 ) ) / cost
        else:
            sinv           = 2 * ( q1 * q4 - q2 * q3 )     
            cosv           = 1 - 2 * ( q1*q1 + q3*q3 )  
            sinf           = 0.0
            cosf           = 1.0

        roll               = math.atan2( sinv, cosv )
        pitch              = math.atan2( sint, cost )
        yaw                = math.atan2( sinf, cosf )

        eulers.append( roll )
        eulers.append(-pitch )
        eulers.append(-yaw )  

    return eulers


# -----------------------------------------------------------------------------------
#    Computes quaternions from Milkshape's x, y, z Euler angles.                    |
# -----------------------------------------------------------------------------------
def computequats( eulers ) :

    nfloats                = len( eulers )
    nvecs                  = nfloats / 3 
    quatfloats             = array.array( 'f' )

    for ivec in range( nvecs ) :
        idx                = ivec * 3
        phi_rad            = +eulers[idx+0]                # NOTE: These are the signs that worked in animmerge.py
        theta_rad          = -eulers[idx+1]
        psi_rad            = -eulers[idx+2]
                         
        sphi               = math.sin( phi_rad )
        cphi               = math.cos( phi_rad )
        stheta             = math.sin( theta_rad )
        ctheta             = math.cos( theta_rad )
        spsi               = math.sin( psi_rad )
        cpsi               = math.cos( psi_rad )
                         
#        # 231 method (or 132 if you follow that naming convention).
#        R11                = ctheta * cpsi       
#        R12                = -ctheta * spsi * cphi + stheta * sphi
#        R13                =  ctheta * spsi * sphi + stheta * cphi
#        R21                = spsi
#        R22                = cpsi * cphi
#        R23                = -cpsi * sphi
#        R31                = -stheta * cpsi
#        R32                =  stheta * spsi * cphi + ctheta * sphi
#        R33                = -stheta * spsi * sphi + ctheta * cphi
                         
        # 321 method (or 123 if you follow that naming convention).
        R11                = cpsi * ctheta       
        R12                = cpsi * stheta * sphi - spsi * cphi 
        R13                = cpsi * stheta * cphi + spsi * sphi
        R21                = spsi * ctheta
        R22                = spsi * stheta * sphi + cpsi * cphi
        R23                = spsi * stheta * cphi - cpsi * sphi
        R31                = -stheta
        R32                = ctheta * sphi
        R33                = ctheta * cphi
                         
        q4                 = 0.5 * ( 1 + R11 + R22 + R33 )**0.5
        q1                 = 0.25 * ( R32 - R23 ) / q4
        q2                 = 0.25 * ( R13 - R31 ) / q4
        q3                 = 0.25 * ( R21 - R12 ) / q4

        quatfloats.append( q1 )
        quatfloats.append( q2 )
        quatfloats.append( q3 )
        quatfloats.append( q4 )

    return quatfloats
I am guessing this is the bit you wanted from the code

The animation library has a created date of 16th May in the header...but I guess if youare like me with this sort of thing, that's not necessarily when you last changed it! The e-mail you sent me confirmed that this was the version you did AFTER you had changed the calcs for the extreme rotations to GOM's version. The version you sent me on the 22nd had the scaling function added.

As a last test, I will remove all the tools and bits and pieces, and make 100% I have only the correct versions in place, and double check. At each stage, I will run a test on the model to see if the legs are correct. This should help to pin down if it's an error that was always there, or something that crept in when a feature was added. I will PM you later on when I have some more data for you to work with