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