diff --git a/Modules/Core/SpatialObjects/include/itkGaussianSpatialObject.h b/Modules/Core/SpatialObjects/include/itkGaussianSpatialObject.h
index 9632a5f126d..bb68433ad66 100644
--- a/Modules/Core/SpatialObjects/include/itkGaussianSpatialObject.h
+++ b/Modules/Core/SpatialObjects/include/itkGaussianSpatialObject.h
@@ -29,13 +29,13 @@ namespace itk
  *
  * The Gaussian function G(x) is given by
  * \f[
- * G(\vec{x}) = m e^{-\|\S^{-1} \vec{x}\|^2 / 2},
+ * G(\vec{x}) = m e^{-\|S^{-1} \vec{x}\|^2 / 2},
  * \f]
- * where m is a scaling factor set by SetMaximum(), and \f$\S\f$ is the
+ * where m is a scaling factor set by SetMaximum(), and \f$S\f$ is the
  * (invertible) matrix associated to the IndexToObjectTransform of the object
- * multiplied by the Sigma parameter.  If \f$\S\f$ is symmetric and positive
+ * multiplied by the Sigma parameter.  If \f$S\f$ is symmetric and positive
  * definite, and m is chosen so that the integral of G(x) is 1, then G will
- * denote a normal distribution with mean 0 and covariance matrix \f$\S \times
+ * denote a normal distribution with mean 0 and covariance matrix \f$S \times
  * Sigma\f$.
  * \ingroup ITKSpatialObjects
  */
@@ -117,7 +117,7 @@ class ITK_TEMPLATE_EXPORT GaussianSpatialObject : public SpatialObject<TDimensio
                        unsigned int        depth = 0,
                        const std::string & name = "") const override;
 
-  /** Returns the $sigma = $ \c m_Radius level set of the Gaussian function, as an
+  /** Returns the sigma = \c m_Radius level set of the Gaussian function, as an
    * EllipseSpatialObject. */
   typename EllipseSpatialObject<TDimension>::Pointer
   GetEllipsoid() const;
diff --git a/Modules/Core/Transform/include/itkElasticBodySplineKernelTransform.h b/Modules/Core/Transform/include/itkElasticBodySplineKernelTransform.h
index fc9da5d0b0a..ce425ff1785 100644
--- a/Modules/Core/Transform/include/itkElasticBodySplineKernelTransform.h
+++ b/Modules/Core/Transform/include/itkElasticBodySplineKernelTransform.h
@@ -93,7 +93,6 @@ class ITK_TEMPLATE_EXPORT ElasticBodySplineKernelTransform : public KernelTransf
   using typename Superclass::GMatrixType;
   /** Compute G(x)
    * For the elastic body spline, this is:
-   * \f$ G(x) = [alpha*r(x)^2*I - 3*x*x']*r(x) \f$
    * \f$ G(x) = [\alpha*r(x)^2*I - 3*x*x']*r(x) \f$
    * where
    * \f$\alpha = 12 ( 1 - \nu ) - 1\f$
diff --git a/Modules/Core/Transform/include/itkScaleSkewVersor3DTransform.h b/Modules/Core/Transform/include/itkScaleSkewVersor3DTransform.h
index 1a84f9580b3..24e6ff40720 100644
--- a/Modules/Core/Transform/include/itkScaleSkewVersor3DTransform.h
+++ b/Modules/Core/Transform/include/itkScaleSkewVersor3DTransform.h
@@ -46,7 +46,7 @@ namespace itk
  * where \f$\textbf{R}_v\f$ is the rotation matrix given the versor,
  * \f$S=\left( \begin{array}{ccc}s_0-1 & 0 & 0 \\ 0 & s_1-1 & 0 \\ 0 & 0 & s_2-1 \end{array} \right) \f$
  * , and
- * \f$K=\left( \begin{array}{ccc}0 & k_0 & k_1 \\ k_2 & 0 & k_3 \\ k_4 & k_5 & 0 \end{array} \right)\ \f$.
+ * \f$K=\left( \begin{array}{ccc}0 & k_0 & k_1 \\ k_2 & 0 & k_3 \\ k_4 & k_5 & 0 \end{array} \right) \f$.
  *
  * \ingroup ITKTransform
  */
diff --git a/Modules/Core/Transform/include/itkScaleVersor3DTransform.h b/Modules/Core/Transform/include/itkScaleVersor3DTransform.h
index 65f9d879d86..ec7bd2d90b0 100644
--- a/Modules/Core/Transform/include/itkScaleVersor3DTransform.h
+++ b/Modules/Core/Transform/include/itkScaleVersor3DTransform.h
@@ -30,7 +30,7 @@ namespace itk
  * The transform can be described as:
  * \f$ (\textbf{R}_v + \textbf{S})\textbf{x} \f$ where \f$\textbf{R}_v\f$ is the
  * rotation matrix given the versor, and
- * \f$S=\left( \begin{array}{ccc}s_0-1 & 0 & 0 \\ 0 & s_1-1 & 0 \\ 0 & 0 & s_2-1 \end{array} \right)\ \f$
+ * \f$S=\left( \begin{array}{ccc}s_0-1 & 0 & 0 \\ 0 & s_1-1 & 0 \\ 0 & 0 & s_2-1 \end{array} \right) \f$
  *
  *
  * \note This transform's scale parameters are not related to the
diff --git a/Modules/Filtering/Deconvolution/include/itkInverseDeconvolutionImageFilter.h b/Modules/Filtering/Deconvolution/include/itkInverseDeconvolutionImageFilter.h
index 3ba06d8c6b8..6d2c555b77c 100644
--- a/Modules/Filtering/Deconvolution/include/itkInverseDeconvolutionImageFilter.h
+++ b/Modules/Filtering/Deconvolution/include/itkInverseDeconvolutionImageFilter.h
@@ -32,7 +32,7 @@ namespace itk
  * equivalent to multiplying the Fourier transform of the two images,
  * the inverse filter consists of inverting the multiplication. In
  * other words, this filter computes the following:
- * \f[ hat{F}(\omega) =
+ * \f[ \hat{F}(\omega) =
  *       \begin{cases}
  *         G(\omega) / H(\omega) & \text{if $|H(\omega)| \geq \epsilon$} \\
  *         0                     & \text{otherwise}
diff --git a/Modules/Filtering/FastMarching/include/itkFastMarchingNumberOfElementsStoppingCriterion.h b/Modules/Filtering/FastMarching/include/itkFastMarchingNumberOfElementsStoppingCriterion.h
index 7304e1c275e..88f8fb15bb3 100644
--- a/Modules/Filtering/FastMarching/include/itkFastMarchingNumberOfElementsStoppingCriterion.h
+++ b/Modules/Filtering/FastMarching/include/itkFastMarchingNumberOfElementsStoppingCriterion.h
@@ -31,7 +31,7 @@ namespace itk
  *
  * \note For itk::Image, one element is one pixel. So the number of elements is directly
  * linked to the physical size of the object, i.e.
- * \f$ PhysicalSize = TargetNumberOfElements \cdot \prod_{i=1}{dim} Spacing_{i} \f$
+ * \f$ PhysicalSize = TargetNumberOfElements \cdot \prod_{i=1}^{dim} Spacing_{i} \f$
  *
  * \note For itk::QuadEdgeMesh, one element is one vertex.
  *
diff --git a/Modules/Filtering/ImageNoise/include/itkSaltAndPepperNoiseImageFilter.h b/Modules/Filtering/ImageNoise/include/itkSaltAndPepperNoiseImageFilter.h
index 45c08150ae3..c1014a8219b 100644
--- a/Modules/Filtering/ImageNoise/include/itkSaltAndPepperNoiseImageFilter.h
+++ b/Modules/Filtering/ImageNoise/include/itkSaltAndPepperNoiseImageFilter.h
@@ -44,7 +44,7 @@ namespace itk
  * \begin{cases}
  * M,   & \quad \text{if } U < p/2 \\
  * m,   & \quad \text{if } U > 1 - p/2 \\
- * I_0, & \quad \text{if } p/2 \geq U \leq 1 - p/2
+ * I_0, & \quad \text{if } p/2 \geq U \leq 1 - p/2
  * \end{cases} \f$
  *
  * \par
diff --git a/Modules/Filtering/ImageNoise/include/itkSpeckleNoiseImageFilter.h b/Modules/Filtering/ImageNoise/include/itkSpeckleNoiseImageFilter.h
index 810fbb8f410..648fcf18172 100644
--- a/Modules/Filtering/ImageNoise/include/itkSpeckleNoiseImageFilter.h
+++ b/Modules/Filtering/ImageNoise/include/itkSpeckleNoiseImageFilter.h
@@ -35,14 +35,14 @@ namespace itk
  * It can be modeled as:
  *
  * \par
- * \f$ I = I_0 \ast G \f$
+ * \f$ I = I_0 \ast G \f$
  *
  * \par
  * where \f$ G \f$ is a is a gamma distributed random variable of mean 1 and
  * variance proportional to the noise level:
  *
  * \par
- * \f$ G \sim \Gamma(\frac{1}{\sigma^2}, \sigma^2) \f$
+ * \f$ G \sim \Gamma(\frac{1}{\sigma^2}, \sigma^2) \f$
  *
  * \author Gaetan Lehmann
  *
diff --git a/Modules/Filtering/QuadEdgeMeshFiltering/include/itkSmoothingQuadEdgeMeshFilter.h b/Modules/Filtering/QuadEdgeMeshFiltering/include/itkSmoothingQuadEdgeMeshFilter.h
index 84c5e6f2454..9005730f8c0 100644
--- a/Modules/Filtering/QuadEdgeMeshFiltering/include/itkSmoothingQuadEdgeMeshFilter.h
+++ b/Modules/Filtering/QuadEdgeMeshFiltering/include/itkSmoothingQuadEdgeMeshFilter.h
@@ -33,7 +33,7 @@ namespace itk
  *
  * For one iteration the location of one vertex is computed as follows:
  * \f[
- * \boldsymbol{ v' }_i = v_i + m_RelaxationFactor \cdot \frac{ \sum_j w_{ij} ( \boldsymbol{ v_j } - \boldsymbol{ v_i } )
+ * \boldsymbol{ v' }_i = v_i + m_{RelaxationFactor} \cdot \frac{ \sum_j w_{ij} ( \boldsymbol{ v_j } - \boldsymbol{ v_i } )
  * }{ \sum_j w_{ij} } \f]
  *
  * where \f$ w_{ij} \f$ is computed by the means of the set functor
diff --git a/Modules/Numerics/Optimizers/include/itkLBFGSOptimizer.h b/Modules/Numerics/Optimizers/include/itkLBFGSOptimizer.h
index 238063b0ebc..620a3028767 100644
--- a/Modules/Numerics/Optimizers/include/itkLBFGSOptimizer.h
+++ b/Modules/Numerics/Optimizers/include/itkLBFGSOptimizer.h
@@ -43,9 +43,9 @@ namespace itk
  * The step size \f$ s \f$ is determined through line search with the approach
  * by More and Thuente \cite more1994. This line search approach finds a step size such that
  * \f[
- * \lVert \nabla f(x + s (\nabla^2 f(x_n) )^{-1} \nabla f(x) ) \rVert
+ * \| \nabla f(x + s (\nabla^2 f(x_n) )^{-1} \nabla f(x) ) \|
  *   \le
- * \nu \lVert \nabla f(x) \rVert
+ * \nu \| \nabla f(x) \|
  * \f]
  * The parameter \f$ \nu \f$ is set through SetLineSearchAccuracy() (default 0.9)
  * The default step length, i.e. starting step length for the line search,
@@ -53,7 +53,7 @@ namespace itk
  *
  * The optimization stops when either the gradient satisfies the condition
  * \f[
- * \lVert \nabla f(x) \rVert \le \epsilon \max(1, \lVert X \rVert)
+ * \| \nabla f(x) \| \le \epsilon \max(1, \| X \|)
  * \f]
  * or a maximum number of function evaluations has been reached.
  * The tolerance \f$\epsilon\f$ is set through SetGradientConvergenceTolerance()
diff --git a/Modules/Numerics/Optimizersv4/include/itkLBFGS2Optimizerv4.h b/Modules/Numerics/Optimizersv4/include/itkLBFGS2Optimizerv4.h
index e47e224c97d..68191cb93e3 100644
--- a/Modules/Numerics/Optimizersv4/include/itkLBFGS2Optimizerv4.h
+++ b/Modules/Numerics/Optimizersv4/include/itkLBFGS2Optimizerv4.h
@@ -105,9 +105,9 @@ extern ITKOptimizersv4_EXPORT std::ostream &
  * the approach by More and Thuente  \cite more1994. This line search approach finds a step
  * size such that
  * \f[
- * \lVert \nabla f(x + s (\nabla^2 f(x_n) )^{-1} \nabla f(x) ) \rVert
+ * \| \nabla f(x + s (\nabla^2 f(x_n) )^{-1} \nabla f(x) ) \|
  *   \le
- * \nu \lVert \nabla f(x) \rVert
+ * \nu \| \nabla f(x) \|
  * \f]
  * The parameter \f$\nu\f$ is set through SetLineSearchAccuracy() (default 0.9)
  * and SetGradientLineSearchAccuracy()
@@ -120,7 +120,7 @@ extern ITKOptimizersv4_EXPORT std::ostream &
  *
  * The optimization stops when either the gradient satisfies the condition
  * \f[
- * \lVert \nabla f(x) \rVert \le \epsilon \max(1, \lVert X \rVert)
+ * \| \nabla f(x) \| \le \epsilon \max(1, \| X \|)
  * \f]
  * or a maximum number of function evaluations has been reached.
  * The tolerance \f$\epsilon\f$ is set through SetSolutionAccuracy()
@@ -380,7 +380,7 @@ class ITK_TEMPLATE_EXPORT LBFGS2Optimizerv4Template
    *  problems. Setting this parameter to a positive value activates
    *  Orthant-Wise Limited-memory Quasi-Newton (OWL-QN) method, which
    *  minimizes the objective function F(x) combined with the L1 norm |x|
-   *  of the variables, \f$F(x) + C |x|}. \f$. This parameter is the coefficient
+   *  of the variables, \f$F(x) + C |x| \f$. This parameter is the coefficient
    *  for the |x|, i.e., C. As the L1 norm |x| is not differentiable at
    *  zero, the library modifies function and gradient evaluations from
    *  a client program suitably; a client program thus have only to return
diff --git a/Modules/Numerics/Optimizersv4/include/itkLBFGSOptimizerv4.h b/Modules/Numerics/Optimizersv4/include/itkLBFGSOptimizerv4.h
index 3a9b9d2a624..462434c3ced 100644
--- a/Modules/Numerics/Optimizersv4/include/itkLBFGSOptimizerv4.h
+++ b/Modules/Numerics/Optimizersv4/include/itkLBFGSOptimizerv4.h
@@ -42,9 +42,9 @@ namespace itk
  * The step size \f$ s \f$ is determined through line search with the approach
  * by More and Thuente [4]. This line search approach finds a step size such that
  * \f[
- * \lVert \nabla f(x + s (\nabla^2 f(x_n) )^{-1} \nabla f(x) ) \rVert
+ * \| \nabla f(x + s (\nabla^2 f(x_n) )^{-1} \nabla f(x) ) \|
  *   \le
- * \nu \lVert \nabla f(x) \rVert
+ * \nu \| \nabla f(x) \|
  * \f]
  * The parameter \f$ \nu \f$ is set through SetLineSearchAccuracy() (default 0.9)
  * The default step length, i.e. starting step length for the line search,
@@ -52,7 +52,7 @@ namespace itk
  *
  * The optimization stops when either the gradient satisfies the condition
  * \f[
- * \lVert \nabla f(x) \rVert \le \epsilon \max(1, \lVert X \rVert)
+ * \| \nabla f(x) \| \le \epsilon \max(1, \| X \|)
  * \f]
  * or a maximum number of function evaluations has been reached.
  * The tolerance \f$\epsilon\f$ is set through SetGradientConvergenceTolerance()
diff --git a/Modules/Registration/Metricsv4/include/itkCorrelationImageToImageMetricv4HelperThreader.h b/Modules/Registration/Metricsv4/include/itkCorrelationImageToImageMetricv4HelperThreader.h
index 09659141245..f55725b69b4 100644
--- a/Modules/Registration/Metricsv4/include/itkCorrelationImageToImageMetricv4HelperThreader.h
+++ b/Modules/Registration/Metricsv4/include/itkCorrelationImageToImageMetricv4HelperThreader.h
@@ -29,8 +29,8 @@ namespace itk
  * \brief Helper class for CorrelationImageToImageMetricv4 \c
  * To compute the average pixel intensities of the fixed image and the moving image
  * on the sampled points or inside the virtual image region:
- * \f$ \bar f (CorrelationImageToImageMetricv4::m_AverageFix ) \f$
- * \f$ \bar m (CorrelationImageToImageMetricv4::m_AverageMov ) \f$
+ * \f$ \bar f (CorrelationImageToImageMetricv4::m\_AverageFix ) \f$
+ * \f$ \bar m (CorrelationImageToImageMetricv4::m\_AverageMov ) \f$
  *
  * \ingroup ITKMetricsv4
  */
diff --git a/Modules/Segmentation/KLMRegionGrowing/include/itkKLMSegmentationRegion.h b/Modules/Segmentation/KLMRegionGrowing/include/itkKLMSegmentationRegion.h
index ffb277ee203..b44c5428901 100644
--- a/Modules/Segmentation/KLMRegionGrowing/include/itkKLMSegmentationRegion.h
+++ b/Modules/Segmentation/KLMRegionGrowing/include/itkKLMSegmentationRegion.h
@@ -55,7 +55,7 @@ namespace itk
  * surrounded by four borders.
  *
  * Initial regions of a 8 by 9 image with a 4 by 3 grid partition.
- * \f[\begin{tabular}{|c|c|c|c|c|c|c|c|c|}
+ * \f[\begin{array}{|c|c|c|c|c|c|c|c|c|}
  * \hline
  *   1 & 1 & 1 & 2 & 2 & 2 & 3 & 3 & 3 \\ \hline
  *   1 & 1 & 1 & 2 & 2 & 2 & 3 & 3 & 3 \\ \hline
@@ -65,10 +65,10 @@ namespace itk
  *   7 & 7 & 7 & 8 & 8 & 8 & 9 & 9 & 9 \\ \hline
  *   a & a & a & b & b & b & c & c & c \\ \hline
  *   a & a & a & b & b & b & c & c & c \\ \hline
- *  \end{tabular}\f]
+ *  \end{array}\f]
  *
  * Region borders are shown as "E".
- * \f[\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|}
+ * \f[\begin{array}{|c|c|c|c|c|c|c|c|c|c|c|}
  * \hline
  *   C & C & C &   & C & C & C &   & C & C & C \\ \hline
  *   C & C & C & E & C & C & C & E & C & C & C \\ \hline
@@ -81,7 +81,7 @@ namespace itk
  *     & E &   &   &   & E &   &   &   & E &   \\ \hline
  *   C & C & C &   & C & C & C &   & C & C & C \\ \hline
  *   C & C & C & E & C & C & C & E & C & C & C \\ \hline
- *  \end{tabular}\f]
+ *  \end{array}\f]
  *
  * \ingroup RegionGrowingSegmentation
  * \ingroup ITKKLMRegionGrowing
@@ -206,11 +206,11 @@ class ITKKLMRegionGrowing_EXPORT KLMSegmentationRegion : public SegmentationRegi
    * Lambda set to -1.0.
    *
    * For example, take an image with 3 regions A, B, C
-   * \f[\begin{tabular}{|c|c|}
+   * \f[\begin{array}{|c|c|}
    * \hline
    *   A & A \\ \hline
    *   B & C \\ \hline
-   *  \end{tabular}\f]
+   *  \end{array}\f]
    * where region A has region borders A-B and A-C;
    * region B has region borders A-B and B-C; and
    * region C has region borders A-C and B-C.
diff --git a/Modules/Segmentation/LevelSets/include/itkThresholdSegmentationLevelSetFunction.h b/Modules/Segmentation/LevelSets/include/itkThresholdSegmentationLevelSetFunction.h
index b806da8379e..9d9773c30c2 100644
--- a/Modules/Segmentation/LevelSets/include/itkThresholdSegmentationLevelSetFunction.h
+++ b/Modules/Segmentation/LevelSets/include/itkThresholdSegmentationLevelSetFunction.h
@@ -46,7 +46,7 @@ namespace itk
  *  \f$ U \f$ and lower threshold \f$ L \f$ according to the following formula.
  *
  * \par
- *  \f$  f(x) = \left\{ \begin{array}{ll} g(x) - L & \mbox{if $(g)x < (U-L)/2 + L$} \\ U - g(x) & \mbox{otherwise}
+ *  \f$  f(x) = \left\{ \begin{array}{ll} g(x) - L & \mbox{if g(x) < (U-L)/2 + L} \\ U - g(x) & \mbox{otherwise}
  * \end{array} \right. \f$
  *
  * \sa SegmentationLevelSetImageFunction
diff --git a/Modules/Segmentation/MarkovRandomFieldsClassifiers/include/itkMRFImageFilter.h b/Modules/Segmentation/MarkovRandomFieldsClassifiers/include/itkMRFImageFilter.h
index 32ee3764c3f..76429245bb9 100644
--- a/Modules/Segmentation/MarkovRandomFieldsClassifiers/include/itkMRFImageFilter.h
+++ b/Modules/Segmentation/MarkovRandomFieldsClassifiers/include/itkMRFImageFilter.h
@@ -97,23 +97,23 @@ extern ITKMarkovRandomFieldsClassifiers_EXPORT std::ostream &
  * assigned a weight 1.5.  A weight of 1.3 is assigned to the influence of
  * the north, south, east, west, and diagonal pixels in the previous and next
  * slices.
- * \f[\begin{tabular}{ccc}
- *  \begin{tabular}{|c|c|c|}
+ * \f[\begin{array}{ccc}
+ *  \begin{array}{|c|c|c|}
  *   1.3 & 1.3 & 1.3 \\
  *   1.3 & 1.5 & 1.3 \\
  *   1.3 & 1.3 & 1.3 \\
- *  \end{tabular} &
- *  \begin{tabular}{|c|c|c|}
+ *  \end{array} &
+ *  \begin{array}{|c|c|c|}
  *   1.7 & 1.7 & 1.7 \\
  *   1.7 & 0 & 1.7 \\
  *   1.7 & 1.7 & 1.7 \\
- *  \end{tabular} &
- *  \begin{tabular}{|c|c|c|}
+ *  \end{array} &
+ *  \begin{array}{|c|c|c|}
  *   1.3 & 1.3 & 1.3 \\
  *   1.5 & 1.5 & 1.3 \\
  *   1.3 & 1.3 & 1.3 \\
- *  \end{tabular} \\
- * \end{tabular}\f]
+ *  \end{array} \\
+ * \end{array}\f]
  *
  * The user needs to set the neighborhood size using the SetNeighborhoodRadius
  * function. The details on the semantics of a neighborhood can be found
diff --git a/Utilities/Doxygen/DoxygenConfig.cmake b/Utilities/Doxygen/DoxygenConfig.cmake
index 2d4e62d2739..2c9dc2e2ceb 100644
--- a/Utilities/Doxygen/DoxygenConfig.cmake
+++ b/Utilities/Doxygen/DoxygenConfig.cmake
@@ -74,6 +74,7 @@ set(DOXYGEN_DOCSET_PUBLISHER_NAME "InsightConsortium")
 set(DOXYGEN_ECLIPSE_DOC_ID "org.itk.ITK")
 set(DOXYGEN_ENUM_VALUES_PER_LINE "1")
 set(DOXYGEN_USE_MATHJAX "YES")
+set(DOXYGEN_MATHJAX_VERSION "MathJax_2")
 set(DOXYGEN_MATHJAX_RELPATH "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.7/")
 if(ITK_DOXYGEN_SERVER_BASED_SEARCH)
   set(DOXYGEN_SERVER_BASED_SEARCH "YES")