CCD itself will work like a charm.
If you are working on 3D, first find the rotation you should do in each axis without applying bone limits.
After that, the approach
expectedRotation.X = min(expectedRotation.X, maxLimit.X)
expectedRotation.X = max(expectedRotation.X, minLimit.X)
expectedRotation.Y = min(expectedRotation.Y, maxLimit.Y)
expectedRotation.Y = max(expectedRotation.Y, minLimit.Y)
expectedRotation.Z = min(expectedRotation.Z, maxLimit.Z)
expectedRotation.Z = max(expectedRotation.Z, minLimit.Z)
is wrong. Because, if you can't move in one of the axis further, the other two axis will keep moving and you will get weird results.
The fix:
If any of the 3 axis mismatches with the limit constraints, you must not change the rotation at all.
First convert all the angles into degrees in -180 to 180 format. Then the following will do
vector3df angleDifference = expectedRotation - baseRotation; //baseRotation is just the initial rotation from which the bone limits are calculated.
if(angleDifference.X < boneLimits.minRotation.X || angleDifference.Y < boneLimits.minRotation.Y || angleDifference.Z < boneLimits.minRotation.Z || angleDifference.X > boneLimits.maxRotation.X || angleDifference.Y > boneLimits.maxRotation.Y || angleDifference.Z > boneLimits.maxRotation.Z)
return currentRotation;
return expectedRotation;