As to the significant differences in first order features, see the explanation above.
As to the changes in volume, this is mainly due to the different manner of calculating a surface area. I don’t know what heterogeneityCAD uses, but pyradiomics implements a marching cubes algorithm, creating a surface mesh and summing up the surface area of the individual triangles.
A different approach is to sum the surface area of the border faces of the voxels, but this overestimates surface area, especially in small volumes.
The differences in many of the other features are also due to this difference, as the formulas are dependent on the surface area value.
A possible exception to this is Compactness1; sometimes the exponent of A is flipped (i.e. 2/3 instead of 3/2). Pyradiomics has the 3/2 exponent, as this makes the formula dimensionless.
Differences in GLCM (and most likely the other texture matrices as well) are most likely due to a different method of gray level discretization (similar to the differences in Uniformity and Entropy).