Package SlicerVMTKLevelSet :: Module SlicerVMTKLevelSetLogic
[hide private]
[frames] | no frames]

Source Code for Module SlicerVMTKLevelSet.SlicerVMTKLevelSetLogic

  1  from Slicer import slicer 
  2   
  3  from SlicerVMTKLevelSetContainer import SlicerVMTKLevelSetContainer 
  4   
  5  ### pseudo logic class 
6 -class SlicerVMTKLevelSetLogic(object):
7
8 - def __init__(self,mainGUIClass):
9 10 self._helper = mainGUIClass.GetHelper() 11 12 # advanced settings, hardcoded... 13 self.SigmoidRemapping = 1 14 self.DerivativeSigma = 0.0 15 self.FeatureDerivativeSigma = 0.0 16 self.MaximumRMSError = 1E-20 17 self.IsoSurfaceValue = 0.0
18
19 - def ExecuteCollidingFronts(self,inVolumeNode,lowerThreshold,higherThreshold,sourceSeedsNode,targetSeedsNode):
20 21 self._helper.debug("Starting execution of Colliding Fronts..") 22 23 if not inVolumeNode or not sourceSeedsNode or not targetSeedsNode: 24 self._helper.debug(inVolumeNode) 25 self._helper.debug(lowerThreshold) 26 self._helper.debug(higherThreshold) 27 self._helper.debug(sourceSeedsNode) 28 self._helper.debug(targetSeedsNode) 29 slicer.Application.ErrorMessage("Not enough information!!! Aborting Colliding Fronts..\n") 30 return 31 else: 32 33 sourceSeedIds = slicer.vtkIdList() 34 targetSeedIds = slicer.vtkIdList() 35 36 image = inVolumeNode.GetImageData() 37 38 cast = slicer.vtkImageCast() 39 cast.SetInput(image) 40 cast.SetOutputScalarTypeToFloat() 41 cast.Update() 42 image = cast.GetOutput() 43 44 rasPt = sourceSeedsNode.GetNthFiducialXYZ(0) 45 self._helper.debug(rasPt) 46 ijkPt = self._helper.ConvertRAS2IJK(rasPt) 47 self._helper.debug(ijkPt) 48 sourceSeedIds.InsertNextId(image.ComputePointId(int(ijkPt[0]),int(ijkPt[1]),int(ijkPt[2]))) 49 50 rasPt = targetSeedsNode.GetNthFiducialXYZ(0) 51 self._helper.debug(rasPt) 52 ijkPt = self._helper.ConvertRAS2IJK(rasPt) 53 self._helper.debug(ijkPt) 54 targetSeedIds.InsertNextId(image.ComputePointId(int(ijkPt[0]),int(ijkPt[1]),int(ijkPt[2]))) 55 56 scalarRange = image.GetScalarRange() 57 self._helper.debug("CF: after converting seeds") 58 59 threshold = slicer.vtkImageThreshold() 60 threshold.SetInput(image) 61 threshold.ThresholdBetween(lowerThreshold, higherThreshold) 62 threshold.ReplaceInOff() 63 threshold.ReplaceOutOn() 64 threshold.SetOutValue(scalarRange[0] - scalarRange[1]) 65 threshold.Update() 66 67 self._helper.debug("CF: after thresholding") 68 69 scalarRange = threshold.GetOutput().GetScalarRange() 70 71 thresholdedImage = threshold.GetOutput() 72 73 shiftScale = slicer.vtkImageShiftScale() 74 shiftScale.SetInput(thresholdedImage) 75 shiftScale.SetShift(-scalarRange[0]) 76 shiftScale.SetScale(1/(scalarRange[1]-scalarRange[0])) 77 shiftScale.Update() 78 79 speedImage = shiftScale.GetOutput() 80 81 self._helper.debug("CF: after shiftScale") 82 83 collidingFronts = slicer.vtkvmtkCollidingFrontsImageFilter() 84 collidingFronts.SetInput(speedImage) 85 collidingFronts.SetSeeds1(sourceSeedIds) 86 collidingFronts.SetSeeds2(targetSeedIds) 87 collidingFronts.ApplyConnectivityOn() 88 collidingFronts.StopOnTargetsOn() 89 collidingFronts.Update() 90 91 self._helper.debug("CF: after CF") 92 93 subtract = slicer.vtkImageMathematics() 94 subtract.SetInput(collidingFronts.GetOutput()) 95 subtract.SetOperationToAddConstant() 96 subtract.SetConstantC(-10.0 * collidingFronts.GetNegativeEpsilon()) 97 subtract.Update() 98 99 self._helper.debug("CF: after substract") 100 101 matrix = slicer.vtkMatrix4x4() 102 inVolumeNode.GetIJKToRASMatrix(matrix) 103 104 outVolumeData = slicer.vtkImageData() 105 outVolumeData.DeepCopy(subtract.GetOutput()) 106 outVolumeData.Update() 107 108 # volume calculated... 109 110 outVolumeNode = slicer.vtkMRMLScalarVolumeNode() 111 outVolumeNode.SetAndObserveImageData(outVolumeData) 112 outVolumeNode.SetIJKToRASMatrix(matrix) 113 114 outputContainer = SlicerVMTKLevelSetContainer(outVolumeNode,collidingFronts.GetNegativeEpsilon()) 115 116 self._helper.debug("Colliding Fronts done...") 117 118 return outputContainer
119 120 # run completed 121
122 - def ExecuteFastMarching(self,inVolumeNode,lowerThreshold,higherThreshold,sourceSeedsNode,targetSeedsNode):
123 124 self._helper.debug("Starting execution of Fast Marching...") 125 126 if not inVolumeNode or not sourceSeedsNode or not targetSeedsNode: 127 self._helper.debug(inVolumeNode) 128 self._helper.debug(lowerThreshold) 129 self._helper.debug(higherThreshold) 130 self._helper.debug(sourceSeedsNode) 131 self._helper.debug(targetSeedsNode) 132 slicer.Application.ErrorMessage("Not enough information!!! Aborting Fast Marching..\n") 133 return 134 else: 135 136 sourceSeedIds = slicer.vtkIdList() 137 targetSeedIds = slicer.vtkIdList() 138 139 image = inVolumeNode.GetImageData() 140 141 cast = slicer.vtkImageCast() 142 cast.SetInput(image) 143 cast.SetOutputScalarTypeToFloat() 144 cast.Update() 145 image = cast.GetOutput() 146 147 for i in range(sourceSeedsNode.GetNumberOfFiducials()): 148 rasPt = sourceSeedsNode.GetNthFiducialXYZ(i) 149 ijkPt = self._helper.ConvertRAS2IJK(rasPt) 150 sourceSeedIds.InsertNextId(image.ComputePointId(int(ijkPt[0]),int(ijkPt[1]),int(ijkPt[2]))) 151 152 for i in range(targetSeedsNode.GetNumberOfFiducials()): 153 rasPt = targetSeedsNode.GetNthFiducialXYZ(i) 154 ijkPt = self._helper.ConvertRAS2IJK(rasPt) 155 targetSeedIds.InsertNextId(image.ComputePointId(int(ijkPt[0]),int(ijkPt[1]),int(ijkPt[2]))) 156 157 scalarRange = image.GetScalarRange() 158 159 threshold = slicer.vtkImageThreshold() 160 threshold.SetInput(image) 161 threshold.ThresholdBetween(lowerThreshold,higherThreshold) 162 threshold.ReplaceInOff() 163 threshold.ReplaceOutOn() 164 threshold.SetOutValue(scalarRange[0] - scalarRange[1]) 165 threshold.Update() 166 167 scalarRange = threshold.GetOutput().GetScalarRange() 168 169 thresholdedImage = threshold.GetOutput() 170 171 shiftScale = slicer.vtkImageShiftScale() 172 shiftScale.SetInput(thresholdedImage) 173 shiftScale.SetShift(-scalarRange[0]) 174 shiftScale.SetScale(1/(scalarRange[1]-scalarRange[0])) 175 shiftScale.Update() 176 177 speedImage = shiftScale.GetOutput() 178 179 fastMarching = slicer.vtkvmtkFastMarchingUpwindGradientImageFilter() 180 fastMarching.SetInput(speedImage) 181 fastMarching.SetSeeds(sourceSeedIds) 182 fastMarching.GenerateGradientImageOff() 183 fastMarching.SetTargetOffset(100.0) 184 if targetSeedIds.GetNumberOfIds() > 0: 185 fastMarching.SetTargets(targetSeedIds) 186 fastMarching.SetTargetReachedModeToOneTarget() 187 else: 188 fastMarching.SetTargetReachedModeToNoTargets() 189 fastMarching.Update() 190 191 subtract = slicer.vtkImageMathematics() 192 subtract.SetInput(fastMarching.GetOutput()) 193 subtract.SetOperationToAddConstant() 194 subtract.SetConstantC(-fastMarching.GetTargetValue()) 195 subtract.Update() 196 197 matrix = slicer.vtkMatrix4x4() 198 inVolumeNode.GetIJKToRASMatrix(matrix) 199 200 outVolumeData = slicer.vtkImageData() 201 outVolumeData.DeepCopy(subtract.GetOutput()) 202 outVolumeData.Update() 203 204 # volume calculated... 205 206 outVolumeNode = slicer.vtkMRMLScalarVolumeNode() 207 outVolumeNode.SetAndObserveImageData(outVolumeData) 208 outVolumeNode.SetIJKToRASMatrix(matrix) 209 210 outputContainer = SlicerVMTKLevelSetContainer(outVolumeNode,0.0) 211 212 self._helper.debug("Fast Marching done...") 213 214 return outputContainer
215 216 # run completed 217 218 219
220 - def ExecuteThreshold(self,inVolumeNode,lowerThreshold,higherThreshold):
221 222 223 self._helper.debug("Starting execution of Threshold...") 224 225 226 if not inVolumeNode: 227 slicer.Application.ErrorMessage("No input volume found. Aborting Threshold..\n") 228 return 229 else: 230 231 image = inVolumeNode.GetImageData() 232 233 cast = slicer.vtkImageCast() 234 cast.SetInput(image) 235 cast.SetOutputScalarTypeToFloat() 236 cast.Update() 237 image = cast.GetOutput() 238 239 scalarRange = image.GetScalarRange() 240 241 242 threshold = slicer.vtkImageThreshold() 243 threshold.SetInput(image) 244 245 threshold.ThresholdBetween(lowerThreshold,higherThreshold) 246 247 threshold.ReplaceInOff() 248 threshold.ReplaceOutOn() 249 threshold.SetInValue(-1.0) 250 threshold.SetOutValue(1.0) 251 threshold.Update() 252 253 matrix = slicer.vtkMatrix4x4() 254 inVolumeNode.GetIJKToRASMatrix(matrix) 255 256 outVolumeData = slicer.vtkImageData() 257 outVolumeData.DeepCopy(threshold.GetOutput()) 258 outVolumeData.Update() 259 260 # volume calculated... 261 262 outVolumeNode = slicer.vtkMRMLScalarVolumeNode() 263 outVolumeNode.SetAndObserveImageData(outVolumeData) 264 outVolumeNode.SetIJKToRASMatrix(matrix) 265 266 outputContainer = SlicerVMTKLevelSetContainer(outVolumeNode,0.0) 267 268 self._helper.debug("Threshold done...") 269 270 return outputContainer
271 272 # run completed 273
274 - def ExecuteIsosurface(self,inVolumeNode,value):
275 276 277 self._helper.debug("Starting execution of Isosurface...") 278 279 if not inVolumeNode: 280 slicer.Application.ErrorMessage("No input volume found. Aborting Isosurface..\n") 281 return 282 else: 283 284 image = inVolumeNode.GetImageData() 285 286 cast = slicer.vtkImageCast() 287 cast.SetInput(image) 288 cast.SetOutputScalarTypeToFloat() 289 cast.Update() 290 image = cast.GetOutput() 291 292 imageMathematics = slicer.vtkImageMathematics() 293 imageMathematics.SetInput(image) 294 imageMathematics.SetConstantK(-1.0) 295 imageMathematics.SetOperationToMultiplyByK() 296 imageMathematics.Update() 297 298 subtract = slicer.vtkImageMathematics() 299 subtract.SetInput(imageMathematics.GetOutput()) 300 subtract.SetOperationToAddConstant() 301 subtract.SetConstantC(value) 302 subtract.Update() 303 304 matrix = slicer.vtkMatrix4x4() 305 inVolumeNode.GetIJKToRASMatrix(matrix) 306 307 outVolumeData = slicer.vtkImageData() 308 outVolumeData.DeepCopy(subtract.GetOutput()) 309 310 outVolumeData.Update() 311 312 # volume calculated... 313 314 outVolumeNode = slicer.vtkMRMLScalarVolumeNode() 315 outVolumeNode.SetAndObserveImageData(outVolumeData) 316 outVolumeNode.SetIJKToRASMatrix(matrix) 317 318 outputContainer = SlicerVMTKLevelSetContainer(outVolumeNode,10) 319 320 self._helper.debug("Isosurface done...") 321 322 return outputContainer
323 324 # run completed 325 326
327 - def ExecuteSeeds(self,inVolumeNode,sourceSeedsNode):
328 329 330 self._helper.debug("Starting execution of seeds...") 331 332 if not inVolumeNode or not sourceSeedsNode: 333 slicer.Application.ErrorMessage("Not enough information. Aborting Seeds..\n") 334 return 335 else: 336 337 image = inVolumeNode.GetImageData() 338 339 cast = slicer.vtkImageCast() 340 cast.SetInput(image) 341 cast.SetOutputScalarTypeToFloat() 342 cast.Update() 343 image = cast.GetOutput() 344 345 seeds = sourceSeedsNode 346 347 initialLevelSets = slicer.vtkImageData() 348 initialLevelSets.DeepCopy(image) 349 initialLevelSets.Update() 350 351 levelSetsInputScalars = initialLevelSets.GetPointData().GetScalars() 352 levelSetsInputScalars.FillComponent(0,1.0) 353 354 for i in range(seeds.GetNumberOfFiducials()): 355 rasPt = seeds.GetNthFiducialXYZ(i) 356 ijkPt = self._helper.ConvertRAS2IJK(rasPt) 357 id = image.ComputePointId(int(ijkPt[0]),int(ijkPt[1]),int(ijkPt[2])) 358 levelSetsInputScalars.SetComponent(id,0,-1.0) 359 360 dilateErode = slicer.vtkImageDilateErode3D() 361 dilateErode.SetInput(initialLevelSets) 362 dilateErode.SetDilateValue(-1.0) 363 dilateErode.SetErodeValue(1.0) 364 dilateErode.SetKernelSize(3,3,3) 365 dilateErode.Update() 366 367 matrix = slicer.vtkMatrix4x4() 368 inVolumeNode.GetIJKToRASMatrix(matrix) 369 370 outVolumeData = slicer.vtkImageData() 371 outVolumeData.DeepCopy(dilateErode.GetOutput()) 372 outVolumeData.Update() 373 374 # volume calculated... 375 376 outVolumeNode = slicer.vtkMRMLScalarVolumeNode() 377 outVolumeNode.SetAndObserveImageData(outVolumeData) 378 outVolumeNode.SetIJKToRASMatrix(matrix) 379 380 outputContainer = SlicerVMTKLevelSetContainer(outVolumeNode,0.0) 381 382 self._helper.debug("Seeds done...") 383 384 return outputContainer
385 386 # run completed 387 388 389 390
391 - def BuildGradientBasedFeatureImage(self,imageData):
392 393 cast = slicer.vtkImageCast() 394 cast.SetInput(imageData) 395 cast.SetOutputScalarTypeToFloat() 396 cast.Update() 397 398 if (self.DerivativeSigma > 0.0): 399 gradientMagnitude = slicer.vtkvmtkGradientMagnitudeRecursiveGaussianImageFilter() 400 gradientMagnitude.SetInput(cast.GetOutput()) 401 gradientMagnitude.SetSigma(self.DerivativeSigma) 402 gradientMagnitude.SetNormalizeAcrossScale(0) 403 gradientMagnitude.Update() 404 else: 405 gradientMagnitude = slicer.vtkvmtkGradientMagnitudeImageFilter() 406 gradientMagnitude.SetInput(cast.GetOutput()) 407 gradientMagnitude.Update() 408 409 featureImage = None 410 if self.SigmoidRemapping==1: 411 scalarRange = gradientMagnitude.GetOutput().GetPointData().GetScalars().GetRange() 412 inputMinimum = scalarRange[0] 413 inputMaximum = scalarRange[1] 414 alpha = - (inputMaximum - inputMinimum) / 6.0 415 beta = (inputMaximum + inputMinimum) / 2.0 416 sigmoid = slicer.vtkvmtkSigmoidImageFilter() 417 sigmoid.SetInput(gradientMagnitude.GetOutput()) 418 sigmoid.SetAlpha(alpha) 419 sigmoid.SetBeta(beta) 420 sigmoid.SetOutputMinimum(0.0) 421 sigmoid.SetOutputMaximum(1.0) 422 sigmoid.Update() 423 featureImage = sigmoid.GetOutput() 424 else: 425 boundedReciprocal = slicer.vtkvmtkBoundedReciprocalImageFilter() 426 boundedReciprocal.SetInput(gradientMagnitude.GetOutput()) 427 boundedReciprocal.Update() 428 featureImage = boundedReciprocal.GetOutput() 429 430 featureImageOutput = slicer.vtkImageData() 431 featureImageOutput.DeepCopy(featureImage) 432 featureImageOutput.Update() 433 434 return featureImageOutput
435 436
437 - def ExecuteGeodesic(self,origInVolumeNode,inVolumeNode,numberOfIterations,propagationScaling,curvatureScaling,advectionScaling):
438 439 self._helper.debug("Starting execution of Geodesic..") 440 441 if not inVolumeNode or not origInVolumeNode: 442 slicer.Application.ErrorMessage("Not enough information!!! Aborting Geodesic..\n") 443 return 444 else: 445 446 origImage = origInVolumeNode.GetImageData() 447 image = inVolumeNode.GetImageData() 448 449 levelSets = slicer.vtkvmtkGeodesicActiveContourLevelSetImageFilter() 450 levelSets.SetFeatureImage(self.BuildGradientBasedFeatureImage(origImage)) 451 levelSets.SetDerivativeSigma(self.FeatureDerivativeSigma) 452 levelSets.SetAutoGenerateSpeedAdvection(1) 453 levelSets.SetPropagationScaling(propagationScaling) 454 levelSets.SetCurvatureScaling(curvatureScaling) 455 levelSets.SetAdvectionScaling(advectionScaling) 456 457 458 levelSets.SetInput(image) 459 levelSets.SetNumberOfIterations(numberOfIterations) 460 levelSets.SetIsoSurfaceValue(self.IsoSurfaceValue) 461 levelSets.SetMaximumRMSError(self.MaximumRMSError) 462 levelSets.SetInterpolateSurfaceLocation(1) 463 levelSets.SetUseImageSpacing(1) 464 levelSets.Update() 465 466 467 matrix = slicer.vtkMatrix4x4() 468 inVolumeNode.GetIJKToRASMatrix(matrix) 469 470 outVolumeData = slicer.vtkImageData() 471 outVolumeData.DeepCopy(levelSets.GetOutput()) 472 outVolumeData.Update() 473 474 outVolumeNode = slicer.vtkMRMLScalarVolumeNode() 475 outVolumeNode.SetAndObserveImageData(outVolumeData) 476 outVolumeNode.SetIJKToRASMatrix(matrix) 477 478 outputContainer = SlicerVMTKLevelSetContainer(outVolumeNode,0.0) 479 480 return outputContainer
481 482
483 - def ExecuteCurves(self,origInVolumeNode,inVolumeNode,numberOfIterations,propagationScaling,curvatureScaling,advectionScaling):
484 485 self._helper.debug("Starting execution of Curves..") 486 487 if not inVolumeNode or not origInVolumeNode: 488 slicer.Application.ErrorMessage("Not enough information!!! Aborting Geodesic..\n") 489 return 490 else: 491 492 origImage = origInVolumeNode.GetImageData() 493 image = inVolumeNode.GetImageData() 494 495 levelSets = slicer.vtkvmtkCurvesLevelSetImageFilter() 496 levelSets.SetFeatureImage(self.BuildGradientBasedFeatureImage(origImage)) 497 levelSets.SetDerivativeSigma(self.FeatureDerivativeSigma) 498 levelSets.SetAutoGenerateSpeedAdvection(1) 499 levelSets.SetPropagationScaling(propagationScaling) 500 levelSets.SetCurvatureScaling(curvatureScaling) 501 levelSets.SetAdvectionScaling(advectionScaling) 502 503 504 levelSets.SetInput(image) 505 levelSets.SetNumberOfIterations(numberOfIterations) 506 levelSets.SetIsoSurfaceValue(self.IsoSurfaceValue) 507 levelSets.SetMaximumRMSError(self.MaximumRMSError) 508 levelSets.SetInterpolateSurfaceLocation(1) 509 levelSets.SetUseImageSpacing(1) 510 levelSets.Update() 511 512 513 matrix = slicer.vtkMatrix4x4() 514 inVolumeNode.GetIJKToRASMatrix(matrix) 515 516 outVolumeData = slicer.vtkImageData() 517 outVolumeData.DeepCopy(levelSets.GetOutput()) 518 outVolumeData.Update() 519 520 outVolumeNode = slicer.vtkMRMLScalarVolumeNode() 521 outVolumeNode.SetAndObserveImageData(outVolumeData) 522 outVolumeNode.SetIJKToRASMatrix(matrix) 523 524 outputContainer = SlicerVMTKLevelSetContainer(outVolumeNode,0.0) 525 526 return outputContainer
527 528 529 ### 530 ### marching cubes algorithm, takes imageData and returns polyData 531 ###
532 - def MarchingCubes(self,image,matrix,threshold):
533 534 transformIJKtoRAS = slicer.vtkTransform() 535 transformIJKtoRAS.SetMatrix(matrix) 536 537 marchingCubes = slicer.vtkMarchingCubes() 538 marchingCubes.SetInput(image) 539 marchingCubes.SetValue(0,threshold) 540 marchingCubes.ComputeScalarsOn() 541 marchingCubes.ComputeGradientsOn() 542 marchingCubes.ComputeNormalsOn() 543 marchingCubes.GetOutput().ReleaseDataFlagOn() 544 marchingCubes.Update() 545 546 if transformIJKtoRAS.GetMatrix().Determinant() < 0: 547 self.debug("Determinant smaller than 0") 548 reverser = slicer.vtkReverseSense() 549 reverser.SetInput(marchingCubes.GetOutput()) 550 reverser.ReverseNormalsOn() 551 reverser.GetOutput().ReleaseDataFlagOn() 552 reverser.Update() 553 correctedOutput = reverser.GetOutput() 554 else: 555 correctedOutput = marchingCubes.GetOutput() 556 557 transformer = slicer.vtkTransformPolyDataFilter() 558 transformer.SetInput(correctedOutput) 559 transformer.SetTransform(transformIJKtoRAS) 560 transformer.GetOutput().ReleaseDataFlagOn() 561 transformer.Update() 562 563 normals = slicer.vtkPolyDataNormals() 564 normals.ComputePointNormalsOn() 565 normals.SetInput(transformer.GetOutput()) 566 normals.SetFeatureAngle(60) 567 normals.SetSplitting(1) 568 normals.GetOutput().ReleaseDataFlagOn() 569 normals.Update() 570 571 stripper = slicer.vtkStripper() 572 stripper.SetInput(normals.GetOutput()) 573 stripper.GetOutput().ReleaseDataFlagOff() 574 stripper.Update() 575 stripper.GetOutput().Update() 576 577 result = slicer.vtkPolyData() 578 result.DeepCopy(stripper.GetOutput()) 579 result.Update() 580 581 return result
582