1 from Slicer import slicer
2
3 from SlicerVMTKLevelSetContainer import SlicerVMTKLevelSetContainer
4
5
7
9
10 self._helper = mainGUIClass.GetHelper()
11
12
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
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
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
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
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
217
218
219
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
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
273
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
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
325
326
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
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
387
388
389
390
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
531
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