Percentage of tissues within resampled voxels

I need to simulate MRI signals and to quantify the contribution of different tissues (materials) to the signal generated by each single MRI voxel. This is my pipeline to calculate the partial tissue volumes, but I hope you can propose faster and more efficient solutions.

Step 1: Extraction of the edges of MRI voxels
I executed a dummy EM simulation with a pretty dense grid to resolve all tissues over the region of interest (ROI). Then using the resampling function in the Python API, I resampled the fields on a grid providing the desired number of MRI voxels (e.g. 64x64x64 voxels) and I extracted the new grid. This permits me to have information about dimensions, edge coordinates and voxel centers of all my MRI voxels.

Step2: Creation of Masks for Each Tissue in the ROI
I created a masks (mask_tissue) for each tissue within the ROI using the masking filter tool.

Step3: Creation of masks for each single MRI voxel
Using the modeling functions in the python API, I created an entity for each MRI voxel. With this entity object I can create a mask (mask_voxel) using the masking filter function.

Step4: Calculation of tissue partial volumes within MRI voxels
For each voxel and for each tissue I calculate the product masks_tissue x mask_voxel. The partial volume of the ith tissue into the jth MRI voxel is: sum(masks_tissue x mask_voxel)/sum(mask_voxel).

The pipeline works fine but is very slow. Would you suggest alternative options?

Do you know which of the four steps is slowing down the pipeline?
Because you could always choose to go for a more dedicated programming language rather than python to do part of the processing.