diff --git a/applications/pctbinning/pctbinning.py b/applications/pctbinning/pctbinning.py
index ce291ee0..80f286e0 100644
--- a/applications/pctbinning/pctbinning.py
+++ b/applications/pctbinning/pctbinning.py
@@ -22,6 +22,7 @@ def build_parser():
     parser.add_argument(
         "-c", "--count", help="Image of count of proton pairs per pixel"
     )
+    parser.add_argument("--variance", help="Image of variance of WEPL values per pixel")
     parser.add_argument(
         "--scatwepl", help="Image of scattering WEPL of proton pairs per pixel"
     )
@@ -33,6 +34,11 @@ def build_parser():
         action="store_true",
         default=False,
     )
+    parser.add_argument(
+        "--particle",
+        help="Particle used for imaging (proton or alpha)",
+        default="proton",
+    )
     parser.add_argument(
         "-v", "--verbose", help="Verbose execution", action="store_true", default=False
     )
@@ -88,6 +94,12 @@ def build_parser():
         action="store_true",
         default=False,
     )
+    parser.add_argument(
+        "--fillvariance",
+        help="Fill holes, i.e. pixels that were not hit by protons (for variance projections)",
+        action="store_true",
+        default=False,
+    )
     parser.add_argument(
         "--trackerresolution", help="Tracker resolution in mm", type=float
     )
@@ -159,6 +171,8 @@ def process(args_info: argparse.Namespace):
     projection.SetRobust(args_info.robust)
     projection.SetComputeScattering(bool(args_info.scatwepl))
     projection.SetComputeNoise(bool(args_info.noise))
+    projection.SetComputeVariance(bool(args_info.variance))
+    projection.SetParticle(args_info.particle)
 
     if args_info.quadricIn:
         # quadric = object surface
@@ -190,15 +204,14 @@ def process(args_info: argparse.Namespace):
 
     projection.Update()
 
-    filler = pct.SmallHoleFiller[OutputImageType]()
+    hole_filter = pct.HoleFillingImageFilter[OutputImageType].New()
     if args_info.fill:
-        filler.SetImage(projection.GetOutput())
-        filler.SetHolePixel(0.0)
-        filler.Fill()
+        hole_filter.SetInput(projection.GetOutput())
+        hole_filter.Update()
 
     cii = itk.ChangeInformationImageFilter[OutputImageType].New()
     if args_info.fill:
-        cii.SetInput(filler.GetOutput())
+        cii.SetInput(hole_filter.GetOutput())
     else:
         cii.SetInput(projection.GetOutput())
     cii.ChangeOriginOn()
@@ -217,6 +230,25 @@ def process(args_info: argparse.Namespace):
         # Write
         itk.imwrite(projection.GetCount(), args_info.count)
 
+    if args_info.variance:
+        VarianceImageType = itk.Image[itk.F, 3]
+        variance_hole_filter = pct.HoleFillingImageFilter[VarianceImageType].New()
+        if args_info.fillvariance:
+            variance_hole_filter.SetInput(projection.GetVariance())
+            variance_hole_filter.Update()
+        cii_var = itk.ChangeInformationImageFilter[VarianceImageType].New()
+        if args_info.fillvariance:
+            cii_var.SetInput(variance_hole_filter.GetOutput())
+        else:
+            cii_var.SetInput(projection.GetVariance())
+        cii_var.ChangeOriginOn()
+        cii_var.ChangeDirectionOn()
+        cii_var.ChangeSpacingOn()
+        cii_var.SetOutputDirection(projection.GetOutput().GetDirection())
+        cii_var.SetOutputOrigin(projection.GetOutput().GetOrigin())
+        cii_var.SetOutputSpacing(projection.GetOutput().GetSpacing())
+        itk.imwrite(cii_var.GetOutput(), args_info.variance)
+
     if args_info.scatwepl:
         # Write
         itk.imwrite(projection.GetAngle(), args_info.scatwepl)