por_dis_coor

import numpy as np
import matplotlib.pyplot as plt
import os

def calculate_porosity_per_slice(raw_file, sample_size, z_min):
    
    #calculates the porosity of each slice along the length of a sample from a raw file.
    #raw_file: Path to the raw file.
    #sample_size: Size of the 3D sample.
    #z_min: Starting z-coordinate of the sample.

    with open(raw_file, 'rb') as f:
        data = np.fromfile(f, dtype=np.uint8)
        data = data.reshape(sample_size)

    #create an empty list to store the porosities
    porosities = []

    #calculate the porosity for each slice along the length (z-axis)
    for z in range(sample_size[0]):  # Adjusted for (z, y, x) format
        slice_data = data[z, :, :]  # Adjusted for (z, y, x) format

        #count the number of pore voxels (value 0) and solid voxels (value 255)
        pore_voxels = np.count_nonzero(slice_data == 0)
        solid_voxels = np.count_nonzero(slice_data == 255)

        #calculate porosity
        porosity = pore_voxels / (pore_voxels + solid_voxels)
        porosities.append(porosity)

    return porosities

#define input parameters
sample_size = (350, 200, 200)  # Adjusted for (z, y, x) format
output_prefix = "sample"

#list of sample numbers that should be processed
high_value_samples = [1, 4, 6, 10, 11, 16, 17, 20, 21, 22, 23]

#calculate the porosity for each sample in high_value_samples
for i in high_value_samples:
    raw_file = f"{output_prefix}_{i}.raw"
    coordinates_file = f"{output_prefix}_coordinates.txt"
    if os.path.exists(raw_file) and os.path.exists(coordinates_file):
        #read the coordinates from the text file
        with open(coordinates_file, 'r') as f:
            coordinates = [eval(line.strip()) for line in f.readlines()]
        x_min, y_min, z_min = coordinates[i]  # Adjusted for (x, y, z) format

        porosities = calculate_porosity_per_slice(raw_file, sample_size, z_min)

        #convert slice number to length in micrometers, starting from the top
        length_in_micrometers = [(z + z_min) * 1.5 for z in range(len(porosities))]

        #plot porosity distribution along the length for each sample
        plt.figure()
        plt.plot(length_in_micrometers, porosities)
        plt.xlabel('Length (micrometers)')
        plt.ylabel('Porosity')
        plt.title(f'Porosity Distribution along the Length for Sample {i}')
        plt.show()