I am interested in understanding how the values in GLCM and GLSZM are computed for each gray level. I know that there are C++ functions calculate_glcm and calculate_glszm that compute these matrices. Using these functions and adding some more lines of code, can I, for example, track step-by-step how size zones are calculated for each gray level in the GLSZM? Is it possible?

The actual algortihm filling the matrices of the various feature classes is programmed into the C-extension of PyRadiomics. In the earliest versions, we did have python implementations, but due to unsatisfactory performance we have since moved to C-only calculation of the matrices.

Depending on your specific aim, there are several possibilities. If you need to influence the process or extract intermediate results, you can revisit the history of the repository, extracting the python implementations of the matrix calucaltion (though this does not support voxel-based radiomics).
Moreover, the python algorithms were slightly different to make optimum use of optimized numpy calls.

However, if you are seeking a conceptual view of how the matrices are filled, I suggest to start with the C-source code included in the repository. Here you have 2 options. Option 1 is to read the current version, which is more complex to allow for voxel-based calculation and to support both 2D and 3D input. Alternatively, you can read the C-extensions of an earlier version, which required 3D input and could therefore have less complex algorithms. The last version of PyRadiomics with this earlier version of the C-extensions can be found here.

In the cmatrices.c file, there is a function calculate_glszm which, as I understand, calculates the size zones for each gray-level and this is stored in temp_data. Later, the function fill_glszm is used to create glszm. I couldn’t locate how fill_glszm is called in the C implementation. Could you please let me know how that is done?

Yes, when comparing time required for matrix calculation in C and in Python only, there was a 300-1000 fold faster calculation in C

note: this only applies to the actual calculation of the matrices. Though overall the time gain was less dramatic, it was still quite significant (in the order of: “instead of a week it just took a couple of minutes”)