from glob import glob
import pydicom
import json
import os.path as osp
import numpy as np

dicom_dir = "D:/3Dsyokudo"

out_name = "3D.raw"
out_base_path = "../Assets/"
dicom_file_list = glob(osp.join(dicom_dir, "*"))
dicom_file_list = sorted(dicom_file_list, key=lambda x: float(pydicom.read_file(x).SliceLocation), reverse=False)
sample_dcm1, sample_dcm2 = pydicom.read_file(dicom_file_list[0]), pydicom.read_file(dicom_file_list[1])

height, width = pydicom.read_file(dicom_file_list[0]).pixel_array.shape
depth = len(dicom_file_list)

x_scale_f, y_scale_f = sample_dcm1.PixelSpacing
z_scale_f = abs(sample_dcm1.ImagePositionPatient[2] - sample_dcm2.ImagePositionPatient[2])

stack_img = np.stack([pydicom.read_file(x).pixel_array + pydicom.read_file(x).RescaleIntercept for x in dicom_file_list], axis=2)
row_format = stack_img.transpose((1, 0, 2)).reshape(-1, order="F")

# ウィンドウ使うバージョン
window_max = 330
window_min = -270

scaled_row = row_format.copy().astype(np.float)
scaled_row[scaled_row < window_min] = window_min
scaled_row[window_max < scaled_row] = window_max

scaled_row -= np.mean(scaled_row)
scaled_row = scaled_row / (np.max(np.abs(scaled_row)) + 1e-5) * (2 ** 16)
scaled_row -= np.mean(scaled_row)
scaled_row = np.clip(scaled_row, 0, 2 ** 16).astype(np.uint16)
scaled_row.tofile(osp.join(out_base_path, out_name))

to_json = {
    "x_scale": width * x_scale_f,
    "y_scale": height * y_scale_f,
    "z_scale": depth * z_scale_f
}

with open(osp.join(out_base_path, "setting.json"), "w") as f:
    json.dump(to_json, f, indent=4)