stim_pixe_tomography/scripts/generate_voxelized_worm_phantom.py

0001 import sys
0002 import struct
0003 import math
0004
0005 # tables for super resolution (vol_work) and regular resolution (vol_result)
0006 vol_work = []
0007 vol_result = []
0008
0009 ############################################################################
0010 # size image 128x128 -> 76.464 microm -> resol x,y : 0,597375 microm      #
0011 # nb slices 128 -> 201,217 microm -> resol z : 1,57200781 microm            #
0012 # slice of interest : 11 -> 18,07 microm                                      #
0013 ############################################################################
0014
0015 sizex = 128
0016 sizey = 128
0017 sizez = 128
0018 double_sizez = 256  # for computation only
0019 resolx = 0.597375
0020 resoly = 0.597375
0021 resolz = 1.57200781
0022 superres = 8
0023
0024 slice = 11
0025
0026 # ellipsoide 1 (semi-axes, center, rotation z, value) - skin
0027 a1 = 20.61
0028 b1 = 21.42
0029 c1 = 187.82
0030 x1 = 0.0
0031 y1 = 0.0
0032 z1 = 0.0
0033 rz1 = 0.0
0034 value1 = 49.73
0035
0036 # ellipsoide 2 (semi-axes, center, rotation z, value) - body
0037 a2 = 18.61
0038 b2 = 19.01
0039 c2 = 186.64
0040 x2 = -0.39
0041 y2 = 0.0
0042 z2 = 0.0
0043 rz2 = 0.0
0044 value2 = 40.85
0045
0046 # ellipsoide 3 (semi-axes, center, rotation z, value) - core 1
0047 a3 = 1.95
0048 b3 = 3.23
0049 c3 = 4.32
0050 x3 = 1.97
0051 y3 = -7.09
0052 z3 = 18.07
0053 rz3 = 0.0
0054 value3 = 66.26
0055
0056 # ellipsoide 4 (semi-axes, center, rotation z, value) - core 2
0057 a4 = 2.08
0058 b4 = 2.46
0059 c4 = 4.32
0060 x4 = 8.27
0061 y4 = -3.15
0062 z4 = 18.07
0063 rz4 = 0.0
0064 value4 = 60.23
0065
0066 # ellipsoide 5 (semi-axes, center, rotation z, value) - intestine
0067 a5 = 3.67
0068 b5 = 16.48
0069 c5 = 28.68
0070 x5 = 1.25
0071 y5 = 0.61
0072 z5 = 0
0073 rz5 = -58.99
0074 value5 = 54.06
0075
0076 # ellipsoide 6 ((emi-axes, center, rotation z, value) - region titane
0077 a6 = 1.62
0078 b6 = 1.95
0079 c6 = 1.62
0080 x6 = 6.25
0081 y6 = 3.61
0082 z6 = 18.07
0083 rz6 = 0.0
0084 value6 = 75.14
0085
0086 # size in super resolution by voxel
0087 super_sizex = sizex * superres
0088 super_sizey = sizey * superres
0089 super_sizez = double_sizez * superres  # we will save only half of the z
0090 center_shift = sizex / 2  # for x and y axis
0091 center_shift_z = double_sizez / 2
0092
0093
0094 def make_ellipse_center(a, b, c, x, y, z, rz, value):
0095     asample = (a / resolx) * superres
0096     bsample = (b / resoly) * superres
0097     csample = (c / resolz) * superres
0098     # if (rz != 0.0):
0099     rzsample = math.radians(rz)  # angle in radians
0100     xsample = (x / resolx) * superres
0101     ysample = (y / resoly) * superres
0102     zsample = (z / resolz) * superres
0103     center_shift_sample = center_shift * superres  # translation of center for x and y axis
0104     center_shift_z_sample = center_shift_z * superres  # translation of center for z axis
0105
0106     print(asample, bsample, csample, xsample, ysample, zsample)
0107     number = 0  ## debug
0108
0109     # a loop in axes of voxels
0110     for k in range(0, super_sizez):
0111         for j in range(0, super_sizey):
0112             for i in range(0, super_sizex):
0113                 ii = i + 0.5
0114                 jj = j + 0.5
0115                 kk = k + 0.5
0116                 res = pow(((ii - xsample - center_shift_sample) * math.cos(rzsample) + (
0117                             jj - ysample - center_shift_sample) * math.sin(rzsample)), 2) / (asample * asample) + \
0118                       pow(((ii - xsample - center_shift_sample) * math.sin(rzsample) - (
0119                                   jj - ysample - center_shift_sample) * math.cos(rzsample)), 2) / (bsample * bsample) + \
0120                       ((kk - zsample - center_shift_z_sample) * (kk - zsample - center_shift_z_sample)) / (
0121                                   csample * csample)  # z-axis correction done
0122                 # print(res)
0123                 # if the voxel belongs to the ellipsoide, set the value
0124                 if (res <= 1.0):
0125                     vol_work[i + j * super_sizex + k * super_sizex * super_sizey] = value
0126                     number += 1
0127
0128     print(number)
0129
0130
0131 # initialisation
0132 def initialize():
0133     for k in range(0, double_sizez):
0134         for j in range(0, sizey):
0135             for i in range(0, sizex):
0136                 vol_result.append(0.0)
0137     for k in range(0, super_sizez):
0138         for j in range(0, super_sizey):
0139             for i in range(0, super_sizex):
0140                 vol_work.append(0.0)
0141
0142
0143 # Calculate the vol_result based on vol_work
0144 def undersample():
0145     x = 0
0146     y = 0
0147     z = 0
0148     for k in range(0, super_sizez, superres):
0149         for j in range(0, super_sizey, superres):
0150             for i in range(0, super_sizex, superres):
0151                 # on se place sur v1 et on recupere les valeurs de densite des 8 voxels du voisinage qui vont correspondre a 1 voxel de l'image finale
0152                 # v1 = vol_work[i+j*super_sizex+k*(super_sizex*super_sizey)]
0153                 # v2 = vol_work[i+j*super_sizex+k*(super_sizex*super_sizey) + 1]
0154                 # v3 = vol_work[i+j*super_sizex+k*(super_sizex*super_sizey) + super_sizex]
0155                 # v4 = vol_work[i+j*super_sizex+k*(super_sizex*super_sizey) + super_sizex + 1]
0156                 # v5 = vol_work[i+j*super_sizex+k*(super_sizex*super_sizey) + super_sizex*super_sizey]
0157                 # v6 = vol_work[i+j*super_sizex+k*(super_sizex*super_sizey) + 1 + super_sizex*super_sizey]
0158                 # v7 = vol_work[i+j*super_sizex+k*(super_sizex*super_sizey) + super_sizex + super_sizex*super_sizey]
0159                 # v8 = vol_work[i+j*super_sizex+k*(super_sizex*super_sizey) + super_sizex + 1 + super_sizex*super_sizey]
0160                 # vol_result[x+y*sizex+z*(sizex*sizey)] = (v1+v2+v3+v4+v5+v6+v7+v8)/8.0
0161                 # print ("***",i,j,k)
0162                 total = 0
0163                 for kk in range(0, superres):
0164                     for jj in range(0, superres):
0165                         for ii in range(0, superres):
0166                             total = total + vol_work[
0167                                 i + ii + (j + jj) * super_sizex + (k + kk) * (super_sizex * super_sizey)]
0168                 vol_result[x + y * sizex + z * (sizex * sizey)] = total / (superres * superres * superres)
0169                 # print ("###",vol_result[x+y*sizex+z*(sizex*sizey)])
0170                 x += 1
0171             y += 1
0172             x = 0
0173         z += 1
0174         y = 0
0175
0176
0177 # Save the total volume (including the negative parts of the ellipsoids)
0178 def save_whole_work(file):
0179     fd = open(file, "wb")
0180     for i in range(0, len(vol_work)):
0181         fd.write(struct.pack("f", vol_work[i]))
0182
0183
0184 # Save half the volume (including only the positive parts of the ellipsoids)
0185 def save_half_work(file):
0186     fd = open(file, "wb")
0187     for i in range(int(len(vol_work) / 2), len(vol_work)):
0188         fd.write(struct.pack("f", vol_work[i]))
0189
0190
0191 # save one slice in super resolution vol_work
0192 #  0<=slice<super_sizez
0193 def save_workslice(file, slice):
0194     fd = open(file, "wb")
0195     for j in range(0, super_sizey):
0196         for i in range(0, super_sizex):
0197             fd.write(struct.pack("f", vol_work[i + j * super_sizex + slice * super_sizex * super_sizey]))
0198
0199
0200 # Save the total result volume (including the negative parts of the ellipsoids)
0201 def save_whole_result(file):
0202     fd = open(file, "wb")
0203     for i in range(0, len(vol_result)):
0204         fd.write(struct.pack("f", vol_result[i]))
0205
0206
0207 # Save half of the result volume (including only the positive parts of the ellipsoids)
0208 def save_half_result(file):
0209     fd = open(file, "wb")
0210     for i in range(int(len(vol_result) / 2), len(vol_result)):
0211         fd.write(struct.pack("f", vol_result[i]))
0212
0213
0214 # save one slie in regular resolution vol_result
0215 #  0<=slice<128
0216 def save_slice(file, slice):
0217     fd = open(file, "wb")
0218     for j in range(0, sizey):
0219         for i in range(0, sizex):
0220             fd.write(struct.pack("f", vol_result[i + j * sizex + slice * sizex * sizey]))
0221
0222
0223 print("------intialisation------")
0224 initialize()
0225 print("------end intialisation------")
0226
0227 print("------begin ellipse 1------------")
0228 make_ellipse_center(a1, b1, c1, x1, y1, z1, rz1, value1)
0229 print("------end ellipse 1--------------")
0230
0231 print("------begin ellipse 2------------")
0232 make_ellipse_center(a2, b2, c2, x2, y2, z2, rz2, value2)
0233 print("------end ellipse 2--------------")
0234
0235 print("------begin ellipse 3------------")
0236 make_ellipse_center(a3, b3, c3, x3, y3, z3, rz3, value3)
0237 print("------end ellipse 3--------------")
0238
0239 print("------begin ellipse 4------------")
0240 make_ellipse_center(a4, b4, c4, x4, y4, z4, rz4, value4)
0241 print("------end ellipse 4--------------")
0242
0243 print("------begin ellipse 5------------")
0244 make_ellipse_center(a5, b5, c5, x5, y5, z5, rz5, value5)
0245 print("------end ellipse 5--------------")
0246
0247 print("------begin ellipse 6------------")
0248 make_ellipse_center(a6, b6, c6, x6, y6, z6, rz6, value6)
0249 print("------end ellipse 6--------------")
0250
0251 print("------begin undersample------")
0252 undersample()
0253 print("------end undersample------")
0254
0255 print("The length of vol_work: ", len(vol_work))
0256 print("The length of vol_result: ", len(vol_result))
0257
0258 # print ("------begin save slice 139 ------")
0259 save_slice("./slice_139.dat", 139)
0260 # print ("------end save slice 139 ------")
0261
0262
0263 save_half_result("./vol_result.dat")
0264 # save_whole_work("./vol_work.dat")
0265 print("------end save work------")