Difference between revisions of "User talk:Haehn"

From NAMIC Wiki
Jump to: navigation, search
 
(26 intermediate revisions by the same user not shown)
Line 1: Line 1:
 
==my notes..==
 
==my notes..==
 +
 +
* centerline prepare
 +
<pre>
 +
        # new
 +
        subdiv = slicer.vtkLinearSubdivisionFilter()
 +
        subdiv.SetInput(surfaceTriangulator.GetOutput())
 +
        subdiv.SetNumberOfSubdivisions(1)
 +
        subdiv.Update()
 +
 +
        smooth = slicer.vtkWindowedSincPolyDataFilter()
 +
        smooth.SetInput(subdiv.GetOutput())
 +
        smooth.SetNumberOfIterations(20)
 +
        smooth.SetPassBand(0.1)
 +
        smooth.SetBoundarySmoothing(1)
 +
        smooth.Update()
 +
 +
        normals = slicer.vtkPolyDataNormals()
 +
        normals.SetInput(smooth.GetOutput())
 +
        normals.SetAutoOrientNormals(1)
 +
        normals.SetFlipNormals(0)
 +
        normals.SetConsistency(1)
 +
        normals.SplittingOff()
 +
        normals.Update()
 +
</pre>
  
 
* this code works locally:
 
* this code works locally:
Line 28: Line 52:
 
${OUTPUT}=$const(/Volumes/DATA2/dhaehn/gwe-results/${FILE}-${TYPE}/${FILE}-${TYPE}-out-${SIGMA_MIN}-${SIGMA_MAX}-N${SIGMA_STEPS}-${ALPHA}-${BETA}-${GAMMA})
 
${OUTPUT}=$const(/Volumes/DATA2/dhaehn/gwe-results/${FILE}-${TYPE}/${FILE}-${TYPE}-out-${SIGMA_MIN}-${SIGMA_MAX}-N${SIGMA_STEPS}-${ALPHA}-${BETA}-${GAMMA})
  
−
${SLICER_HOME}/Slicer3 --no_splash --evalpython "import sys; from Slicer import slicer; sys.path.append(str(slicer.Application.GetExtensionsInstallPath())+'/'+str(slicer.Application.GetSvnRevision())+'/VMTKVesselEnhancement/VMTKVesselEnhancement'); from SlicerVMTKVesselEnhancementGUI import *; hiddengui = VMTKVesselEnhancement(); from SlicerVMTKVesselEnhancementLogic import *; vesselness=SlicerVMTKVesselEnhancementLogic(hiddengui); volNode=slicer.VolumesGUI.GetLogic().AddArchetypeVolume('${INPUT}','Vessels',0); matrix = slicer.vtkMatrix4x4(); volNode.GetIJKToRASMatrix(matrix); from datetime import *; s=datetime.now(); outVolumeData = vesselness.ApplyFrangiVesselness(volNode.GetImageData(),${SIGMA_MIN},${SIGMA_MAX},${SIGMA_STEPS},${ALPHA},${BETA},${GAMMA}); e=datetime.now(); delta=e-s; outputNode=slicer.MRMLScene.CreateNodeByClass('vtkMRMLScalarVolumeNode'); volumeNode=slicer.MRMLScene.AddNode(outputNode); volumeNode.SetAndObserveImageData(outVolumeData); volumeNode.SetIJKToRASMatrix(matrix); volumeNode.SetModifiedSinceRead(1); slicer.VolumesGUI.GetLogic().SaveArchetypeVolume('${OUTPUT}--'+str(delta).replace(':','-')+'.nrrd',volumeNode);Slicer.tk.eval('exit 0');"
+
if [ -e ${OUTPUT}*.nrrd ]; then echo '${OUTPUT}*.nrrd exists!'; else ${SLICER_HOME}/bin/Slicer3-real --no_splash --evalpython "import sys; from Slicer import slicer; sys.path.append(str(slicer.Application.GetExtensionsInstallPath())+'/'+str(slicer.Application.GetSvnRevision())+'/VMTKVesselEnhancement/VMTKVesselEnhancement'); from SlicerVMTKVesselEnhancementGUI import *; hiddengui = VMTKVesselEnhancement(); from SlicerVMTKVesselEnhancementLogic import *; vesselness=SlicerVMTKVesselEnhancementLogic(hiddengui); volNode=slicer.VolumesGUI.GetLogic().AddArchetypeVolume('${INPUT}','Vessels',0); matrix = slicer.vtkMatrix4x4(); volNode.GetIJKToRASMatrix(matrix); from datetime import *; s=datetime.now(); outVolumeData = vesselness.ApplyFrangiVesselness(volNode.GetImageData(),${SIGMA_MIN},${SIGMA_MAX},${SIGMA_STEPS},${ALPHA},${BETA},${GAMMA}); e=datetime.now(); delta=e-s; outputNode=slicer.MRMLScene.CreateNodeByClass('vtkMRMLScalarVolumeNode'); volumeNode=slicer.MRMLScene.AddNode(outputNode); volumeNode.SetAndObserveImageData(outVolumeData); volumeNode.SetIJKToRASMatrix(matrix); volumeNode.SetModifiedSinceRead(1); slicer.VolumesGUI.GetLogic().SaveArchetypeVolume('${OUTPUT}--'+str(delta).replace(':','-')+'.nrrd',volumeNode);"; fi;
 
</pre>
 
</pre>
  
−
* remote levelset:
+
* remote levelset plus centerlines and export:
 
<pre>
 
<pre>
 
${SIGMA_MIN}=$const(1.0)
 
${SIGMA_MIN}=$const(1.0)
Line 39: Line 63:
 
${BETA}=$range(500.0,1500.0,500.0)
 
${BETA}=$range(500.0,1500.0,500.0)
 
${GAMMA}=$range(300.0,1000.0,100.0)
 
${GAMMA}=$range(300.0,1000.0,100.0)
 +
${FILE}=$const(01)
 +
${TYPE}=$const(RCA)
 +
${INPUT}=$const(/Volumes/DATA2/dhaehn/cutted/Image${FILE}_${TYPE}.nrrd)
 +
${VESSELNESS_NAME}=$const(${FILE}-${TYPE}-out-${SIGMA_MIN}-${SIGMA_MAX}-N${SIGMA_STEPS}-${ALPHA}-${BETA}-${GAMMA})
 +
${INPUT_VESSELNESS}=$const(/Volumes/DATA2/dhaehn/gwe-results/${FILE}-${TYPE}/${VESSELNESS_NAME}*.nrrd)
 +
${VESSEL}=$const(0)
 +
${SEED_X}=$const(-50.1738)
 +
${SEED_Y}=$const(-58.7432)
 +
${SEED_Z}=$const(77.2189)
 +
${CL_START_X}=$const(-61.3221)
 +
${CL_START_Y}=$const(-70.5911)
 +
${CL_START_Z}=$const(86.7588)
 +
${CL_END_X}=$const(-110.032)
 +
${CL_END_Y}=$const(-123.933)
 +
${CL_END_Z}=$const(34.8108)
 +
${LOW_THRESHOLD}=$range(0.4,0.6,0.1)
 +
${HIGH_THRESHOLD}=$const(1.0)
 +
${PROPAGATION}=$range(-10,10,10)
 +
${CURVATURE}=$range(50,80,10)
 +
${ADVECTION}=$const(100)
 +
${ITERATIONS}=$const(10)
 +
${OUTPUT_A}=$const(/Volumes/DATA2/dhaehn/gwe-results-centerlines/${FILE}-${TYPE}/)
 +
${OUTPUT_B}=$const(-V${VESSEL}-FM-${LOW_THRESHOLD}-${HIGH_THRESHOLD}--)
 +
${OUTPUT_C}=$const(-GAC-${PROPAGATION}-${CURVATURE}-${ADVECTION}-N${ITERATIONS}--)
 +
 +
if [ -e ${OUTPUT_A}${VESSELNESS_NAME}*${OUTPUT_B}*${OUTPUT_C}*.dat ]; then echo '${OUTPUT_A}${VESSELNESS_NAME}*${OUTPUT_B}*${OUTPUT_C}*.dat exists!'; else
 +
${SLICER_HOME}/bin/Slicer3-real --no_splash --evalpython "import sys; from Slicer import slicer;
 +
 +
sys.path.append(str(slicer.Application.GetExtensionsInstallPath())+'/'+str(slicer.Application.GetSvnRevision())+'/VMTKEasyLevelSetSegmentation/VMTKEasyLevelSetSegmentation');
 +
sys.path.append(str(slicer.Application.GetExtensionsInstallPath())+'/'+str(slicer.Application.GetSvnRevision())+'/VMTKCenterlines/VMTKCenterlines');
 +
 +
 +
from VMTKEasyLevelSetSegmentationGUI import *;
 +
 +
hiddengui = VMTKEasyLevelSetSegmentationGUI();
 +
 +
from VMTKEasyLevelSetSegmentationLogic import *;
 +
 +
levelset = VMTKEasyLevelSetSegmentationLogic(hiddengui);
 +
 +
from VMTKCenterlinesGUI import *;
 +
 +
hidden2gui = VMTKCenterlinesGUI();
 +
 +
from VMTKCenterlinesLogic import *;
 +
 +
centerlines = VMTKCenterlinesLogic(hidden2gui);
 +
 +
volNode=slicer.VolumesGUI.GetLogic().AddArchetypeVolume('${INPUT}','Original',0);
 +
 +
import glob;
 +
 +
vesselnessFile=glob.glob('${INPUT_VESSELNESS}');
 +
vesselnessNode=slicer.VolumesGUI.GetLogic().AddArchetypeVolume(vesselnessFile[0],'Vessels',0);
 +
 +
matrix = slicer.vtkMatrix4x4(); volNode.GetIJKToRASMatrix(matrix);
 +
rasMatrix = slicer.vtkMatrix4x4(); volNode.GetRASToIJKMatrix(rasMatrix);
 +
image = volNode.GetImageData();
 +
vesselImage = vesselnessNode.GetImageData();
 +
 +
seedList = slicer.vtkIdList();
 +
targetList = slicer.vtkIdList();
 +
seedPt = [${SEED_X},${SEED_Y},${SEED_Z},1];
 +
 +
ijkSeedPt = rasMatrix.MultiplyPoint(*seedPt);
 +
seedList.InsertNextId(vesselImage.ComputePointId(int(ijkSeedPt[0]),int(ijkSeedPt[1]),int(ijkSeedPt[2])));
 +
 +
from datetime import *;
 +
 +
fm_s = datetime.now();
 +
fm = levelset.ExecuteFM(vesselImage,${LOW_THRESHOLD},${HIGH_THRESHOLD},seedList,targetList);
 +
fm_e = datetime.now();
 +
delta_fm=fm_e-fm_s;
 +
gac_s=datetime.now();
 +
gac = levelset.ExecuteGAC(image,fm,${ITERATIONS},${PROPAGATION},${CURVATURE},${ADVECTION},'geodesic');
 +
gac_e=datetime.now();
 +
delta_gac=gac_e-gac_s;
 +
mc_s=datetime.now();
 +
mc = levelset.MarchingCubes(gac,matrix,0.0);
 +
mc_e=datetime.now();
 +
delta_mc=mc_e-mc_s;
 +
 +
clSeeds = slicer.vtkIdList();
 +
clTargets = slicer.vtkIdList();
 +
clSeedPt = [${CL_START_X},${CL_START_Y},${CL_START_Z},1];
 +
clTargetPt = [${CL_END_X},${CL_END_Y},${CL_END_Z},1];
 +
 +
pointLocator=slicer.vtkPointLocator();
 +
pointLocator.SetDataSet(mc);
 +
pointLocator.BuildLocator();
 +
 +
id=pointLocator.FindClosestPoint(int(clSeedPt[0]),int(clSeedPt[1]),int(clSeedPt[2]));
 +
clSeeds.InsertNextId(id);
 +
id=pointLocator.FindClosestPoint(int(clTargetPt[0]),int(clTargetPt[1]),int(clTargetPt[2]));
 +
clTargets.InsertNextId(id);
 +
 +
cl_s=datetime.now();
 +
prepared = centerlines.prepareModel(mc);
 +
centerline = centerlines.computeCenterlines(prepared,clSeeds,clTargets);
 +
cl_e=datetime.now();
 +
delta_cl=cl_e-cl_s;
 +
vesselnessFile=str(vesselnessFile[0]).replace('.nrrd','').split('/');
 +
outputFile='${OUTPUT_A}'+vesselnessFile[len(vesselnessFile)-1]+'${OUTPUT_B}'+str(delta_fm).replace(':','-')+'${OUTPUT_C}'+str(delta_gac).replace(':','-')+'-MC--'+str(delta_mc).replace(':','-')+'-CL--'+str(delta_cl).replace(':','-')+'.dat'; centerlines.Export(centerline,outputFile,0,0,1);"; fi
 +
</pre>
 +
 +
* remote Frangi 2nd Run
 +
<pre>
 +
${SIGMA_MIN}=$const(1.0)
 +
${SIGMA_MAX}=$const(5.0)
 +
${SIGMA_STEPS}=$const(10)
 +
${ALPHA}=$range(0.1,0.3,0.1)
 +
${BETA}=$const(500.0)
 +
${GAMMA}=$range(200.0,600.0,100.0)
 +
${FILE}=$const(00)
 +
${TYPE}=$const(RCA,LCA)
 +
${INPUT}=$const(/Volumes/DATA2/dhaehn/cutted/Image${FILE}_${TYPE}.nrrd)
 +
${OUTPUT}=$const(/Volumes/DATA2/dhaehn/gwe-results-snd/${FILE}-${TYPE}/${FILE}-${TYPE}-out-${SIGMA_MIN}-${SIGMA_MAX}-N${SIGMA_STEPS}-${ALPHA}-${BETA}-${GAMMA})
 +
${SLICER}=$const(/Users/dhaehn/current_slicer)
 +
 +
if [ -e ${OUTPUT}*.nrrd ]; then echo '${OUTPUT}*.nrrd exists!'; else ${SLICER}/bin/Slicer3-real --no_splash --evalpython "import sys; from Slicer import slicer; sys.path.append(str(slicer.Application.GetExtensionsInstallPath())+'/'+str(slicer.Application.GetSvnRevision())+'/VMTKVesselEnhancement/VMTKVesselEnhancement'); from SlicerVMTKVesselEnhancementGUI import *; hiddengui = VMTKVesselEnhancement(); from SlicerVMTKVesselEnhancementLogic import *; vesselness=SlicerVMTKVesselEnhancementLogic(hiddengui); volNode=slicer.VolumesGUI.GetLogic().AddArchetypeVolume('${INPUT}','Vessels',0); matrix = slicer.vtkMatrix4x4(); volNode.GetIJKToRASMatrix(matrix); from datetime import *; s=datetime.now(); outVolumeData = vesselness.ApplyFrangiVesselness(volNode.GetImageData(),${SIGMA_MIN},${SIGMA_MAX},${SIGMA_STEPS},${ALPHA},${BETA},${GAMMA}); e=datetime.now(); delta=e-s; outputNode=slicer.MRMLScene.CreateNodeByClass('vtkMRMLScalarVolumeNode'); volumeNode=slicer.MRMLScene.AddNode(outputNode); volumeNode.SetAndObserveImageData(outVolumeData); volumeNode.SetIJKToRASMatrix(matrix); volumeNode.SetModifiedSinceRead(1); slicer.VolumesGUI.GetLogic().SaveArchetypeVolume('${OUTPUT}--'+str(delta).replace(':','-')+'.nrrd',volumeNode);"; fi;
 +
</pre>
 +
 +
* remote levelset plus centerlines and export 2nd run:
 +
<pre>
 +
${SIGMA_MIN}=$const(1.0)
 +
${SIGMA_MAX}=$const(5.0)
 +
${SIGMA_STEPS}=$const(10)
 +
${ALPHA}=$range(0.1,0.3,0.1)
 +
${BETA}=$const(500.0)
 +
${GAMMA}=$range(200.0,600.0,100.0)
 
${FILE}=$const(00)
 
${FILE}=$const(00)
 
${TYPE}=$const(RCA)
 
${TYPE}=$const(RCA)
 
${INPUT}=$const(/Volumes/DATA2/dhaehn/cutted/Image${FILE}_${TYPE}.nrrd)
 
${INPUT}=$const(/Volumes/DATA2/dhaehn/cutted/Image${FILE}_${TYPE}.nrrd)
−
${VESSELNESS_NAME}=$const(${FILE}-${TYPE}/${FILE}-${TYPE}-out-${SIGMA_MIN}-${SIGMA_MAX}-N${SIGMA_STEPS}-${ALPHA}-${BETA}-${GAMMA})
+
${VESSELNESS_NAME}=$const(${FILE}-${TYPE}-out-${SIGMA_MIN}-${SIGMA_MAX}-N${SIGMA_STEPS}-${ALPHA}-${BETA}-${GAMMA})
−
${INPUT_VESSELNESS}=$const(/Volumes/DATA2/dhaehn/gwe-results/${VESSELNESS_NAME}*.nrrd)
+
${INPUT_VESSELNESS}=$const(/Volumes/DATA2/dhaehn/gwe-results-snd/${FILE}-${TYPE}/${VESSELNESS_NAME}*.nrrd)
−
${SEED_X}=$const(1)
+
${VESSEL}=$const(0)
−
${SEED_Y}=$const(1)
+
${LOW_THRESHOLD}=$range(0.4,0.6,0.1)
−
${SEED_Z}=$const(1)
 
−
${LOW_THRESHOLD}=$const(0.1)
 
 
${HIGH_THRESHOLD}=$const(1.0)
 
${HIGH_THRESHOLD}=$const(1.0)
−
${PROPAGATION}=$const(0)
+
${PROPAGATION}=$range(-15,15,15)
−
${CURVATURE}=$const(70)
+
${CURVATURE}=$range(40,90,10)
 
${ADVECTION}=$const(100)
 
${ADVECTION}=$const(100)
 
${ITERATIONS}=$const(10)
 
${ITERATIONS}=$const(10)
−
${OUTPUT_A}=$const(/Volumes/DATA2/dhaehn/gwe-results-polydata/${VESSELNESS_NAME}-FM-${LOW_THRESHOLD}-${HIGH_THRESHOLD}--)
+
${OUTPUT_A}=$const(/Volumes/DATA2/dhaehn/gwe-results-centerlines-snd/${FILE}-${TYPE}${VESSEL}/)
−
${OUTPUT_B}=$const(-GAC-${PROPAGATION}-${CURVATURE}-${ADVECTION}-N${ITERATIONS}--)
+
${OUTPUT_B}=$const(-V${VESSEL}-FM-${LOW_THRESHOLD}-${HIGH_THRESHOLD}--)
 +
${OUTPUT_C}=$const(-GAC-${PROPAGATION}-${CURVATURE}-${ADVECTION}-N${ITERATIONS}--)
 +
${SLICER}=$const(/Users/dhaehn/current_slicer)
  
−
{SLICER_HOME}/bin/Slicer3-real --no_splash --evalpython "import sys; from Slicer import slicer;  
+
if [ -e ${OUTPUT_A}${VESSELNESS_NAME}*${OUTPUT_B}*${OUTPUT_C}*.dat ]; then echo '${OUTPUT_A}${VESSELNESS_NAME}*${OUTPUT_B}*${OUTPUT_C}*.dat exists!'; else
 +
if [ -e ${OUTPUT_A}ended_${VESSELNESS_NAME}${OUTPUT_B}${OUTPUT_C}.txt ]; then echo '${OUTPUT_A}ended_${VESSELNESS_NAME}${OUTPUT_B}${OUTPUT_C}.dat exists!'; else
 +
${SLICER}/bin/Slicer3-real --no_splash --evalpython "import sys; from Slicer import slicer;  
  
 
sys.path.append(str(slicer.Application.GetExtensionsInstallPath())+'/'+str(slicer.Application.GetSvnRevision())+'/VMTKEasyLevelSetSegmentation/VMTKEasyLevelSetSegmentation');  
 
sys.path.append(str(slicer.Application.GetExtensionsInstallPath())+'/'+str(slicer.Application.GetSvnRevision())+'/VMTKEasyLevelSetSegmentation/VMTKEasyLevelSetSegmentation');  
 +
sys.path.append(str(slicer.Application.GetExtensionsInstallPath())+'/'+str(slicer.Application.GetSvnRevision())+'/VMTKCenterlines/VMTKCenterlines');
 +
  
 
from VMTKEasyLevelSetSegmentationGUI import *;  
 
from VMTKEasyLevelSetSegmentationGUI import *;  
Line 68: Line 226:
 
levelset = VMTKEasyLevelSetSegmentationLogic(hiddengui);
 
levelset = VMTKEasyLevelSetSegmentationLogic(hiddengui);
  
−
volNode=slicer.VolumesGUI.GetLogic().AddArchetypeVolume('${INPUT}','Originals',0);
+
from VMTKCenterlinesGUI import *;
 +
 
 +
hidden2gui = VMTKCenterlinesGUI();
 +
 
 +
from VMTKCenterlinesLogic import *;
 +
 
 +
centerlines = VMTKCenterlinesLogic(hidden2gui);
 +
 
 +
volNode=slicer.VolumesGUI.GetLogic().AddArchetypeVolume('${INPUT}','Original',0);
  
 
import glob;
 
import glob;
Line 82: Line 248:
 
seedList = slicer.vtkIdList();
 
seedList = slicer.vtkIdList();
 
targetList = slicer.vtkIdList();
 
targetList = slicer.vtkIdList();
−
seedPt = [${SEED_X},${SEED_Y},${SEED_Z},1];
+
 
 +
f=open('/Volumes/DATA2/dhaehn/CTA_training_data/dataset${FILE}/vessel${VESSEL}/pointA.txt','r');
 +
line=f.readline();
 +
coord=line.strip('\r\n').split(' ',3);
 +
coord[0]=float(coord[0])*-1;
 +
coord[1]=float(coord[1])*-1;
 +
coord[2]=float(coord[2]);
 +
coord.append(1);
 +
f.close();
 +
seedPt = coord;
  
 
ijkSeedPt = rasMatrix.MultiplyPoint(*seedPt);
 
ijkSeedPt = rasMatrix.MultiplyPoint(*seedPt);
 
seedList.InsertNextId(vesselImage.ComputePointId(int(ijkSeedPt[0]),int(ijkSeedPt[1]),int(ijkSeedPt[2])));
 
seedList.InsertNextId(vesselImage.ComputePointId(int(ijkSeedPt[0]),int(ijkSeedPt[1]),int(ijkSeedPt[2])));
 +
 +
f=open('/Volumes/DATA2/dhaehn/CTA_training_data/dataset${FILE}/vessel${VESSEL}/pointB.txt','r');
 +
line=f.readline();
 +
coord=line.strip('\r\n').split(' ',3);
 +
coord[0]=float(coord[0])*-1;
 +
coord[1]=float(coord[1])*-1;
 +
coord[2]=float(coord[2]);
 +
coord.append(1);
 +
f.close();
 +
seed2Pt = coord;
 +
 +
ijkSeed2Pt = rasMatrix.MultiplyPoint(*seed2Pt);
 +
seedList.InsertNextId(vesselImage.ComputePointId(int(ijkSeed2Pt[0]),int(ijkSeed2Pt[1]),int(ijkSeed2Pt[2])));
  
 
from datetime import *;
 
from datetime import *;
Line 102: Line 290:
 
delta_mc=mc_e-mc_s;
 
delta_mc=mc_e-mc_s;
  
−
outputNode=slicer.MRMLScene.CreateNodeByClass('vtkMRMLModelNode');  
+
clSeeds = slicer.vtkIdList();
−
modelNode=slicer.MRMLScene.AddNode(outputNode);  
+
clTargets = slicer.vtkIdList();
−
modelNode.SetAndObservePolyData(mc); modelNode.SetModifiedSinceRead(1);  
+
 
 +
f=open('/Volumes/DATA2/dhaehn/CTA_training_data/dataset${FILE}/vessel${VESSEL}/pointS.txt','r');
 +
line=f.readline();
 +
coord=line.strip('\r\n').split(' ',3);
 +
coord[0]=float(coord[0])*-1;
 +
coord[1]=float(coord[1])*-1;
 +
coord[2]=float(coord[2]);
 +
coord.append(1);
 +
f.close();
 +
clSeedPt = coord;
 +
 
 +
f=open('/Volumes/DATA2/dhaehn/CTA_training_data/dataset${FILE}/vessel${VESSEL}/pointE.txt','r');
 +
line=f.readline();
 +
coord=line.strip('\r\n').split(' ',3);
 +
coord[0]=float(coord[0])*-1;
 +
coord[1]=float(coord[1])*-1;
 +
coord[2]=float(coord[2]);
 +
coord.append(1);
 +
f.close();
 +
clTargetPt = coord;
 +
 
 +
pointLocator=slicer.vtkPointLocator();
 +
pointLocator.SetDataSet(mc);
 +
pointLocator.BuildLocator();
 +
 
 +
id=pointLocator.FindClosestPoint(int(clSeedPt[0]),int(clSeedPt[1]),int(clSeedPt[2]));
 +
clSeeds.InsertNextId(id);
 +
id=pointLocator.FindClosestPoint(int(clTargetPt[0]),int(clTargetPt[1]),int(clTargetPt[2]));
 +
clTargets.InsertNextId(id);
  
−
slicer.ModelsGUI.GetLogic().SaveModel('${OUTPUT_A}'+str(delta_fm).replace(':','-')+'${OUTPUT_B}'+str(delta_gac).replace(':','-')+'-MC--'+str(delta_mc).replace(':','-')+'.vtk',modelNode);"
+
cl_s=datetime.now();
 +
prepared = centerlines.prepareModel(mc);
 +
centerline = centerlines.computeCenterlines(prepared,clSeeds,clTargets);
 +
cl_e=datetime.now();
 +
delta_cl=cl_e-cl_s;
 +
vesselnessFile=str(vesselnessFile[0]).replace('.nrrd','').split('/');
 +
outputFile='${OUTPUT_A}${SYSTEM_JOB_NUM}_'+vesselnessFile[len(vesselnessFile)-1]+'${OUTPUT_B}'+str(delta_fm).replace(':','-')+'${OUTPUT_C}'+str(delta_gac).replace(':','-')+'-MC--'+str(delta_mc).replace(':','-')+'-CL--'+str(delta_cl).replace(':','-')+'.dat'; centerlines.Export(centerline,outputFile,0,0,1);"; fi; echo "ended" > ${OUTPUT_A}ended_${VESSELNESS_NAME}${OUTPUT_B}${OUTPUT_C}.txt; fi
 
</pre>
 
</pre>

Latest revision as of 11:24, 29 April 2010

Home < User talk:Haehn

my notes..

  • centerline prepare
        # new
        subdiv = slicer.vtkLinearSubdivisionFilter()
        subdiv.SetInput(surfaceTriangulator.GetOutput())
        subdiv.SetNumberOfSubdivisions(1)
        subdiv.Update()

        smooth = slicer.vtkWindowedSincPolyDataFilter()
        smooth.SetInput(subdiv.GetOutput())
        smooth.SetNumberOfIterations(20)
        smooth.SetPassBand(0.1)
        smooth.SetBoundarySmoothing(1)
        smooth.Update()

        normals = slicer.vtkPolyDataNormals()
        normals.SetInput(smooth.GetOutput())
        normals.SetAutoOrientNormals(1)
        normals.SetFlipNormals(0)
        normals.SetConsistency(1)
        normals.SplittingOff()
        normals.Update()
  • this code works locally:

"import sys; from Slicer import slicer; volNode=slicer.VolumesGUI.GetLogic().AddArchetypeVolume('/home/haehn/data/Image00_LCA.nrrd','Liver',0); matrix = slicer.vtkMatrix4x4(); volNode.GetIJKToRASMatrix(matrix);

outputNode=slicer.MRMLScene.CreateNodeByClass('vtkMRMLScalarVolumeNode'); volumeNode=slicer.MRMLScene.AddNode(outputNode); volumeNode.SetAndObserveImageData(volNode.GetImageData()); volumeNode.SetIJKToRASMatrix(matrix); volumeNode.SetModifiedSinceRead(1); slicer.VolumesGUI.GetLogic().SaveArchetypeVolume('/home/haehn/data/test.nrrd',volumeNode);"

  • remote Frangi:
${SIGMA_MIN}=$const(1.0)
${SIGMA_MAX}=$const(5.0)
${SIGMA_STEPS}=$const(10)
${ALPHA}=$range(0.1,0.5,0.1)
${BETA}=$range(500.0,1500.0,500.0)
${GAMMA}=$range(300.0,1000.0,100.0)
${FILE}=$const(00)
${TYPE}=$const(RCA)
${INPUT}=$const(/Volumes/DATA2/dhaehn/cutted/Image${FILE}_${TYPE}.nrrd)
${OUTPUT}=$const(/Volumes/DATA2/dhaehn/gwe-results/${FILE}-${TYPE}/${FILE}-${TYPE}-out-${SIGMA_MIN}-${SIGMA_MAX}-N${SIGMA_STEPS}-${ALPHA}-${BETA}-${GAMMA})

if [ -e ${OUTPUT}*.nrrd ]; then echo '${OUTPUT}*.nrrd exists!'; else ${SLICER_HOME}/bin/Slicer3-real --no_splash --evalpython "import sys; from Slicer import slicer; sys.path.append(str(slicer.Application.GetExtensionsInstallPath())+'/'+str(slicer.Application.GetSvnRevision())+'/VMTKVesselEnhancement/VMTKVesselEnhancement'); from SlicerVMTKVesselEnhancementGUI import *; hiddengui = VMTKVesselEnhancement(); from SlicerVMTKVesselEnhancementLogic import *; vesselness=SlicerVMTKVesselEnhancementLogic(hiddengui); volNode=slicer.VolumesGUI.GetLogic().AddArchetypeVolume('${INPUT}','Vessels',0); matrix = slicer.vtkMatrix4x4(); volNode.GetIJKToRASMatrix(matrix); from datetime import *; s=datetime.now(); outVolumeData = vesselness.ApplyFrangiVesselness(volNode.GetImageData(),${SIGMA_MIN},${SIGMA_MAX},${SIGMA_STEPS},${ALPHA},${BETA},${GAMMA}); e=datetime.now(); delta=e-s; outputNode=slicer.MRMLScene.CreateNodeByClass('vtkMRMLScalarVolumeNode'); volumeNode=slicer.MRMLScene.AddNode(outputNode); volumeNode.SetAndObserveImageData(outVolumeData); volumeNode.SetIJKToRASMatrix(matrix); volumeNode.SetModifiedSinceRead(1); slicer.VolumesGUI.GetLogic().SaveArchetypeVolume('${OUTPUT}--'+str(delta).replace(':','-')+'.nrrd',volumeNode);"; fi;
  • remote levelset plus centerlines and export:
${SIGMA_MIN}=$const(1.0)
${SIGMA_MAX}=$const(5.0)
${SIGMA_STEPS}=$const(10)
${ALPHA}=$range(0.1,0.5,0.1)
${BETA}=$range(500.0,1500.0,500.0)
${GAMMA}=$range(300.0,1000.0,100.0)
${FILE}=$const(01)
${TYPE}=$const(RCA)
${INPUT}=$const(/Volumes/DATA2/dhaehn/cutted/Image${FILE}_${TYPE}.nrrd)
${VESSELNESS_NAME}=$const(${FILE}-${TYPE}-out-${SIGMA_MIN}-${SIGMA_MAX}-N${SIGMA_STEPS}-${ALPHA}-${BETA}-${GAMMA})
${INPUT_VESSELNESS}=$const(/Volumes/DATA2/dhaehn/gwe-results/${FILE}-${TYPE}/${VESSELNESS_NAME}*.nrrd)
${VESSEL}=$const(0)
${SEED_X}=$const(-50.1738)
${SEED_Y}=$const(-58.7432)
${SEED_Z}=$const(77.2189)
${CL_START_X}=$const(-61.3221)
${CL_START_Y}=$const(-70.5911)
${CL_START_Z}=$const(86.7588)
${CL_END_X}=$const(-110.032)
${CL_END_Y}=$const(-123.933)
${CL_END_Z}=$const(34.8108)
${LOW_THRESHOLD}=$range(0.4,0.6,0.1)
${HIGH_THRESHOLD}=$const(1.0)
${PROPAGATION}=$range(-10,10,10)
${CURVATURE}=$range(50,80,10)
${ADVECTION}=$const(100)
${ITERATIONS}=$const(10)
${OUTPUT_A}=$const(/Volumes/DATA2/dhaehn/gwe-results-centerlines/${FILE}-${TYPE}/)
${OUTPUT_B}=$const(-V${VESSEL}-FM-${LOW_THRESHOLD}-${HIGH_THRESHOLD}--)
${OUTPUT_C}=$const(-GAC-${PROPAGATION}-${CURVATURE}-${ADVECTION}-N${ITERATIONS}--)

if [ -e ${OUTPUT_A}${VESSELNESS_NAME}*${OUTPUT_B}*${OUTPUT_C}*.dat ]; then echo '${OUTPUT_A}${VESSELNESS_NAME}*${OUTPUT_B}*${OUTPUT_C}*.dat exists!'; else 
${SLICER_HOME}/bin/Slicer3-real --no_splash --evalpython "import sys; from Slicer import slicer; 

sys.path.append(str(slicer.Application.GetExtensionsInstallPath())+'/'+str(slicer.Application.GetSvnRevision())+'/VMTKEasyLevelSetSegmentation/VMTKEasyLevelSetSegmentation'); 
sys.path.append(str(slicer.Application.GetExtensionsInstallPath())+'/'+str(slicer.Application.GetSvnRevision())+'/VMTKCenterlines/VMTKCenterlines'); 


from VMTKEasyLevelSetSegmentationGUI import *; 

hiddengui = VMTKEasyLevelSetSegmentationGUI(); 

from VMTKEasyLevelSetSegmentationLogic import *; 

levelset = VMTKEasyLevelSetSegmentationLogic(hiddengui);

from VMTKCenterlinesGUI import *;

hidden2gui = VMTKCenterlinesGUI();

from VMTKCenterlinesLogic import *;

centerlines = VMTKCenterlinesLogic(hidden2gui);

volNode=slicer.VolumesGUI.GetLogic().AddArchetypeVolume('${INPUT}','Original',0);

import glob;

vesselnessFile=glob.glob('${INPUT_VESSELNESS}');
vesselnessNode=slicer.VolumesGUI.GetLogic().AddArchetypeVolume(vesselnessFile[0],'Vessels',0);

matrix = slicer.vtkMatrix4x4(); volNode.GetIJKToRASMatrix(matrix);
rasMatrix = slicer.vtkMatrix4x4(); volNode.GetRASToIJKMatrix(rasMatrix);
image = volNode.GetImageData();
vesselImage = vesselnessNode.GetImageData();

seedList = slicer.vtkIdList();
targetList = slicer.vtkIdList();
seedPt = [${SEED_X},${SEED_Y},${SEED_Z},1];

ijkSeedPt = rasMatrix.MultiplyPoint(*seedPt);
seedList.InsertNextId(vesselImage.ComputePointId(int(ijkSeedPt[0]),int(ijkSeedPt[1]),int(ijkSeedPt[2])));

from datetime import *;

fm_s = datetime.now();
fm = levelset.ExecuteFM(vesselImage,${LOW_THRESHOLD},${HIGH_THRESHOLD},seedList,targetList);
fm_e = datetime.now();
delta_fm=fm_e-fm_s;
gac_s=datetime.now();
gac = levelset.ExecuteGAC(image,fm,${ITERATIONS},${PROPAGATION},${CURVATURE},${ADVECTION},'geodesic');
gac_e=datetime.now();
delta_gac=gac_e-gac_s;
mc_s=datetime.now();
mc = levelset.MarchingCubes(gac,matrix,0.0);
mc_e=datetime.now();
delta_mc=mc_e-mc_s;

clSeeds = slicer.vtkIdList();
clTargets = slicer.vtkIdList();
clSeedPt = [${CL_START_X},${CL_START_Y},${CL_START_Z},1];
clTargetPt = [${CL_END_X},${CL_END_Y},${CL_END_Z},1];

pointLocator=slicer.vtkPointLocator();
pointLocator.SetDataSet(mc);
pointLocator.BuildLocator();

id=pointLocator.FindClosestPoint(int(clSeedPt[0]),int(clSeedPt[1]),int(clSeedPt[2]));
clSeeds.InsertNextId(id);
id=pointLocator.FindClosestPoint(int(clTargetPt[0]),int(clTargetPt[1]),int(clTargetPt[2]));
clTargets.InsertNextId(id);

cl_s=datetime.now();
prepared = centerlines.prepareModel(mc);
centerline = centerlines.computeCenterlines(prepared,clSeeds,clTargets);
cl_e=datetime.now();
delta_cl=cl_e-cl_s;
vesselnessFile=str(vesselnessFile[0]).replace('.nrrd','').split('/'); 
outputFile='${OUTPUT_A}'+vesselnessFile[len(vesselnessFile)-1]+'${OUTPUT_B}'+str(delta_fm).replace(':','-')+'${OUTPUT_C}'+str(delta_gac).replace(':','-')+'-MC--'+str(delta_mc).replace(':','-')+'-CL--'+str(delta_cl).replace(':','-')+'.dat'; centerlines.Export(centerline,outputFile,0,0,1);"; fi
  • remote Frangi 2nd Run
${SIGMA_MIN}=$const(1.0)
${SIGMA_MAX}=$const(5.0)
${SIGMA_STEPS}=$const(10)
${ALPHA}=$range(0.1,0.3,0.1)
${BETA}=$const(500.0)
${GAMMA}=$range(200.0,600.0,100.0)
${FILE}=$const(00)
${TYPE}=$const(RCA,LCA)
${INPUT}=$const(/Volumes/DATA2/dhaehn/cutted/Image${FILE}_${TYPE}.nrrd)
${OUTPUT}=$const(/Volumes/DATA2/dhaehn/gwe-results-snd/${FILE}-${TYPE}/${FILE}-${TYPE}-out-${SIGMA_MIN}-${SIGMA_MAX}-N${SIGMA_STEPS}-${ALPHA}-${BETA}-${GAMMA})
${SLICER}=$const(/Users/dhaehn/current_slicer)

if [ -e ${OUTPUT}*.nrrd ]; then echo '${OUTPUT}*.nrrd exists!'; else ${SLICER}/bin/Slicer3-real --no_splash --evalpython "import sys; from Slicer import slicer; sys.path.append(str(slicer.Application.GetExtensionsInstallPath())+'/'+str(slicer.Application.GetSvnRevision())+'/VMTKVesselEnhancement/VMTKVesselEnhancement'); from SlicerVMTKVesselEnhancementGUI import *; hiddengui = VMTKVesselEnhancement(); from SlicerVMTKVesselEnhancementLogic import *; vesselness=SlicerVMTKVesselEnhancementLogic(hiddengui); volNode=slicer.VolumesGUI.GetLogic().AddArchetypeVolume('${INPUT}','Vessels',0); matrix = slicer.vtkMatrix4x4(); volNode.GetIJKToRASMatrix(matrix); from datetime import *; s=datetime.now(); outVolumeData = vesselness.ApplyFrangiVesselness(volNode.GetImageData(),${SIGMA_MIN},${SIGMA_MAX},${SIGMA_STEPS},${ALPHA},${BETA},${GAMMA}); e=datetime.now(); delta=e-s; outputNode=slicer.MRMLScene.CreateNodeByClass('vtkMRMLScalarVolumeNode'); volumeNode=slicer.MRMLScene.AddNode(outputNode); volumeNode.SetAndObserveImageData(outVolumeData); volumeNode.SetIJKToRASMatrix(matrix); volumeNode.SetModifiedSinceRead(1); slicer.VolumesGUI.GetLogic().SaveArchetypeVolume('${OUTPUT}--'+str(delta).replace(':','-')+'.nrrd',volumeNode);"; fi;
  • remote levelset plus centerlines and export 2nd run:
${SIGMA_MIN}=$const(1.0)
${SIGMA_MAX}=$const(5.0)
${SIGMA_STEPS}=$const(10)
${ALPHA}=$range(0.1,0.3,0.1)
${BETA}=$const(500.0)
${GAMMA}=$range(200.0,600.0,100.0)
${FILE}=$const(00)
${TYPE}=$const(RCA)
${INPUT}=$const(/Volumes/DATA2/dhaehn/cutted/Image${FILE}_${TYPE}.nrrd)
${VESSELNESS_NAME}=$const(${FILE}-${TYPE}-out-${SIGMA_MIN}-${SIGMA_MAX}-N${SIGMA_STEPS}-${ALPHA}-${BETA}-${GAMMA})
${INPUT_VESSELNESS}=$const(/Volumes/DATA2/dhaehn/gwe-results-snd/${FILE}-${TYPE}/${VESSELNESS_NAME}*.nrrd)
${VESSEL}=$const(0)
${LOW_THRESHOLD}=$range(0.4,0.6,0.1)
${HIGH_THRESHOLD}=$const(1.0)
${PROPAGATION}=$range(-15,15,15)
${CURVATURE}=$range(40,90,10)
${ADVECTION}=$const(100)
${ITERATIONS}=$const(10)
${OUTPUT_A}=$const(/Volumes/DATA2/dhaehn/gwe-results-centerlines-snd/${FILE}-${TYPE}${VESSEL}/)
${OUTPUT_B}=$const(-V${VESSEL}-FM-${LOW_THRESHOLD}-${HIGH_THRESHOLD}--)
${OUTPUT_C}=$const(-GAC-${PROPAGATION}-${CURVATURE}-${ADVECTION}-N${ITERATIONS}--)
${SLICER}=$const(/Users/dhaehn/current_slicer)

if [ -e ${OUTPUT_A}${VESSELNESS_NAME}*${OUTPUT_B}*${OUTPUT_C}*.dat ]; then echo '${OUTPUT_A}${VESSELNESS_NAME}*${OUTPUT_B}*${OUTPUT_C}*.dat exists!'; else 
if [ -e ${OUTPUT_A}ended_${VESSELNESS_NAME}${OUTPUT_B}${OUTPUT_C}.txt ]; then echo '${OUTPUT_A}ended_${VESSELNESS_NAME}${OUTPUT_B}${OUTPUT_C}.dat exists!'; else 
${SLICER}/bin/Slicer3-real --no_splash --evalpython "import sys; from Slicer import slicer; 

sys.path.append(str(slicer.Application.GetExtensionsInstallPath())+'/'+str(slicer.Application.GetSvnRevision())+'/VMTKEasyLevelSetSegmentation/VMTKEasyLevelSetSegmentation'); 
sys.path.append(str(slicer.Application.GetExtensionsInstallPath())+'/'+str(slicer.Application.GetSvnRevision())+'/VMTKCenterlines/VMTKCenterlines'); 


from VMTKEasyLevelSetSegmentationGUI import *; 

hiddengui = VMTKEasyLevelSetSegmentationGUI(); 

from VMTKEasyLevelSetSegmentationLogic import *; 

levelset = VMTKEasyLevelSetSegmentationLogic(hiddengui);

from VMTKCenterlinesGUI import *;

hidden2gui = VMTKCenterlinesGUI();

from VMTKCenterlinesLogic import *;

centerlines = VMTKCenterlinesLogic(hidden2gui);

volNode=slicer.VolumesGUI.GetLogic().AddArchetypeVolume('${INPUT}','Original',0);

import glob;

vesselnessFile=glob.glob('${INPUT_VESSELNESS}');
vesselnessNode=slicer.VolumesGUI.GetLogic().AddArchetypeVolume(vesselnessFile[0],'Vessels',0);

matrix = slicer.vtkMatrix4x4(); volNode.GetIJKToRASMatrix(matrix);
rasMatrix = slicer.vtkMatrix4x4(); volNode.GetRASToIJKMatrix(rasMatrix);
image = volNode.GetImageData();
vesselImage = vesselnessNode.GetImageData();

seedList = slicer.vtkIdList();
targetList = slicer.vtkIdList();

f=open('/Volumes/DATA2/dhaehn/CTA_training_data/dataset${FILE}/vessel${VESSEL}/pointA.txt','r');
line=f.readline();
coord=line.strip('\r\n').split(' ',3);
coord[0]=float(coord[0])*-1;
coord[1]=float(coord[1])*-1;
coord[2]=float(coord[2]);
coord.append(1);
f.close();
seedPt = coord;

ijkSeedPt = rasMatrix.MultiplyPoint(*seedPt);
seedList.InsertNextId(vesselImage.ComputePointId(int(ijkSeedPt[0]),int(ijkSeedPt[1]),int(ijkSeedPt[2])));

f=open('/Volumes/DATA2/dhaehn/CTA_training_data/dataset${FILE}/vessel${VESSEL}/pointB.txt','r');
line=f.readline();
coord=line.strip('\r\n').split(' ',3);
coord[0]=float(coord[0])*-1;
coord[1]=float(coord[1])*-1;
coord[2]=float(coord[2]);
coord.append(1);
f.close();
seed2Pt = coord;

ijkSeed2Pt = rasMatrix.MultiplyPoint(*seed2Pt);
seedList.InsertNextId(vesselImage.ComputePointId(int(ijkSeed2Pt[0]),int(ijkSeed2Pt[1]),int(ijkSeed2Pt[2])));

from datetime import *;

fm_s = datetime.now();
fm = levelset.ExecuteFM(vesselImage,${LOW_THRESHOLD},${HIGH_THRESHOLD},seedList,targetList);
fm_e = datetime.now();
delta_fm=fm_e-fm_s;
gac_s=datetime.now();
gac = levelset.ExecuteGAC(image,fm,${ITERATIONS},${PROPAGATION},${CURVATURE},${ADVECTION},'geodesic');
gac_e=datetime.now();
delta_gac=gac_e-gac_s;
mc_s=datetime.now();
mc = levelset.MarchingCubes(gac,matrix,0.0);
mc_e=datetime.now();
delta_mc=mc_e-mc_s;

clSeeds = slicer.vtkIdList();
clTargets = slicer.vtkIdList();

f=open('/Volumes/DATA2/dhaehn/CTA_training_data/dataset${FILE}/vessel${VESSEL}/pointS.txt','r');
line=f.readline();
coord=line.strip('\r\n').split(' ',3);
coord[0]=float(coord[0])*-1;
coord[1]=float(coord[1])*-1;
coord[2]=float(coord[2]);
coord.append(1);
f.close();
clSeedPt = coord;

f=open('/Volumes/DATA2/dhaehn/CTA_training_data/dataset${FILE}/vessel${VESSEL}/pointE.txt','r');
line=f.readline();
coord=line.strip('\r\n').split(' ',3);
coord[0]=float(coord[0])*-1;
coord[1]=float(coord[1])*-1;
coord[2]=float(coord[2]);
coord.append(1);
f.close();
clTargetPt = coord;

pointLocator=slicer.vtkPointLocator();
pointLocator.SetDataSet(mc);
pointLocator.BuildLocator();

id=pointLocator.FindClosestPoint(int(clSeedPt[0]),int(clSeedPt[1]),int(clSeedPt[2]));
clSeeds.InsertNextId(id);
id=pointLocator.FindClosestPoint(int(clTargetPt[0]),int(clTargetPt[1]),int(clTargetPt[2]));
clTargets.InsertNextId(id);

cl_s=datetime.now();
prepared = centerlines.prepareModel(mc);
centerline = centerlines.computeCenterlines(prepared,clSeeds,clTargets);
cl_e=datetime.now();
delta_cl=cl_e-cl_s;
vesselnessFile=str(vesselnessFile[0]).replace('.nrrd','').split('/'); 
outputFile='${OUTPUT_A}${SYSTEM_JOB_NUM}_'+vesselnessFile[len(vesselnessFile)-1]+'${OUTPUT_B}'+str(delta_fm).replace(':','-')+'${OUTPUT_C}'+str(delta_gac).replace(':','-')+'-MC--'+str(delta_mc).replace(':','-')+'-CL--'+str(delta_cl).replace(':','-')+'.dat'; centerlines.Export(centerline,outputFile,0,0,1);"; fi; echo "ended" > ${OUTPUT_A}ended_${VESSELNESS_NAME}${OUTPUT_B}${OUTPUT_C}.txt; fi
Retrieved from "https://www.na-mic.org/w/index.php?title=User_talk:Haehn&oldid=51674"