stim_pixe_tomography/scripts/generate_voxelized_sphere_phantom.py

0001 import sys
0002 import struct
0003 import math
0004 import numpy as np
0005
0006 # lists for super resolution (vol_work) and regular resolution (vol_result)
0007 vol_work = []
0008 vol_result = []
0009
0010 ############################################################################
0011 # size image 128x128 -> 500 microm -> resol x,y : 3.90625 microm      #
0012 # nb slices 128 -> 500 microm -> resol z : 3.90625 microm            #
0013 ############################################################################
0014
0015 sizex = 128
0016 sizey = 128
0017 sizez = 128
0018 resolx = 3.90625
0019 resoly = 3.90625
0020 resolz = 3.90625
0021 superres = 8  # factor of super resolution
0022
0023 # sphere 1 (radius) - outer sphere
0024 r1 = 196
0025
0026 # sphere 2 (radius) - inner sphere
0027 r2 = 171
0028
0029 # density values for sphere 1 and sphere 2
0030 type = 2  # type for constructing a STIM or PIXE phantom
0031 # type = 1, STIM phantom, density value for STIM in 0.01 g/cm3
0032 # type =2, PIXE phantom, density value for PIXE in 0.000001 g/cm3, namely microgram/cm3
0033 value1 = 0
0034 value2 = 0
0035
0036 if type == 1:
0037     value1 = 108
0038     value2 = 0
0039 elif type == 2:
0040     value1 = 54000
0041     value2 = 0
0042
0043 # center of two spheres
0044 x0 = 0.0
0045 y0 = 0.0
0046 z0 = -1.953125
0047
0048 # size in super resolution by voxel
0049 super_sizex = sizex * superres
0050 super_sizey = sizey * superres
0051 super_sizez = sizez * superres
0052 center_shift = sizex / 2  # translation of half scan
0053
0054
0055 def make_sphere_center(r, x, y, z, value):
0056     rsample = (r / resolx) * superres  # radius in super resolution by voxel
0057     xsample = (x / resolx) * superres
0058     ysample = (y / resoly) * superres
0059     zsample = (z / resolz) * superres
0060     center_shift_sample = center_shift * superres  # translation of the center for x (i) and y (j) axis
0061
0062     print(rsample, xsample, ysample, zsample)
0063     number = 0
0064
0065     for k in range(0, super_sizez):
0066         for j in range(0, super_sizey):
0067             for i in range(0, super_sizex):
0068                 ii = i + 0.5
0069                 jj = j + 0.5
0070                 kk = k + 0.5
0071                 res = pow((ii - xsample - center_shift_sample), 2) / (rsample * rsample) + \
0072                       pow((jj - ysample - center_shift_sample), 2) / (rsample * rsample) + \
0073                       pow((kk - zsample - center_shift_sample), 2) / (rsample * rsample)  # z-axis correction done
0074                 # if the point (ii, jj, kk) is in the sphere, we attribute the voxel (i,j, k) value
0075                 if (res <= 1.0):
0076                     vol_work[i + j * super_sizex + k * super_sizex * super_sizey] = value
0077                     number += 1
0078
0079     print(number)
0080
0081
0082 # initialisation for two tables vol_result (128*128*128), vol_work (128*superres)*(128*superres)*(128*superres)
0083 def initialize():
0084     for k in range(0, sizez):
0085         for j in range(0, sizey):
0086             for i in range(0, sizex):
0087                 vol_result.append(0.0)
0088     for k in range(0, super_sizez):
0089         for j in range(0, super_sizey):
0090             for i in range(0, super_sizex):
0091                 vol_work.append(0.0)
0092
0093
0094 # Calculate the vol_result based on vol_work
0095 def undersample():
0096     x = 0
0097     y = 0
0098     z = 0
0099     for k in range(0, super_sizez, superres):
0100         for j in range(0, super_sizey, superres):
0101             for i in range(0, super_sizex, superres):
0102                 # print ("***",i,j,k)
0103                 total = 0
0104                 for kk in range(0, superres):
0105                     for jj in range(0, superres):
0106                         for ii in range(0, superres):
0107                             total = total + vol_work[
0108                                 i + ii + (j + jj) * super_sizex + (k + kk) * (super_sizex * super_sizey)]
0109                 vol_result[x + y * sizex + z * (sizex * sizey)] = total / (superres * superres * superres)
0110                 # print ("###",vol_result[x+y*sizex+z*(sizex*sizey)])
0111                 x += 1
0112             y += 1
0113             x = 0
0114         z += 1
0115         y = 0
0116
0117
0118 # save the total volume of super resolution
0119 def save_whole_work(file):
0120     fd = open(file, "wb")
0121     for i in range(0, len(vol_work)):
0122         fd.write(struct.pack("f", vol_work[i]))
0123     fd.close()
0124
0125
0126 # save one slice in super resolution vol_work
0127 #  0<=slice<super_sizez (128*8=1024)
0128 def save_workslice(file, slice):
0129     fd = open(file, "wb")
0130     for j in range(0, super_sizey):
0131         for i in range(0, super_sizex):
0132             fd.write(struct.pack("f", vol_work[i + j * super_sizex + slice * super_sizex * super_sizey]))
0133     fd.close()
0134
0135
0136 # save the total result volume
0137 def save_whole_result(file):
0138     fd = open(file, "wb")
0139     for i in range(0, len(vol_result)):
0140         fd.write(struct.pack("f", vol_result[i]))
0141     fd.close()
0142
0143
0144 # save one slie in regular resolution vol_result
0145 #  0<=slice<128
0146 def save_slice(file, slice):
0147     fd = open(file, "wb")
0148     for j in range(0, sizey):
0149         for i in range(0, sizex):
0150             fd.write(struct.pack("f", vol_result[i + j * sizex + slice * sizex * sizey]))
0151
0152     fd.close()
0153
0154
0155 print("------intialisation------")
0156 initialize()
0157 print("------end intialisation------")
0158
0159 print("------begin sphere 1------")
0160 make_sphere_center(r1, x0, y0, z0, value1)
0161 print("------end sphere 1------")
0162
0163 print("------begin sphere 2------")
0164 make_sphere_center(r2, x0, y0, z0, value2)
0165 print("------end sphere 2------")
0166
0167 print("------begin undersample------")
0168 undersample()
0169 print("------end undersample------")
0170
0171 print("The length of vol_work: ", len(vol_work))
0172 print("The length of vol_result: ", len(vol_result))
0173
0174 # print ("------begin save slice 63 ------")
0175 # save_slice("./slice_63.dat",63)
0176 # print ("------end save slice 63 ------")
0177
0178 save_whole_result("./vol_result.dat")
0179 save_whole_work("./vol_work.dat")
0180 print("------end save work------")