I am trying to obtain the grayvalues attributed to a specific material ID so I can then create a histogram to visualize the distribution of grayvalues.

I tried running the code found here (https://geodict-userguide.math2market.de/2025/automation_api-results.html?anchor=getvolumefield) in a GeoApp.



However, I received this error.


The logic behind running this code was to get all the statistics attributed to a specific Material ID, and then index out the grayvalues.

Hi,

I assume, your gray value image does not have index 1. If only one image is loaded it has always index 0. If more images are loaded, have a look at the Color By combo-box in the Volume Field tab. Counting starts with 0. The fields have the order in which they are loaded to GeoDict.


Another thing I noted in your code example, is that the indentation does not seem to make sense within the loop. In the example in the User Guide it was meant like the following. I marked the significant differences to your code with "!!!".

VolumeInfo  = gd.getVolumeFieldsInfo() # get list of dictionaries of loaded Volume Fields and assign it to variable VolumeInfo
nx, ny, nz  = gd.getVolDimensions()    # get the number of voxels in X-, Y- and Z-direction and assign the numbers to variables
Structure   = gd.getStructure()        # assign 3D numpy array describing the loaded structure to the variable Structure

# !!! In the following command ensure to have the right index for your gray value field !!!
Name        = VolumeInfo[0]['name']    # assign the name of the volume field to the variable Name.
VolumeField = gd.getVolumeField(Name)  # assign the numpy array describing the volume field to the variable VolumeField
Statistic   = []                       # Create empty list to store the statistical values

for k in range(nz): # loop over all Z-layers
 value_sum   = 0    # creating variable value_sum to sum up the result values
 value_count = 0    # creating variable value_count to count the summands

 for j in range(ny):  # loop over all Y-layers

  for i in range(nx): # loop over all X-layers

    if Structure[i][j][k] == 0:                       # condition: if the kth Z-value in the jth Y-column of the ith X-layeris pore space (ID 0), the following indented section is executed
       value_sum   = value_sum + VolumeField[i][j][k] # add all pore space result values of the kth Z-layer to the sum value_sum
       value_count = value_count + 1                  # count the summands of value_sum

 # !!! The following section should have the right indentation according to what should be plotted. Since we want to have one value for each Z-layer it belongs to the outer loop (for k in range nz) and thus needs to have the same indentation level as the head of the second loop, here 1 space. !!!
 meanVal = value_sum / value_count # compute mean value of all pore space result values in the kth Z-layer and assign it to the variable meanVal
 Statistic.append(meanVal)         # append the mean value of the Z-layer to the Statistic list

# !!! The following section needs to be outside of the loop to plot all data in one plot !!!
gDlg  = gd.makeGraphDialog()                  # create a graph dialog object
graph = gDlg.addGraph(Name, "Z layers", Name) # add a graph object with the name of the volume field as title and Y-axis legend and Z-layers as X-axis legend 
Z_layers = list(range(1, nz + 1))             # writes the Z-layer numbers 1, 2, …, nz-1, nz into a list named Z-layers
graph.addData(Z_layers, Statistic, Name)      # add a single dataset with the Z-layers as X-values, the mean result values as Y-values and the name of the volume field as legend to this graph
gDlg.run()                                    # show graph dialog

This will output something like this:

Hello,

Thank you for the helpful updates.

I was able to create histograms. However, for both Material ID 1and 2, the histogram outputs were the same (see below image).



The graph output a count of voxels equivalent to the total amount of voxels present in the structure, with an attributed grayvalue of 0.

The code I ran can be seen below.

.

I switched the values between 1 and 2 on line 16 (if Structure [j][k] ==2:).

Is there something I'm missing in regards to specifying the material ID or referencing the grayvalues? Thank you.

Hi,

your code worked more or less for me. In your plot it seems, as if no solid material ID was found. Did you load the segmented structure file with material ID 1 and 2?

By the way, it would be great if you would copy paste your code here instead of using screenshots. Then, it would be easier for me to reproduce the problem ;) Use the "#" symbol above the editor to enter the code beginning and ending markers

However I did not understand, why you want to divide the X-values by 2? Then the histogram gets moved to the left. I assume, that you wanted to account for the problem of the different lengths of the arrays count and bin_edges. For that I think it would be better to just remove the last entry of bin_edges. Even better would be to use matplotlib instead to plot the histogram. You can find an example for a very basic plot below. Matplotlib offers many options to change the plotting style as desired.

Another improvement you find in my code example below is to get the gray values with numpy.where instead of a for loop. This is much faster for large gray value images.

import numpy as np
import matplotlib.pyplot as plt

VolumeInfo  = gd.getVolumeFieldsInfo() # get list of dictionaries of loaded Volume Fields and assign it to variable VolumeInfo
nx, ny, nz  = gd.getVolDimensions()    # get the number of voxels in X-, Y- and Z-direction and assign the numbers to variables
Structure   = gd.getStructure()        # assign 3D numpy array describing the loaded structure to the variable Structure

Name        = VolumeInfo[0]['name']    # assign the name of the volume field to the variable Name.
VolumeField = gd.getVolumeField(Name)  # assign the numpy array describing the volume field to the variable VolumeField
Statistic   = []                       # Create empty list to store the statistical values

grayvalues = np.empty((nx,ny,nz),dtype=np.float32) # create 3D numpy array with the same shape as the gray value image
material_ID = np.where(Structure[:,:,:]==1)        # find out where the structure has ID 1
grayvalues[material_ID]=VolumeField[material_ID]   # take the gray values where the structure has ID 1 and save them in the new field

flat      = grayvalues.flatten()

# ---------- plot with graph dialog GeoPy API ------------------
gDlg      = gd.makeGraphDialog()                   
graph     = gDlg.addGraph('Histogram', Name, "Voxel Count")
bin_count = 256
counts, bin_edges = np.histogram(flat,bins = bin_count, range=(0,255))

x_values = bin_edges[1:-1]   # crop the first value, because the number of zeros is not interesting, crop the last value to match the number of entries in counts
y_values = counts[1:]        # crop the first entry, because the number of zeros is not interesting
x_list   = x_values.tolist()
y_list   = y_values.tolist()

graph.addData(x_list,y_list, Name)     
gDlg.run()                                   

# ---------- plot with matplotlib ------------------
plt.hist(flat, bins=bin_edges)
plt.title('Histogram')
plt.xlabel('Gray Values')
plt.ylabel('Voxel Count')
plt.savefig('HistogramPlot.png')