the quaternion stabilization is indeed currently missing a MAX_SUM_ERROR, so it will simply wrap around (since the quaternion is normalized) and look like in your plots.
Although this should not happen in flight with a normal quad (I have never seen that happen), you can of course force this on the ground/desk and weird things will happen...
The first thing that should be added is a configurable error integration saturation limit (MAX_SUM_ERROR).
This should to be done for each axis separately. Currently I can't think of a way to do this in quaternion formulation directly, so we probably need to convert it to 3 separate angles relative to the NED plane, limit it and convert back to a quaternion....
A fast but very hacky way would be to convert the quaternion to axis/angle, limit the angle and convert it back. But you then get coupling between the axis and start capping the sum_error for e.g. roll other axis if yaw sum_err get's too big.
Ideas for doing this properly and efficiently are welcome :-)
Secondly it would of course be good to have a proper anti-windup as well (stop the integrator when the output command is saturated), but that means that you would need to be able to extract the saturation information for each axis from motor_mixing and feed that back to stabilization.