소스 검색

Modified notebooks

Andrej 2 년 전
부모
커밋
034e09bc30
2개의 변경된 파일34개의 추가작업 그리고 413개의 파일을 삭제
  1. 19 340
      pythonScripts/clusterAnalysis.ipynb
  2. 15 73
      pythonScripts/segmentation.ipynb

+ 19 - 340
pythonScripts/clusterAnalysis.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [
     {
@@ -25,6 +25,7 @@
     "import re\n",
     "import config\n",
     "import getData\n",
+    "import analysis\n",
     "\n",
     "sys.path.append(os.path.join(os.path.expanduser('~'),'software','src','nixSuite','wrapper'))\n",
     "import nixWrapper\n",
@@ -32,336 +33,7 @@
     "import labkeyInterface\n",
     "import labkeyFileBrowser\n",
     "import labkeyDatabaseBrowser\n",
-    "\n",
-    "def calculateCenters(db,setup,nclass,nr):\n",
-    "    #path of scripts is ../scripts\n",
-    "    #one up\n",
-    "    baseDir=os.path.dirname(os.getcwd())\n",
-    "    #scriptPath\n",
-    "    scriptPath=os.path.join(baseDir,'scripts','generateCenters.sh')\n",
-    "    rows=getData.getPatients(db,setup)\n",
-    "    for r in rows:\n",
-    "        #download \n",
-    "        code=config.getCode(r,setup)\n",
-    "        for nc in nclass:\n",
-    "            print('centers: {}/{},{}'.format(code,nc,nr))\n",
-    "            print(subprocess.run([scriptPath,code,str(nc),str(nr)], \\\n",
-    "                             check=True, stdout=subprocess.PIPE).stdout)\n",
-    "            \n",
-    "def doAnalysis(db,setup,nclass,nr,mode):\n",
-    "    baseDir=os.path.dirname(os.getcwd())#one up\n",
-    "    #this is if IVF is inferred together with k1 fits\n",
-    "    #runScript='doAnalysis.sh'\n",
-    "    #this is if IVF is taken from previous cluster fits\n",
-    "    #runScript='doAnalysisIVF.sh'\n",
-    "    \n",
-    "    \n",
-    "    if mode=='global':\n",
-    "        runScript='doAnalysis.sh'\n",
-    "        analysisType=''\n",
-    "    if mode=='IVF':\n",
-    "        runScript='doAnalysisIVF.sh'\n",
-    "        analysisType='IVF_'\n",
-    "    \n",
-    "    try:\n",
-    "        print('Setting runScript to {}'.format(runScript))\n",
-    "    except NameError:\n",
-    "        print('Mode can be one of (global,IVF))')\n",
-    "        return\n",
-    "        \n",
-    "    #either doAnalysis or doAnalysisIVF.sh\n",
-    "    scriptPath=os.path.join(baseDir,'scripts',runScript)\n",
-    "\n",
-    "    rows=getData.getPatients(db,setup)\n",
-    "    for r in rows:\n",
-    "        dataPath=config.getLocalDir(r,setup)\n",
-    "        code=config.getCode(r,setup)\n",
-    "        for nc in nclass:\n",
-    "            for rId in numpy.arange(nr):\n",
-    "                print('{} [{} {}/{}]'.format(code,nc,rId+1,nr))\n",
-    "                #this is a duplicate of path generated in fitCenters.m\n",
-    "                fName=os.path.join(dataPath,'{}_{}_{}_{}fitParFinal.txt'.format(code,nc,rId+1,analysisType))\n",
-    "                if os.path.isfile(fName):\n",
-    "                    print('Skipping; {} available.'.format(fName))\n",
-    "                    continue\n",
-    "                subprocess.run([scriptPath,code,str(nc),str(rId+1)],\\\n",
-    "                         check=True, stdout=subprocess.PIPE)\n",
-    "\n",
-    "def doPixelAnalysis(db,setup,sigma2,mode='IVF'):\n",
-    "    \n",
-    "    baseDir=os.path.dirname(os.getcwd())#one up\n",
-    "                \n",
-    "    #in global mode, IVF parameters are inferred together with fits to classes\n",
-    "    #this is essentially repeat of the above, except that classes are taken as\n",
-    "    #time-response curves from pixels in the sigma2-defined neighborhood of\n",
-    "    #target pixels\n",
-    "    if mode=='global':\n",
-    "        runScript='doPixelAnalysis.sh'\n",
-    "        analysisType=''\n",
-    "    #in IVF mode, the parameters of input function are taken from the cluster fit \n",
-    "    #(doAnalysis above, with mode=general). The rest is the same as for global mode\n",
-    "    if mode=='IVF':\n",
-    "        runScript='doPixelIVFAnalysis.sh'\n",
-    "        analysisType='IVF_'\n",
-    "   \n",
-    "    #either doPixelAnalysis or doPixelAnalysisIVF.sh\n",
-    "    scriptPath=os.path.join(baseDir,'scripts',runScript)\n",
-    "\n",
-    "\n",
-    "    try:\n",
-    "        print('Setting runScript to {}'.format(runScript))\n",
-    "    except NameError:\n",
-    "        print('Mode can be one of (global,IVF))')\n",
-    "        return\n",
-    "    \n",
-    "    rows=getData.getPatients(db,setup)\n",
-    "    for r in rows:\n",
-    "        dataPath=config.getLocalDir(r,setup)\n",
-    "        code=config.getCode(r,setup)\n",
-    "        sName=os.path.join(dataPath,'{}_Segmentation.txt'.format(code))\n",
-    "        sFile=os.path.join(dataPath,sName)\n",
-    "        x=numpy.loadtxt(sFile)\n",
-    "        nc=x.shape[0]\n",
-    "        for s2 in sigma2:\n",
-    "            sigmaCode='{:.2f}'.format(s2)\n",
-    "            sigmaCode=re.sub('\\.','p',sigmaCode)\n",
-    "            fName=os.path.join(dataPath,'{}_{}_{}_Pixel{}_fitParFinal.txt'.format(code,nc,sigmaCode,analysisType))\n",
-    "            if os.path.isfile(fName):\n",
-    "                print('Skipping; {} available.'.format(fName))\n",
-    "                continue\n",
-    "            subprocess.run([scriptPath,code,str(s2)], \\\n",
-    "                             check=True, stdout=subprocess.PIPE)\n",
-    "            \n",
-    "            \n",
-    "def getIWeights(r,setup,nclass,realizationId,ic):\n",
-    "    locDir=config.getLocalDir(r,setup)\n",
-    "    code=config.getCode(r,setup)\n",
-    "    fname='{}_{}_{}_center{}_centerWeigth.nrrd'.\\\n",
-    "            format(code,nclass,realizationId+1,ic+1)\n",
-    "    uFile=os.path.join(locDir,fname)\n",
-    "            \n",
-    "    imU=SimpleITK.ReadImage(uFile)\n",
-    "    return SimpleITK.GetArrayFromImage(imU)\n",
-    "\n",
-    "def getGaussianWeight(nU,pt,sigma2,na):\n",
-    "    #gaussian weighted summation of surrounding pixels\n",
-    "    \n",
-    "    #find point closest to the target point\n",
-    "    idx=[int(x) for x in pt]\n",
-    "    #running offset        \n",
-    "    fidx=numpy.zeros(3)\n",
-    "    #half of the neighborhood\n",
-    "    na2=0.5*(na-1)\n",
-    "    \n",
-    "    fs=0\n",
-    "    fWeight=0\n",
-    "    for i in numpy.arange(na):\n",
-    "        fidx[0]=idx[0]+(i-na2)\n",
-    "        for j in numpy.arange(na):\n",
-    "            fidx[1]=idx[1]+(j-na2)\n",
-    "            for k in numpy.arange(na):\n",
-    "                fidx[2]=idx[2]+(k-na2)\n",
-    "                fidx=[int(x) for x in fidx]\n",
-    "                fd=numpy.array(fidx)-numpy.array(pt)\n",
-    "                dist2=numpy.sum(fd*fd)\n",
-    "                fw=numpy.exp(-0.5*dist2/sigma2);\n",
-    "                fs+=fw\n",
-    "                fWeight+=fw*nU[tuple(fidx)]\n",
-    "                #print('\\t{}/{}: {}/{:.2f} {:.2g} {:.3g} {:.2g}'.format(idx,fidx,numpy.sqrt(dist2),fw,nU[tuple(fidx)],fs,fWeight))\n",
-    "    fWeight/=fs\n",
-    "    return fWeight\n",
-    "            \n",
-    "\n",
-    "\n",
-    "#gets weights by class for a particular realization and sigma2 averaging\n",
-    "def getWeights(db,r,setup,nclass,realizationId,sigma2,na):\n",
-    "    #for w1, classes are in 0 to nclass-1 space\n",
-    "    #na is the size of the neighborhood\n",
-    "    idFilter=config.getIdFilter(r,setup)\n",
-    "    visitFilter=config.getVisitFilter(r,setup)\n",
-    "    code=config.getCode(r,setup)\n",
-    "    rows=getData.getSegmentation(db,setup,[idFilter,visitFilter])\n",
-    "    pts={r['regionId']:[float(x) for x in [r['x'],r['y'],r['z']]] for r in rows}\n",
-    "    \n",
-    "    w={region:numpy.zeros(nclass) for region in pts}\n",
-    "    na2=0.5*(na-1)\n",
-    "    \n",
-    "    for c in numpy.arange(nclass):\n",
-    "        nU=getIWeights(r,setup,nclass,realizationId,c)\n",
-    "        \n",
-    "        for region in w:\n",
-    "            #print(pts[region])\n",
-    "            #print('{} {}/{} {}'.format(code,c,nclass,region))\n",
-    "            #gaussian weighted summation of surrounding pixels\n",
-    "            w[region][c]=getGaussianWeight(nU,pts[region],sigma2,na)\n",
-    "            \n",
-    "            \n",
-    "    return w\n",
-    "\n",
-    "\n",
-    "#gets fitPar for a particular realization in [0..nr-1] range\n",
-    "def getFitPar(r,setup,nclass,realizationId,mode=''):\n",
-    "    #fitGoodnes A tau alpha delay [k1 BVF k2 delay]xNcenters\n",
-    "    allowedModes=['','IVF','Pixel','PixelIVF']\n",
-    "    if mode not in allowedModes:\n",
-    "        print('Mode should be one of {}'.format(allowedModes))\n",
-    "        return []\n",
-    "    \n",
-    "    if mode=='PixelIVF':\n",
-    "        #4th parameter is sigma, not realizationId\n",
-    "        rCode='{:.2f}'.format(realizationId)\n",
-    "        rCode=re.sub('\\.','p',rCode)\n",
-    "    else:\n",
-    "        #add one to match matlab 1..N array indexing\n",
-    "        rCode='{}'.format(realizationId+1)\n",
-    "    \n",
-    "    if len(mode)>0:\n",
-    "        mode=mode+'_'\n",
-    "    \n",
-    "    code=config.getCode(r,setup)\n",
-    "    fname='{}_{}_{}_{}fitParFinal.txt'.format(code,nclass,rCode,mode)\n",
-    "    locDir=config.getLocalDir(r,setup)\n",
-    "    uFile=os.path.join(locDir,fname)\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(w,fitPar):\n",
-    "    #calculateK1 for region weights\n",
-    "    #return the k1 that belongs to the \n",
-    "    #maximum class in region (M) and \n",
-    "    #a weighted average (W)\n",
-    "    k1={region:{'M':0,'W':0} for region in w}\n",
-    "    for region in w:\n",
-    "        cmax=numpy.argmax(w[region])\n",
-    "        k1[region]['M']=getK1(fitPar,cmax)\n",
-    "        fs=0\n",
-    "        for c in numpy.arange(len(w[region])):\n",
-    "            fs+=w[region][c]*getK1(fitPar,c)\n",
-    "        k1[region]['W']=fs\n",
-    "    return k1\n",
-    "\n",
-    "def getPatientValuesByNclass(db,r,setup,nclass,nrealizations,sigma2,na):\n",
-    "    #summary script\n",
-    "    #db is for database; needs segmentations\n",
-    "    #r,setup identify patient\n",
-    "    #nclass and nrealizations select strategy\n",
-    "    #sigma2 is for combining output from adjacent pixels\n",
-    "    #na is neighborhood size where smoothing/combination is done\n",
-    "    k1={}\n",
-    "    for rId in numpy.arange(nrealizations):\n",
-    "        w=getWeights(db,r,setup,nclass,rId,sigma2,na)\n",
-    "        fitPar=getFitPar(r,setup,nclass,rId,'IVF')\n",
-    "        qk1=calculateK1(w,fitPar)\n",
-    "        for region in w:\n",
-    "            for type in qk1[region]:\n",
-    "                try:\n",
-    "                    k1[region][type].append(qk1[region][type])\n",
-    "                except KeyError:\n",
-    "                    k1={region:{type:[] for type in qk1[region]} for region in w}\n",
-    "                    print(type)\n",
-    "                    k1[region][type].append(qk1[region][type])\n",
-    "        print('[{}] {}/{}'.format(nclass,rId+1,nrealizations))\n",
-    "    return k1   \n",
-    "\n",
-    "def getSummaryPatientValuesByNclass(db,r,setup,nclass,nrealizations,sigma2,na):\n",
-    "    #summary script, same arguments as above\n",
-    "    #also return deviation over realization\n",
-    "    k1=getPatientValuesByNclass(db,r,setup,nclass,nrealizations,sigma2,na)\n",
-    "    avgType=['M','W']\n",
-    "    summaryK1={type:{region:{\n",
-    "        'mean':numpy.mean(k1[region][type]), \n",
-    "        'std':numpy.std(k1[region][type]), \n",
-    "        'median':numpy.median(k1[region][type])} for region in k1}\n",
-    "               for type in avgType}\n",
-    "    \n",
-    "    return summaryK1\n",
-    "\n",
-    "def fullSummary(db,setup,classes,nr,sigma2,na):\n",
-    "    rows=getData.getPatients(db,setup)\n",
-    "    return \\\n",
-    "        {config.getCode(r,setup):\\\n",
-    "         {c:getSummaryPatientValuesByNclass(db,r,setup,c,nr,sigma2,na) for c in classes} for r in rows}\n",
-    "          \n",
-    "def storeSummary(db,setup,summary,sigma2,na):\n",
-    "    #dsM=db.selectRows(project,'study','Imaging',[])\n",
-    "    for rCode in summary:\n",
-    "        r=config.decode(rCode,setup)\n",
-    "        idFilter=config.getIdFilter(r,setup)\n",
-    "        visitFilter=config.getVisitFilter(r,setup)\n",
-    "        for c in summary[rCode]:\n",
-    "            cFilter={'variable':'nclass','value':str(c),'oper':'eq'}\n",
-    "            for t in summary[rCode][c]:\n",
-    "                tFilter={'variable':'option','value':t,'oper':'eq'}\n",
-    "                for rId in summary[rCode][c][t]:\n",
-    "                    rFilter={'variable':'regionId','value':str(rId),'oper':'eq'}\n",
-    "                    rows=getData.getSummary(db,setup,[idFilter,visitFilter,cFilter,tFilter,rFilter])\n",
-    "                    if len(rows)>0:\n",
-    "                        qrow=rows[0]\n",
-    "                        mode='update'\n",
-    "                    else:\n",
-    "                        qrow={qr:r[qr] for qr in r}\n",
-    "                        qrow['nclass']=c\n",
-    "                        qrow['option']=t\n",
-    "                        qrow['regionId']=rId\n",
-    "                        seqNum=config.getTargetSeqNum(r,setup)\n",
-    "                        qrow['SequenceNum']=100+seqNum+c+0.001*rId\n",
-    "                        if t=='M':\n",
-    "                            qrow['SequenceNum']+=0.0001\n",
-    "                        mode='insert'\n",
-    "                    for v in summary[rCode][c][t][rId]:\n",
-    "                        qrow[v]=summary[rCode][c][t][rId][v]\n",
-    "                    qrow['sigma2']=sigma2\n",
-    "                    qrow['na']=na\n",
-    "                    getData.updateSummary(db,setup,mode,[qrow])\n",
-    "                    \n",
-    "def summaryPixelIVF(db,setup,sigma2):\n",
-    "    #for second type of analysis (pixel based regions)\n",
-    "    rows=getData.getPatients(db,setup)\n",
-    "    return \\\n",
-    "        {config.getCode(r,setup):\\\n",
-    "         {s2:getPixelIVF(db,r,setup,s2) for s2 in sigma2} for r in rows}\n",
-    "\n",
-    "    \n",
-    "def storeIVF(db,setup,summary):\n",
-    "    for rCode in summary:\n",
-    "        r=config.decode(rCode,setup)\n",
-    "        idFilter=config.getIdFilter(r,setup)\n",
-    "        visitFilter=config.getVisitFilter(r,setup)\n",
-    "        for s2 in summary[rCode]:\n",
-    "            sigmaFilter={'variable':'sigma2','value':str(s2),'oper':'eq'}\n",
-    "            nr=len(summary[rCode][s2])\n",
-    "            for rId in summary[rCode][s2]:\n",
-    "                rFilter={'variable':'regionId','value':str(rId),'oper':'eq'}\n",
-    "                typeFilter={'variable':'option','value':'D','oper':'eq'}\n",
-    "                rows=getData.getSummary(db,setup,[idFilter,visitFilter,sigmaFilter,rFilter,typeFilter])\n",
-    "                if len(rows)>0:\n",
-    "                    qrow=rows[0]\n",
-    "                    mode='update'\n",
-    "                else:\n",
-    "                    qrow={qr:r[qr] for qr in r}\n",
-    "                    qrow['sigma2']=s2\n",
-    "                    qrow['regionId']=rId\n",
-    "                    seqNum=config.getTargetSeqNum(r,setup)\n",
-    "                    qrow['SequenceNum']=140+seqNum+0.01*rId+0.001*s2\n",
-    "                    mode='insert'\n",
-    "                qrow['mean']=summary[rCode][s2][rId]\n",
-    "                qrow['na']=7\n",
-    "                qrow['nclass']=nr\n",
-    "                qrow['option']='D'\n",
-    "                getData.updateSummary(db,setup,mode,[qrow])\n",
-    "                \n",
-    "def getPixelIVF(db,r,setup,sigma2):\n",
-    "    idFilter=config.getIdFilter(r,setup)\n",
-    "    visitFilter=config.getVisitFilter(r,setup)\n",
-    "    rows=getData.getSegmentation(db,setup,[idFilter,visitFilter])\n",
-    "    nclassIVF=len(rows)\n",
-    "    fitPar=getFitPar(r,setup,nclassIVF,sigma2,'PixelIVF')\n",
-    "    k1={r['regionId']:getK1(fitPar,r['regionId']) for r in rows}\n",
-    "    return k1"
+    "\n"
    ]
   },
   {
@@ -412,33 +84,40 @@
     "#takes multiple hours for nr=20, nclass=[10,20,30]\n",
     "#runs matlab - maybe the first time licensing will be checked and error will be returned\n",
     "#just run once\n",
-    "#calculateCenters(db,setup,nclass,nr)\n",
+    "#analysis.calculateCenters(db,setup,nclass,nr)\n",
     "\n",
     "#3.) now do analysis - global gets the average IVF fit (see point one of the slides for both branches of analysis)\n",
     "#this again takes multiple hours for four images\n",
     "#from what I see this could be shortend (since only global parameters are sought)\n",
-    "#doAnalysis(db,setup,nclass,nr,'global')\n",
+    "#analysis.doAnalysis(db,setup,nclass,nr,'global')\n",
     "\n",
     "#4.) do cluster analysis (bullets two and three under ML analysis)\n",
-    "#doAnalysis(db,setup,nclass,nr,'IVF')\n",
+    "#analysis.doAnalysis(db,setup,nclass,nr,'IVF')\n",
     "\n",
     "#5.) sort out the segmentation points -> segmentation.ipynb and generate Segmentation.txt files read by pixelAnalysis\n",
     "\n",
     "#6.) run doPixelAnalysisIVF to calculate parameters from locations\n",
     "sigma2=[0.1,1,4]\n",
-    "#doPixelAnalysis(db,setup,sigma2)\n",
+    "#analysis.doPixelAnalysis(db,setup,sigma2)\n",
     "\n",
     "#7.) store summary to database\n",
     "na=3\n",
     "s2=0.1\n",
-    "#summary=fullSummary(db,setup,nclass,nr,s2,na)\n",
-    "#storeSummary(db,setup,summary,s2,na)\n",
+    "#summary=analysis.fullSummary(db,setup,nclass,nr,s2,na)\n",
+    "#analysis.storeSummary(db,setup,summary,s2,na)\n",
     "\n",
-    "#summaryIVF=summaryPixelIVF(db,setup,sigma2)\n",
-    "#storeIVF(db,setup,summaryIVF)\n",
+    "summaryIVF=analysis.summaryPixelIVF(db,setup,sigma2)\n",
+    "analysis.storeIVF(db,setup,summaryIVF)\n",
     "\n",
-    "#8.) Plot regions through chart wizard on labkey\n"
+    "#8.) Plot regions through chart wizard on labkey, see also Summary R Report for error-plots on multi-class analysis\n"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {

파일 크기가 너무 크기때문에 변경 상태를 표시하지 않습니다.
+ 15 - 73
pythonScripts/segmentation.ipynb


이 변경점에서 너무 많은 파일들이 변경되어 몇몇 파일들은 표시되지 않았습니다.