dkazanc / TomoPhantom

Software to generate 2D/3D/4D analytical phantoms and their Radon transforms (parallel beam) for image processing
https://dkazanc.github.io/TomoPhantom/
Apache License 2.0
116 stars 53 forks source link

TomoP3D.object on Windows crashes for large phantom sizes >1500^3 #112

Closed dkazanc closed 2 years ago

dkazanc commented 2 years ago

As reported by Oliver Helps

Kernel crashes in the TomoP3D.object function with anything above roughly 1100 ("Windows fatal exception: Windows fatal exception: access violation access violation)

See below the code I have used:

import numpy as np
import time
import tomophantom
from tomophantom import TomoP3D
from tomophantom.TomoP3D import Objects3D
import matplotlib.pyplot as plt

Horiz_det = 1500
Vert_det = Horiz_det
sliceSel = int(0.5*Horiz_det)

angles_num = int(0.5*np.pi*Horiz_det); # angles number
angles = np.linspace(0.0,179.9,angles_num,dtype='float32') # in degrees
angles_rad = angles*(np.pi/180.0)

#generate base cylinder
obj3D1 = {'Obj': Objects3D.ELLIPCYLINDER,
      'C0' : 0.5,#amplitude - replace for linear attenuation coefficient at the given energy bin
      'x0' : 0,
      'y0' : 0,
      'z0' : 0.0,
      'a'  : 0.9,
      'b'  : 0.9,
      'c'  : 0.5,
      'phi1'  : 0.0}

obj3D2 = {'Obj': Objects3D.ELLIPCYLINDER,
      'C0' : -0.5,#amplitude - replace for linear attenuation coefficient at the given energy bin
      'x0' : 0,
      'y0' : 0,
      'z0' : 0.0,
      'a'  : 0.3,
      'b'  : 0.3,
      'c'  : 0.2,
      'phi1'  : 0.0}

myObjects = [obj3D1, obj3D2]

t1 = time.perf_counter()
print('Generating 3D object of size ' + str(Horiz_det))
Object3D = TomoP3D.Object(Horiz_det, myObjects)
t2 = time.perf_counter()
finish_time = t2-t1
print('time to generate 3D object: ' + str(finish_time))

#ignore intenties but this is the shape of the reference sample that has been made
plt.figure(dpi=300)
plt.imshow(Object3D[sliceSel,:,:],vmin=0, vmax=1, cmap='gray')
plt.title('3D Object, axial view')

print('Generating sino data')
t1 = time.perf_counter()
tempProjData3D = TomoP3D.ObjectSino(Horiz_det, Horiz_det, Vert_det, angles, myObjects)  
t2 = time.perf_counter()
finish_time = t2-t1
print('time to generate sino data: ' + str(finish_time))
dkazanc commented 2 years ago

The issue has been confirmed on a Windows machine. Could be related to this and this