{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#functions used in final analysis\n", "def getWeights(locDir,p,nclass,realizationId):\n", " #for w1, classes are in 0 to nclass-1 space\n", " \n", " segmFile=os.path.join(locDir,p,p+'_Segm.nrrd')\n", " imS=SimpleITK.ReadImage(segmFile)\n", " nS=SimpleITK.GetArrayFromImage(imS)\n", " \n", " unit=numpy.ones(nS.shape)\n", " w={region:numpy.zeros(nclass) for region in numpy.arange(1,10)}\n", " w1={region:{} for region in w}\n", " \n", " sUnit={region:numpy.sum(unit[nS==region]) for region in w}\n", " \n", " for c in numpy.arange(nclass):\n", " uFile=os.path.join(locDir,p,'{}_{}_{}_center{}_centerWeigth.nrrd'.format(p,nclass,realizationId,c+1))\n", " imU=SimpleITK.ReadImage(uFile)\n", " nU=SimpleITK.GetArrayFromImage(imU)\n", " \n", " for region in w:\n", " v=nU[(nS==region) ]\n", " #print('{}: {}'.format(region,v.shape))\n", " w[region][c]=numpy.sum(v)/sUnit[region]\n", " w1[region][c]=v\n", " #print('[{}/{}]: {}'.format(region,j,w[region,j]))\n", " return w,w1\n", "\n", "def getFuzziness(locDir,p,nclass,realizationId):\n", " \n", " segmFile=os.path.join(locDir,p,p+'_Segm.nrrd')\n", " imS=SimpleITK.ReadImage(segmFile)\n", " nS=SimpleITK.GetArrayFromImage(imS)\n", " \n", " w={region:0 for region in numpy.arange(1,10)}\n", " nU={}\n", " for c in numpy.arange(nclass):\n", " uFile=os.path.join(locDir,p,'{}_{}_{}_center{}_centerWeigth.nrrd'.format(p,nclass,realizationId,c+1))\n", " imU=SimpleITK.ReadImage(uFile)\n", " fy=SimpleITK.GetArrayFromImage(imU)\n", " \n", " for region in w:\n", " fy1=fy[(nS==region)]\n", " try:\n", " nU[region][:,c]=fy1\n", " except KeyError:\n", " nU[region]=numpy.zeros((len(fy1),nclass))\n", " nU[region][:,c]=fy1\n", " for u in nU:\n", " w[u]=numpy.mean(nU[u].max(1))\n", " return w\n", " \n", "\n", "def getFitPar(locDir,p,nclass,realizationId):\n", " #fitGoodnes A tau alpha delay [k1 BVF k2 delay]xNcenters\n", " uFile=os.path.join(locDir,p,'{}_{}_{}_fitParFinal.txt'.format(p,nclass,realizationId))\n", " return numpy.genfromtxt(uFile,delimiter='\\t')\n", "\n", "def getK1(fitPar,iclass):\n", " #fitGoodnes A tau alpha delay [k1 BVF k2 delay]xNcenters\n", " #0 to nclass-1 space\n", " return fitPar[4*iclass+5]\n", "\n", "def calculateK1(w1,fitPar,range=1.0, nbins=50):\n", " #returns k1 by pixel in w1\n", " histM={}\n", " histW={}\n", " classWeightM={}\n", " classWeightW={}\n", " fuzzy={}\n", " std={'W':{},'M':{}}\n", " median={'W':{},'M':{}}\n", " f90p={'W':{},'M':{}}\n", " f95p={'W':{},'M':{}}\n", " mean={'W':{},'M':{}}\n", " for region in w1:\n", " vTempl=w1[region][list(w1[region].keys())[0]]#should by a numpy array with shape (npx,)\n", " v={v:{} for v in numpy.arange(len(vTempl))}\n", " for c in w1[region]:\n", " for i in v:\n", " v[i][c]=w1[region][c][i]\n", " fuzzy[region]=0\n", " k1W=numpy.zeros(vTempl.shape)\n", " k1M=numpy.zeros(vTempl.shape)\n", " cW=numpy.zeros(len(w1[region]))#number of classes\n", " cM=numpy.zeros(len(w1[region]))#number of classes\n", " for i in v:\n", " maxClass=max(v[i],key=v[i].get)\n", " maxValue=v[i][maxClass]\n", " fuzzy[region]+=maxValue\n", " k1W[i]=0\n", " for c in v[i]:\n", " k1W[i]+=v[i][c]*getK1(fitPar,c)\n", " cW[c]+=v[i][c]\n", " k1M[i]=getK1(fitPar,maxClass)\n", " cM[maxClass]+=1\n", " #print('[{},{}:{}] {},{}'.format(i,max(v[i],key=v[i].get),max(v[i].values()),k1W[i],k1M[i]))\n", " histM[region]=numpy.histogram(60*k1M,bins=nbins,range=(0,range))\n", " histW[region]=numpy.histogram(60*k1W,bins=nbins,range=(0,range))\n", " mean['M'][region]=numpy.mean(k1M)\n", " mean['W'][region]=numpy.mean(k1W)\n", " std['M'][region]=numpy.std(k1M)\n", " std['W'][region]=numpy.std(k1W)\n", " median['M'][region]=numpy.median(k1M)\n", " median['W'][region]=numpy.median(k1W)\n", " f90p['M'][region]=numpy.percentile(k1M,90)\n", " f90p['W'][region]=numpy.percentile(k1W,90)\n", " f95p['M'][region]=numpy.percentile(k1M,95)\n", " f95p['W'][region]=numpy.percentile(k1W,95)\n", " \n", " numpy.set_printoptions(precision=3)\n", " cM/=len(v)\n", " cW/=len(v)\n", " classWeightM[region]=cM\n", " classWeightW[region]=cW\n", " for c1,c2 in zip(cM,cW):\n", " pass\n", " #print('{}/{:.2f}'.format(c1,c2))\n", " fuzzy[region]/=len(v)\n", " \n", " return fuzzy,histM,histW,classWeightM,classWeightW,mean,std,median,f90p,f95p\n", "\n", "def makePlot(h):\n", " width = numpy.diff(h[1])\n", " center = (h[1][:-1] + h[1][1:]) / 2\n", " y=h[0]\n", " return center,y,width\n", "\n", "def median(x,y):\n", " N=numpy.sum(y)\n", " s=0\n", " i=0\n", " while s<0.5*N:\n", " s+=y[i]\n", " i+=1\n", " return x[i]\n", "\n", "def mean(x,y):\n", " return sum(x*y)/sum(y)\n", "\n", "#updates k1 values in ImagingData tables by averaging over realizations\n", "#used calculateK1, getWeights and getFitPar\n", "\n", "def updatePatient(db,fb,project, locDir,patient, nclass):\n", " idFilter={'variable':'PatientId','value':patient,'oper':'eq'}\n", " ds=db.selectRows(project,'study','Imaging',[idFilter])\n", " rows=ds['rows']\n", " #rows=[ds['rows'][1]]\n", " rId=1\n", " hMy={}\n", " hMx={}\n", " hMw={}\n", "\n", " hMean={}\n", " hStd={}\n", " hMedian={}\n", " h90p={}\n", " h95p={}\n", "\n", " selectFlag='M'\n", "\n", " for r in rows:\n", " pId=r['PatientId']\n", " if pId.find('8')==0:\n", " continue\n", " state='MIR'\n", " if r['SequenceNum']>1:\n", " state='OBR'\n", " p=r['aliasID']\n", "\n", " for rId in numpy.arange(1,21):\n", " w,w1=getWeights(locDir,p,nclass,rId)\n", " fitPar=getFitPar(locDir,p,nclass,rId)\n", " fz,hM,hW,cM,cW,mean,qstd,qmedian,f90p,f95p=\\\n", " calculateK1(w1,fitPar,range=0.8)\n", " #x,y,w = makePlot(hW[7])\n", " #matplotlib.pyplot.bar(x,y,width=w)\n", " #matplotlib.pyplot.text(0.3,20,'{} {}'.format(state,rId))\n", " #matplotlib.pyplot.show()\n", " for r in w1:\n", " if selectFlag=='M':\n", " x,y,w = makePlot(hM[r])\n", " else:\n", " x,y,w = makePlot(hW[r])\n", " try:\n", " fy=hMy[r][state]\n", " hMy[r][state]=fy+y\n", " except KeyError:\n", " try:\n", " hMy[r][state]=y\n", " hMx[r][state]=x\n", " hMw[r][state]=w\n", " except KeyError:\n", " hMy[r]={}\n", " hMx[r]={}\n", " hMw[r]={}\n", " hMy[r][state]=y\n", " hMx[r][state]=x\n", " hMw[r][state]=w\n", " try:\n", " hMean[r][state].append(mean[selectFlag][r])\n", " hMedian[r][state].append(qmedian[selectFlag][r])\n", " hStd[r][state].append(qstd[selectFlag][r])\n", " h90p[r][state].append(f90p[selectFlag][r])\n", " h95p[r][state].append(f95p[selectFlag][r])\n", " except KeyError:\n", " try:\n", " hMean[r][state]=[]\n", " hMedian[r][state]=[]\n", " hStd[r][state]=[]\n", " h90p[r][state]=[]\n", " h95p[r][state]=[]\n", " except KeyError:\n", " hMean[r]={}\n", " hMedian[r]={}\n", " hStd[r]={}\n", " h90p[r]={}\n", " h95p[r]={}\n", "\n", "\n", " #generate table\n", " nclassFilter={'variable':'nclass','value':str(nclass),'oper':'eq'}\n", " ds=db.selectRows(project,'study','ImagingData',[idFilter,nclassFilter])\n", "\n", " for qrow in ds['rows']:\n", " r=qrow['regionId']\n", "\n", " ymax=-1e30\n", " s={}\n", " for state in hMy[r]:\n", " y=hMy[r][state]\n", " x=hMx[r][state]\n", " s[state]=sum(x*y)/numpy.sum(y)\n", " s[state]=median(x,y)\n", " if max(y)>ymax:\n", " ymax=max(y)\n", " matplotlib.pyplot.bar(hMx[r][state],hMy[r][state],width=hMw[r][state])\n", "\n", " #fill variables\n", " qrow['k1'+state+'Mean']=numpy.mean(hMean[r][state])\n", " qrow['k1'+state+'MeanStd']=numpy.std(hMean[r][state])\n", " qrow['k1'+state+'Median']=numpy.mean(hMedian[r][state])\n", " qrow['k1'+state+'MedianStd']=numpy.std(hMedian[r][state])\n", " qrow['k1'+state+'Std']=numpy.mean(hStd[r][state])\n", " qrow['k1'+state+'StdStd']=numpy.std(hStd[r][state])\n", " qrow['k1'+state+'90p']=numpy.mean(h90p[r][state])\n", " qrow['k1'+state+'90pStd']=numpy.std(h90p[r][state])\n", " qrow['k1'+state+'95p']=numpy.mean(h90p[r][state])\n", " qrow['k1'+state+'95pStd']=numpy.std(h90p[r][state])\n", "\n", " p=qrow['PatientId']\n", " oDir=os.path.join(locDir,p)\n", " oDirRemote=fb.buildPathURL(project,[p]);\n", " qrow['k1RatioMean']=qrow['k1MIRMean']/qrow['k1OBRMean']\n", " qrow['k1RatioMedian']=qrow['k1MIRMedian']/qrow['k1OBRMedian']\n", " qrow['k1Ratio90p']=qrow['k1MIR90p']/qrow['k1OBR90p']\n", " qrow['k1Ratio95p']=qrow['k1MIR95p']/qrow['k1OBR95p']\n", " db.modifyRows('update',project,'study','ImagingData',[qrow])\n", " matplotlib.pyplot.text(0.4,0.95*ymax,'{} region:{}'.format(qrow['PatientId'], r))\n", " matplotlib.pyplot.text(0.4,0.85*ymax,'k1={:.2f}/min k2={:.2f}/min'.format(s['MIR'], s['OBR']))\n", " matplotlib.pyplot.text(0.4,0.75*ymax,'R={:.2f}'.format(s['OBR']/s['MIR']))\n", " fname='{}_{}_{}_k1Distribution.png'.format(p,nclass,r)\n", " matplotlib.pyplot.savefig(os.path.join(oDir,fname))\n", " fb.writeFileToFile(os.path.join(oDir,fname),oDirRemote+'/'+fname)\n", " matplotlib.pyplot.show()\n", " \n", " \n", "def getReadings(db,fb,project, locDir,patient, nclass, state, region):\n", " readings=[]\n", " idFilter={'variable':'PatientId','value':patient,'oper':'eq'}\n", " if state=='MIR':\n", " seqNum=1\n", " else:\n", " seqNum=2\n", " seqFilter={'variable':'SequenceNum','value':str(seqNum),'oper':'eq'}\n", " ds=db.selectRows(project,'study','Imaging',[idFilter,seqFilter])\n", " r=ds['rows'][0]\n", " p=r['aliasID']\n", " weightOption='M'\n", " for realizationId in numpy.arange(1,21):\n", " w,w1=getWeights(locDir,p,nclass,realizationId)\n", " fitPar=getFitPar(locDir,p,nclass,realizationId)\n", " fz,hM,hW,cM,cW,mean,qstd,qmedian,f90p,f95p=calculateK1(w1,fitPar,range=0.8)\n", " readings.append(qmedian[weightOption][region])\n", " return readings\n", "\n", "import matplotlib.pyplot\n", "import matplotlib.image as mpimg\n", "import json\n", "\n", "#extract patient values for a particular realization,\n", "#return region averages and histograms of k1 in region\n", "#based on pixel selection criteria, which is either 'M' \n", "#for maximum probability or 'W' for probability weighted contribution\n", "#of cluster in each pixel (see calculate K1)\n", "\n", "#fill ImageDataRealization\n", "def getPatient(db,fb,project, locDir,patient, nclass, realizationId,\\\n", " state, region):\n", " avgMode='M'\n", " idFilter={'variable':'PatientId','value':patient,'oper':'eq'}\n", " if state=='MIR':\n", " seqNum=1\n", " else:\n", " seqNum=2\n", " seqFilter={'variable':'SequenceNum','value':str(seqNum),'oper':'eq'}\n", " ds=db.selectRows(project,'study','Imaging',[idFilter,seqFilter])\n", " r=ds['rows'][0]\n", " \n", " p=r['aliasID']\n", " oDirRemote=fb.buildPathURL(project,[p]);\n", " \n", " w,w1=getWeights(locDir,p,nclass,realizationId)\n", " fitPar=getFitPar(locDir,p,nclass,realizationId)\n", " fz,hM,hW,cM,cW,mean,qstd,qmedian,f90p,f95p=calculateK1(w1,fitPar,range=0.8)\n", " #plot most common responses\n", " cmSort=numpy.flip(numpy.argsort(cM[region]))\n", " \n", " #last argument\n", " for i in numpy.arange(3):\n", " idx=cmSort[i]\n", " #print('{}: {}'.format(idx,cM[region][idx]))\n", " fname='{}_{}_{}_centers{}.png'.format(p,nclass,realizationId,idx+1)\n", " fp=os.path.join(locDir,p,fname)\n", " img = mpimg.imread(fp)\n", " #imgplot = matplotlib.pyplot.imshow(img)\n", " #matplotlib.pyplot.show()\n", " rf=oDirRemote+'/'+fname\n", " if not fb.entryExists(rf):\n", " fb.writeFileToFile(fp,rf) \n", " \n", " realizationFilter={'variable':'realizationId','value':str(realizationId),'oper':'eq'}\n", " stateFilter={'variable':'state','value':state,'oper':'eq'}\n", " nclassFilter={'variable':'nclass','value':str(nclass),'oper':'eq'}\n", " regionFilter={'variable':'regionId','value':str(region),'oper':'eq'}\n", " \n", " #fill database\n", " dsQ=db.selectRows(project,'study','ImagingDataRealization',[idFilter,stateFilter,nclassFilter,regionFilter,realizationFilter])\n", " if len(dsQ['rows'])>0:\n", " qrow=dsQ['rows'][0]\n", " mode='update'\n", " else:\n", " qrow={}\n", " qrow['PatientId']=patient\n", " qrow['SequenceNum']=nclass+0.01*region+0.0001*realizationId\n", " if state!='MIR':\n", " qrow['SequenceNum']+=1\n", " qrow['realizationId']=int(realizationId)\n", " qrow['state']=state\n", " qrow['nclass']=int(nclass)\n", " qrow['regionId']=int(region)\n", " mode='insert'\n", " \n", " qrow['k1Mean']=mean[avgMode][region]\n", " qrow['k1Median']=qmedian[avgMode][region]\n", " qrow['k190p']=f90p[avgMode][region]\n", " qrow['k195p']=f95p[avgMode][region]\n", " qrow['k1Std']=qstd[avgMode][region]\n", " qrow['centers']=';'.join([str(cmSort[i]) for i in numpy.arange(3)])\n", " qrow['cW']=';'.join(['{:.2f}'.format(cM[region][cmSort[i]]) for i in numpy.arange(3)])\n", " db.modifyRows(mode,project,'study','ImagingDataRealization',[qrow])\n", " \n", " \n", " fname='{}_{}_{}_histogram{}.png'.format(p,nclass,realizationId,region)\n", " fp=os.path.join(locDir,p,fname)\n", " x,y,w=makePlot(hM[region])\n", " if state=='MIR':\n", " color='tab:blue'\n", " else:\n", " color='tab:orange'\n", " imgplot = matplotlib.pyplot.bar(x,y,width=w,color=color)\n", " matplotlib.pyplot.savefig(fp)\n", " rf=oDirRemote+'/'+fname\n", " if not fb.entryExists(rf):\n", " fb.writeFileToFile(fp,rf) \n", " \n", " #print('Saved to {}'.format(oDirRemote+'/'+fname))\n", " #write file to file\n", " #matplotlib.pyplot.show()\n", " \n", " \n", " jsonHistogram={'x':x.tolist(),'y':y.tolist(),'w':w.tolist()}\n", " fname='{}_{}_{}_histogram{}.json'.format(p,nclass,realizationId,region)\n", " fp=os.path.join(locDir,p,fname)\n", " with open (fp,'w') as f:\n", " json.dump(jsonHistogram,f)\n", " #write to server\n", " rf=oDirRemote+'/'+fname\n", " if not fb.entryExists(rf):\n", " fb.writeFileToFile(fp,rf)\n", " #fitPar?\n", " fname='{}_{}_{}_fitParFinal.txt'.format(p,nclass,realizationId)\n", " fpRemote=oDirRemote+'/'+fname\n", " if not fb.entryExists(fpRemote):\n", " fp=os.path.join(locDir,p,fname)\n", " if os.path.isfile(fp):\n", " fb.writeFileToFile(fp,fpRemote)\n", " \n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "nclass=30\n", "ds=db.selectRows(project,'study','Imaging',[])\n", "patients=list(set([row['PatientId'] for row in ds['rows']]))\n", "idFilter={'variable':'PatientId','value':patients[4],'oper':'eq'}\n", "ds=db.selectRows(project,'study','Imaging',[idFilter])\n", "rows=ds['rows']\n", "#rows=[ds['rows'][1]]\n", "rId=1\n", "k1Mean={}\n", "k1Std={}\n", "k1Median={}\n", "k190p={}\n", "k195p={}\n", "\n", "for r in rows:\n", " pId=r['PatientId']\n", " if pId.find('8')==0:\n", " continue\n", " state='MIR'\n", " if r['SequenceNum']>1:\n", " state='OBR'\n", " p=r['aliasID']\n", " \n", " for rId in numpy.arange(1,21):\n", " w,w1=getWeights(locDir,p,nclass,rId)\n", " fitPar=getFitPar(locDir,p,nclass,rId)\n", " fz,hM,hW,cM,cW=calculateK1(w1,fitPar,range=0.8)\n", " #x,y,w = makePlot(hW[7])\n", " #matplotlib.pyplot.bar(x,y,width=w)\n", " #matplotlib.pyplot.text(0.3,20,'{} {}'.format(state,rId))\n", " #matplotlib.pyplot.show()\n", " for r in w1:\n", " x,y,w = makePlot(hW[r])\n", " try:\n", " k1Mean\n", " fy=hMy[r][state]\n", " hMy[r][state]=fy+y\n", " except KeyError:\n", " try:\n", " hMy[r][state]=y\n", " hMx[r][state]=x\n", " hMw[r][state]=w\n", " except KeyError:\n", " hMy[r]={}\n", " hMx[r]={}\n", " hMw[r]={}\n", " hMy[r][state]=y\n", " hMx[r][state]=x\n", " hMw[r][state]=w\n", " \n", "#generate table\n", "\n", " \n", "for r in hMy:\n", " ymax=-1e30\n", " s={}\n", " for state in hMy[r]:\n", " y=hMy[r][state]\n", " x=hMx[r][state]\n", " s[state]=sum(x*y)/numpy.sum(y)\n", " s[state]=median(x,y)\n", " if max(y)>ymax:\n", " ymax=max(y)\n", " matplotlib.pyplot.bar(hMx[r][state],hMy[r][state],width=hMw[r][state])\n", " \n", " matplotlib.pyplot.text(0.5,0.8*ymax,'{} {:.2f} {:.2f} {:.2f}'.format(r, s['MIR'], s['OBR'], s['OBR']/s['MIR']))\n", " matplotlib.pyplot.show()\n", " \n", " \n", "readings=getReadings(db,fb,project,locDir,'3ZF',30,'OBR',4)\n", "print('{}/{} {}'.format(numpy.mean(readings),numpy.median(readings),readings))\n", "\n", "patients=['2SB','3ZF','5MI','7TM','10KF','11ZM']\n", "states=['MIR','OBR']\n", "classes=[10,20,30]\n", "regions=numpy.arange(1,10)\n", "realizations=numpy.arange(1,21)\n", "for p in patients:\n", " for s in states:\n", " for c in classes:\n", " for realizationId in realizations:\n", " for r in regions:\n", " print('{} {} {} {}/20 {}/9'.format(p,s,c,realizationId,r))\n", " getPatient(db,fb,project,locDir,p,c,realizationId,s,r)\n", " \n", "#update k1 values in ImagingData from getK1 \n", "prows=[row for row in ds['rows'] if row['PatientId']==pId]\n", " print(len(prows))\n", " \n", " for prow in prows:\n", " idFilter={'variable':'PatientId','value':pId,'oper':'eq'}\n", " nclassFilter={'variable':'nclass','value':str(nclass),'oper':'eq'}\n", " dsQ=db.selectRows(project,'study','ImagingData',[idFilter,nclassFilter])\n", " \n", " #skip those already done\n", " #if len(dsQ['rows'])>0:\n", " # continue\n", "\n", " p=prow['aliasID']\n", " visitName='MIR'\n", " if prow['SequenceNum']>1:\n", " visitName='OBR'\n", " \n", " #here the sampling happens\n", " fq=[getK1(p,nclass,visitName,i) for i in numpy.arange(0,nsample)]\n", " qMerge=combine(fq)\n", " \n", " #this is fine also if previous step is omitted\n", " if len(dsQ['rows'])==0:\n", " regionIds=numpy.arange(1,10)\n", " else:\n", " regionIds=[int(r['regionId']) for r in dsQ['rows']]\n", "\n", " \n", " for region in regionIds:\n", " qrow=[r for r in dsQ['rows'] if r['regionId']==region]\n", " \n", " mode='update'\n", " if len(qrow):\n", " qrow=qrow[0]\n", " else:\n", " qrow={}\n", " qrow['PatientId']=prow['PatientId']\n", " qrow['SequenceNum']=str(nclass+0.01*region)\n", " qrow['regionId']=str(region)\n", " mode='insert'\n", " \n", " qrow['nclass']=str(nclass)\n", " vars=['k1','v19']\n", " avgs=['Mean','Std','Median','90p','95p']\n", " for var in vars:\n", " for avg in avgs:\n", " qrow[var+visitName+avg]=str(numpy.mean(qMerge[var+visitName+avg+'_'+str(region)]))\n", " qrow[var+visitName+avg+'Std']=str(numpy.std(qMerge[var+visitName+avg+'_'+str(region)]))\n", " db.modifyRows(mode,project,'study','ImagingData',[qrow])\n", " if i==0:\n", " break\n", " if i==0:" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.10" } }, "nbformat": 4, "nbformat_minor": 4 }