simulPET/phantom.py

import numpy

import csv

def generatePhantom(nx,ny,N):
    im=numpy.zeros([nx,ny])
    file='/home/studen/software/build/observer_tools/fisherStudy/phantoms/phantom_0_rw_0p6.out'
    rows=[]
    with open(file) as fobj:
        d=csv.reader(fobj,delimiter=' ')
        for r in d:
            rows.append(r)
    
    px=0.5*nx
    py=0.5*ny

    phantomSize=25
    sx=phantomSize/px
    sy=phantomSize/py

    fs=0
    for r in rows:
        #ring
        if r[0]=='1':
            fr=float(r[3])
            area=fr*fr*float(r[4])
            fs+=area
            r.append(area)
        #rim
        if r[0]=='2':
            fr2=float(r[2])
            fr1=float(r[1])
            area=float(r[3])*(fr2*fr2-fr1*fr1)
            fs+=area
            r.append(area)

    print("hotspot lookup done")


    for i in numpy.arange(N):
        xe1=numpy.random.random()
        
        for r in rows:
            print(r)
            if r[0]=='1':
                q=r[5]/fs
            if r[0]=='2':
                q=r[4]/fs
            
            if xe1>q:
                xe1-=q
                continue
            if r[0]=='1':
                qr=numpy.sqrt(numpy.random.random())*float(r[3])
                qphi=2*3.14159*numpy.random.random()
                qx=float(r[1])+qr*numpy.cos(qphi)
                qy=float(r[2])+qr*numpy.sin(qphi)
            if r[0]=='2':
                qr=float(r[1])+numpy.sqrt(numpy.random.random())*(float(r[2])-float(r[1]))
                qphi=2*3.14159*numpy.random.random()
                qx=qr*numpy.cos(qphi)
                qy=qr*numpy.sin(qphi)
            
            ix=int((qx+phantomSize)/sx)
            iy=int((qy+phantomSize)/sy)
            im[ix,iy]+=1.0;
            break
        
        print("hotspot lookup done")
    
    print("particle done")
    return im