import SimpleITK import config import os import re import numpy import sklearn.cluster import fitData import getData import geometry def loadTime(r,xsetup): tempDir=config.getTempDir(xsetup) code=config.getCode(r,xsetup) timeFile=os.path.join(tempDir,code,f'{code}_Dummy.csv') if not os.path.isfile(timeFile): timeFile=os.path.join(tempDir,code,f'{code}_Dummy.mcsv') with open(timeFile,'r') as f: lines=[re.sub('\n','',x) for x in f.readlines()] lines=[x for x in lines if x[0]!='#'] v=[[float(x) for x in y.split(',')] for y in lines] t=numpy.array([x[0] for x in v]) dt=numpy.array([x[1] for x in v]) #convert to seconds from miliseconds t*=1e-3 #convert to seconds from miliseconds dt*=1e-3 return t,dt def loadData(r,xsetup,returnGeometry=False): #load data from nrrd t,dt=loadTime(r,xsetup) c1=len(t) nodes=[config.getNodeName(r,xsetup,'NM',i) for i in range(0,c1)] files=[f'{x}.nrrd' for x in nodes] files=[os.path.join(config.getLocalDir(r,xsetup),x) for x in files] filesPresent=[os.path.isfile(x) for x in files] #possible side exit when missing files are encountered xdata=[SimpleITK.ReadImage(x) for x in files] geo=geometry.getGeometry(xdata[0]) xdata=[SimpleITK.GetArrayFromImage(x) for x in xdata] #create new array to hold all data data=numpy.zeros((*xdata[0].shape,len(xdata))) for i in range(len(xdata)): data[...,i]=numpy.array(xdata[i])/dt[i] if returnGeometry: return data,geo return data def getTACAtPixels(data,loc): #data is 4D array, loc are indices as returned by numpy.nonzero() #return nxm array where n is number of time points and m is number of locations #to get TAC for i-th location, do v[:,i] loc1=[loc+(numpy.array([i,i,i,i]),) for i in range(data.shape[3])] v=[data[x] for x in loc1] return numpy.vstack(v) def loadCT(r,xsetup,returnGeometry=False): file='{}.nrrd'.format(config.getNodeName(r,xsetup,'CT')) file=getData.getLocalPath(r,xsetup,file) xd=SimpleITK.ReadImage(file) geo=geometry.getGeometry(xd) xd=SimpleITK.GetArrayFromImage(xd) if returnGeometry: return xd,geo return xd def saveCenters(r,xsetup,data=None,ir=0): #if not data: spect,gSPECT=loadData(r,xsetup,returnGeometry=True) ct,gCT=loadCT(r,xsetup,returnGeometry=True) A=spect.reshape(-1,spect.shape[3]) nclass=xsetup['nclass'][0] #kmeans0 = sklearn.cluster.KMeans(n_clusters=k, random_state=0, n_init="auto").fit(A) #cmeans = sklearn.mixture.GaussianMixture(n_components=k, random_state=0, n_init=1).fit(A) kmeans = sklearn.cluster.BisectingKMeans(n_clusters=nclass, random_state=0, n_init=1).fit(A) centers=kmeans.cluster_centers_ u=kmeans.labels_ u=u.reshape(spect.shape[0:3]) print(u.shape) code=config.getCode(r,xsetup) for i in range(nclass): #ui=(u==i)*numpy.ones(u.shape) #file=getData.getLocalPath(r,xsetup,config.getPattern('centerWeight',code=code,nclass=nclass,ir=ir,ic=i)) #img=SimpleITK.GetImageFromArray(ui) #SimpleITK.WriteImage(img, file) cFile=getData.getLocalPath(r,xsetup,config.getPattern('center',code=code,nclass=nclass,ir=ir,ic=i)) numpy.savetxt(cFile,centers[i:i+1,:],delimiter=',') #write center map as NRRD file in spect geometry: file=getData.getLocalPath(r,xsetup,config.getPattern('centerNRRD',code=code,nclass=nclass,ir=ir,qaName='SPECT')) img=SimpleITK.GetImageFromArray(u) img.SetOrigin(gSPECT.origin) img.SetSpacing(gSPECT.spacing) img.SetDirection(numpy.ravel(gSPECT.direction)) SimpleITK.WriteImage(img, file) #also in CT geometry if True: file1=getData.getLocalPath(r,xsetup,config.getPattern('centerNRRD',code=code,nclass=nclass,ir=ir,qaName='CT')) #method nearest perservse content, which should be id of the k-means group center u1=geometry.toSpace2(u,gSPECT,ct,gCT,method='nearest') img1=SimpleITK.GetImageFromArray(u1) img1.SetOrigin(gCT.origin) img1.SetSpacing(gCT.spacing) img1.SetDirection(numpy.ravel(gCT.direction)) SimpleITK.WriteImage(img1, file1) #write center map as numpy array qFile=getData.getLocalPath(r,xsetup,config.getPattern('centerMap',code=code,nclass=nclass,ir=ir,ic=i)) usave=numpy.zeros(kmeans.labels_.shape[0]+3) usave[0:3]=spect.shape[0:3] usave[3:]=kmeans.labels_ numpy.savetxt(qFile,usave,delimiter=',') def loadCenters(r,xsetup,ir=0): nclass=xsetup['nclass'][0] centers=numpy.array(0) for i in range(nclass): cFile=os.path.join(config.getLocalDir(r,xsetup),config.getCenter(r,xsetup,nclass,ir,i)) #row c=numpy.loadtxt(cFile,delimiter=',') if len(centers.shape)==0: centers=numpy.zeros((nclass,len(c))) centers[i,:]=c return centers def loadCenterMap(r,xsetup,ir=0): nclass=xsetup['nclass'][0] code=config.getCode(r,xsetup) qFile=getData.getLocalPath(r,xsetup,config.getPattern('centerMap',code=code,nclass=nclass,ir=ir)) usave=numpy.loadtxt(qFile,delimiter=',') shape=[int(x) for x in usave[0:3]] u=numpy.reshape(usave[3:],shape) return u def loadCenterMapNRRD(r,xsetup,ir=0): nclass=xsetup['nclass'][0] code=config.getCode(r,xsetup) md=['CT','SPECT'] files={x:getData.getLocalPath(r,xsetup,config.getPattern('centerNRRD',code=code,nclass=nclass,ir=ir,qaName=x)) for x in md} xd={x:SimpleITK.ReadImage(files[x]) for x in files} nd={x:SimpleITK.GetArrayFromImage(xd[x]) for x in xd} return nd['SPECT'],nd['CT'] def saveIVF(r,xsetup,ir=0,nfit=30,nbatch=30,qLambda=0): #fit IVF from centers in realization ir, perform nfit optimized fits where nbatch is used #to find best among nbatch trials (in total, nfit*nbatch fits will be made) #requires saveCenters to be run prior to execution nclass=xsetup['nclass'][0] code=config.getCode(r,xsetup) t,dt=loadTime(r,xsetup) centers=loadCenters(r,xsetup,ir) m,samples=fitData.fitIVFGlobal(t,centers,nfit=nfit,qLambda=qLambda) fm=m*numpy.ones(samples.shape[1]) fw=numpy.vstack((fm,samples)) f=getData.getLocalPath(r,xsetup,config.getPattern('fitIVF',code=code,nclass=nclass,ir=ir,qLambda=qLambda)) print(f'Saving to {f}') numpy.savetxt(f,fw,delimiter=',') def readIVF(r,xsetup,ir=0,qLambda=0): nclass=xsetup['nclass'][0] code=config.getCode(r,xsetup) f=getData.getLocalPath(r,xsetup,config.getPattern('fitIVF',code=code,nclass=nclass,ir=ir,qLambda=qLambda)) fw=numpy.loadtxt(f,delimiter=',') m=int(fw[0,0]) samples=fw[1:,:] return m,samples def saveSamples(r,xsetup,samples,m,name,iseg=0,ir=0,qLambda=0,qLambdaC=0): nclass=xsetup['nclass'][0] code=config.getCode(r,xsetup) fm=numpy.zeros(samples.shape[1]) n=numpy.min(numpy.array([len(m),samples.shape[0]])) for i in range(n): fm[i]=m[i]+1 fw=numpy.vstack((fm,samples)) pat=config.getPattern('fitCompartment',code=code,nclass=nclass,ir=ir,qaName=name,iseg=iseg,qLambda=qLambda,qLambdaC=qLambdaC) f=getData.getLocalPath(r,xsetup,pat) print(f'Saving samples to {f}') numpy.savetxt(f,fw,delimiter=',') def readSamples(r,xsetup,name,iseg=0,ir=0,qLambda=0,qLambdaC=0): m=[] nclass=xsetup['nclass'][0] code=config.getCode(r,xsetup) pat=config.getPattern('fitCompartment',code=code,nclass=nclass,ir=ir,qaName=name,iseg=iseg,qLambda=qLambda,qLambdaC=qLambdaC) f=getData.getLocalPath(r,xsetup,pat) print(f'Reading from {f}') fw=numpy.loadtxt(f,delimiter=',') samples=fw[1:,:] fm=samples[0,:] for i in range(fm.shape[0]): if fm[i]==0: break m.append(fm[i]-1) return m,samples def saveTAC(r,xsetup,tac,name,iseg=0,ir=0,qLambda=0,qLambdaC=0): nclass=xsetup['nclass'][0] code=config.getCode(r,xsetup) pat=config.getPattern('fitCompartment',code=code,nclass=nclass,ir=ir,qaName=name,iseg=iseg,qLambda=qLambda,qLambdaC=qLambdaC) f=getData.getLocalPath(r,xsetup,pat) print(f'Saving samples to {f}') numpy.savetxt(f,tac,delimiter=',') def readTAC(r,xsetup,name,iseg=0,ir=0,qLambda=0,qLambdaC=0): nclass=xsetup['nclass'][0] code=config.getCode(r,xsetup) pat=config.getPattern('fitCompartment',code=code,nclass=nclass,ir=ir,qaName=name,iseg=iseg,qLambda=qLambda,qLambdaC=qLambdaC) f=getData.getLocalPath(r,xsetup,pat) print(f'Reading from {f}') return numpy.loadtxt(f,delimiter=',')