diff --git a/SMDS_3/doc/SMDS_3/PackageDescription.txt b/SMDS_3/doc/SMDS_3/PackageDescription.txt
index 86db00e87e24..91535252f988 100644
--- a/SMDS_3/doc/SMDS_3/PackageDescription.txt
+++ b/SMDS_3/doc/SMDS_3/PackageDescription.txt
@@ -74,6 +74,5 @@
 - `CGAL::IO::output_to_tetgen()`
 - `CGAL::IO::save_binary_file()`
 - `CGAL::IO::load_binary_file()`
-- `CGAL::IO::output_to_medit()` (deprecated)
 */
 
diff --git a/SMDS_3/include/CGAL/IO/File_medit.h b/SMDS_3/include/CGAL/IO/File_medit.h
index 84d11d7d9f54..38296e820244 100644
--- a/SMDS_3/include/CGAL/IO/File_medit.h
+++ b/SMDS_3/include/CGAL/IO/File_medit.h
@@ -790,29 +790,14 @@ output_to_medit(std::ostream& os,
 
 namespace IO {
 
-/**
- * @ingroup PkgSMDS3IOFunctions
- * @deprecated This function is deprecated. Users should instead use `CGAL::IO::write_MEDIT()`
- * @brief outputs a mesh complex to the medit (`.mesh`) file format.
-        See \cgalCite{frey:inria-00069921} for a comprehensive description of this file format.
- * @param os the output stream
- * @param c3t3 the mesh complex
- * @param renumber_subdomain_indices if `true`, labels of cells are renumbered into `[1..nb_of_labels]`
- * @param show_patches if `true`, patches are labeled with different labels than
- *                     cells. If `false`, each surface facet is written twice,
- *                     using the label of each adjacent cell.
- * \see \ref IOStreamMedit
- */
 template <class C3T3>
 void
 output_to_medit(std::ostream& os,
                 const C3T3& c3t3,
                 bool renumber_subdomain_indices, // = false,
                 bool show_patches // = false
-#ifndef DOXYGEN_RUNNING
               , bool all_vertices // = true
               , bool all_cells    // = false
-#endif
 )
 {
   using namespace CGAL::SMDS_3;
@@ -1038,9 +1023,6 @@ bool read_MEDIT(std::istream& in,
 
 } // namespace IO
 
-#ifndef CGAL_NO_DEPRECATED_CODE
-using IO::output_to_medit;
-#endif
 
 } // end namespace CGAL
 
diff --git a/SMDS_3/include/CGAL/SMDS_3/tet_soup_to_c3t3.h b/SMDS_3/include/CGAL/SMDS_3/tet_soup_to_c3t3.h
index c545f9ef70f8..9129628d92a6 100644
--- a/SMDS_3/include/CGAL/SMDS_3/tet_soup_to_c3t3.h
+++ b/SMDS_3/include/CGAL/SMDS_3/tet_soup_to_c3t3.h
@@ -20,6 +20,7 @@
 #include <CGAL/license/SMDS_3.h>
 
 #include <CGAL/assertions.h>
+#include <CGAL/IO/MEDIT.h>
 #include <CGAL/IO/File_medit.h>
 
 #include <boost/unordered_map.hpp>
@@ -585,22 +586,12 @@ bool build_triangulation_from_file(std::istream& is,
   using Facet        = std::array<int, 3>; // 3 = id
   using Tet_with_ref = std::array<int, 4>; // 4 = id
 
-  if(!is)
-    return false;
-
   std::vector<Tet_with_ref> finite_cells;
   std::vector<Subdomain_index> subdomains;
   std::vector<Point_3> points;
   boost::unordered_map<Facet, Surface_patch_index> border_facets;
 
-  int dim;
-  int nv, nf, ntet, ref;
-  std::string word;
-
-  is >> word >> dim; // MeshVersionFormatted 1
-  is >> word >> dim; // Dimension 3
-
-  CGAL_assertion(dim == 3);
+  bool is_CGAL_mesh = false;
 
   if(verbose)
   {
@@ -609,167 +600,11 @@ bool build_triangulation_from_file(std::istream& is,
     std::cout << "Allow non-manifoldness = " << allow_non_manifold << std::endl;
   }
 
-  bool is_CGAL_mesh = false;
-
-  std::string line;
-  while(std::getline(is, line) && line != "End")
-  {
-    // remove trailing whitespace, in particular a possible '\r' from Windows
-    // end-of-line encoding
-    while(!line.empty() && std::isspace(line.back())) {
-      line.pop_back();
-    }
-    if(line.empty())
-      continue;
-
-    // remove whitespaces at the beginning of the line
-    for (std::size_t i=0; i<line.size(); ++i)
-    {
-      if (!std::isspace(line[i]))
-      {
-        if (i!=0)
-          line = line.substr(i);
-        break;
-      }
-    }
-
-    if (line.at(0) == '#' &&
-        line.find("CGAL::Mesh_complex_3_in_triangulation_3") != std::string::npos)
-    {
-      is_CGAL_mesh = true; // with CGAL meshes, domain 0 should be kept
-      continue;
-    }
-
-    // skip non-CGAL comments
-    if (line.at(0)=='#') continue;
-
-    if(line.find("Vertices") != std::string::npos)
-    {
-      is >> nv;
-      if(verbose)
-        std::cerr << "Reading "<< nv << " vertices" << std::endl;
-      for(int i=0; i<nv; ++i)
-      {
-        typename Tr::Geom_traits::FT x,y,z;
-        if(!(is >> x >> y >> z >> ref))
-        {
-          if(verbose)
-            std::cerr << "Issue while reading vertices" << std::endl;
-          return false;
-        }
-        points.emplace_back(x,y,z);
-      }
-    }
-
-    if(line.find("Triangles") != std::string::npos)
-    {
-      bool has_negative_surface_patch_ids = false;
-      Surface_patch_index max_surface_patch_id{0};
-      is >> nf;
-
-      if(verbose)
-        std::cerr << "Reading "<< nf << " triangles" << std::endl;
-
-      for(int i=0; i<nf; ++i)
-      {
-        int n[3];
-        Surface_patch_index surface_patch_id;
-        if(!(is >> n[0] >> n[1] >> n[2] >> surface_patch_id))
-        {
-          if(verbose)
-            std::cerr << "Issue while reading triangles" << std::endl;
-          return false;
-        }
-        has_negative_surface_patch_ids |= (surface_patch_id < 0);
-        max_surface_patch_id = (std::max)(max_surface_patch_id, surface_patch_id);
-        Facet facet;
-        facet[0] = n[0] - 1;
-        facet[1] = n[1] - 1;
-        facet[2] = n[2] - 1;
-
-        if(verbose)
-          std::cout << "Looking at face #" << i << ": " << n[0] << " " << n[1] << " " << n[2] << std::endl;
-
-        CGAL_warning_code(
-        for(int j=0; j<3; ++j)
-          for(int k=0; k<3; ++k)
-            if(j != k)
-              CGAL_warning(n[j] != n[k]);
-        )
-
-        // find the circular permutation that puts the smallest index in the first place.
-        int n0 = (std::min)({facet[0],facet[1], facet[2]});
-        do
-        {
-          std::rotate(std::begin(facet), std::next(std::begin(facet)), std::end(facet));
-        }
-        while(facet[0] != n0);
-
-        border_facets.emplace(facet, surface_patch_id);
-      }
-      if(has_negative_surface_patch_ids)
-      {
-        if(verbose)
-          std::cerr << "Warning: negative surface patch ids" << std::endl;
-        for(auto& facet_and_patch_id  : border_facets) {
-          if(facet_and_patch_id.second < 0)
-            facet_and_patch_id.second = max_surface_patch_id - facet_and_patch_id.second;
-        }
-      }
-    }
-
-    if(line.find("Tetrahedra") != std::string::npos)
-    {
-      is >> ntet;
-
-      if(verbose)
-        std::cerr << "Reading "<< ntet << " tetrahedra" << std::endl;
-
-      for(int i=0; i<ntet; ++i)
-      {
-        int n[4];
-        int reference;
-
-        if(!(is >> n[0] >> n[1] >> n[2] >> n[3] >> reference))
-        {
-          if(verbose)
-            std::cerr << "Issue while reading tetrahedra" << std::endl;
-          return false;
-        }
-
-        if(verbose)
-          std::cout << "Looking at tet #" << i << ": " << n[0] << " " << n[1] << " " << n[2] << " " << n[3] << std::endl;
-
-        CGAL_warning_code(
-        for(int j=0; j<4; ++j)
-          for(int k=0; k<4; ++k)
-            if(j != k)
-              CGAL_warning(n[j] != n[k]);
-        )
-
-        Tet_with_ref t;
-        t[0] = n[0] - 1;
-        t[1] = n[1] - 1;
-        t[2] = n[2] - 1;
-        t[3] = n[3] - 1;
-
-        finite_cells.push_back(t);
-        subdomains.push_back(reference);
-      }
-    }
-  }
-
-  if (verbose)
-  {
-    std::cout << points.size() << " points" << std::endl;
-    std::cout << border_facets.size() << " border facets" << std::endl;
-    std::cout << finite_cells.size() << " cells" << std::endl;
-  }
+  bool ok = CGAL::IO::internal::read_MEDIT(is, points, finite_cells, subdomains, border_facets, true, verbose, is_CGAL_mesh);
 
-  if(finite_cells.empty())
+  if(! ok){
     return false;
-
-  CGAL_assertion(finite_cells.size() == subdomains.size());
+  }
 
   return build_triangulation_with_subdomains_range(tr,
                                                    points, finite_cells, subdomains, border_facets,
diff --git a/Stream_support/doc/Stream_support/PackageDescription.txt b/Stream_support/doc/Stream_support/PackageDescription.txt
index 7875657b1e57..0b7b347e23ae 100644
--- a/Stream_support/doc/Stream_support/PackageDescription.txt
+++ b/Stream_support/doc/Stream_support/PackageDescription.txt
@@ -44,6 +44,10 @@
 /// I/O Functions for the \ref IOStreamXYZ
 /// \ingroup IOstreamFunctions
 
+/// \defgroup PkgStreamSupportIoFuncsMEDIT MEDIT I/O Functions
+/// I/O Functions for the \ref IOStreamMedit
+/// \ingroup IOstreamFunctions
+
 /// \defgroup PkgStreamSupportEnumRef I/O Enums
 /// \ingroup PkgStreamSupportRef
 
@@ -125,5 +129,6 @@ the printing mode.
 - \link PkgStreamSupportIoFuncsWKT I/O for WKT files \endlink
 - \link PkgStreamSupportIoFuncsLAS I/O for LAS files \endlink
 - \link PkgStreamSupportIoFuncsXYZ I/O for XYZ files \endlink
+- \link PkgStreamSupportIoFuncsMEDIT I/O for MEDIT files \endlink
 
 */
diff --git a/Stream_support/include/CGAL/IO/MEDIT.h b/Stream_support/include/CGAL/IO/MEDIT.h
new file mode 100644
index 000000000000..217e74bb1bf3
--- /dev/null
+++ b/Stream_support/include/CGAL/IO/MEDIT.h
@@ -0,0 +1,283 @@
+// Copyright (c) 2015-2020  Geometry Factory
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org)
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s) : Mael Rouxel-Labbé, Andreas Fabri
+
+#ifndef CGAL_IO_MEDIT_H
+#define CGAL_IO_MEDIT_H
+
+#include <CGAL/assertions.h>
+#include <CGAL/Kernel_traits.h>
+#include <iostream>
+#include <string>
+#include <tuple>
+#include <type_traits>
+#include <utility>
+#include <vector>
+#include <array>
+#include <boost/unordered_map.hpp>
+
+namespace CGAL {
+
+////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////
+// Read
+
+namespace IO {
+namespace internal {
+
+template<class PointRange, class CellRange, class SurfacePatchIndex_>
+bool read_MEDIT(std::istream& is,
+                PointRange& points,
+                CellRange& finite_cells,
+                std::vector<int>& subdomains,
+                boost::unordered_map<std::array<int,3>,SurfacePatchIndex_ >& border_facets,
+                bool read_border_facets,
+                bool verbose,
+                bool& is_CGAL_mesh)
+{
+  using Point_3 = typename PointRange::value_type;
+  using FT = typename Kernel_traits<Point_3>::Kernel::FT;
+  using Surface_patch_index = SurfacePatchIndex_;
+  using Facet        = std::array<int, 3>;
+  using Tet_with_ref = std::array<int, 4>;
+
+  if(!is)
+    return false;
+
+  int dim;
+  int nv, nf, ntet, ref;
+  std::string word;
+
+  is >> word >> dim; // MeshVersionFormatted 1
+  is >> word >> dim; // Dimension 3
+
+  CGAL_assertion(dim == 3);
+
+  std::string line;
+  while(std::getline(is, line) && line != "End")
+  {
+    // remove trailing whitespace, in particular a possible '\r' from Windows
+    // end-of-line encoding
+    while(!line.empty() && std::isspace(line.back())) {
+      line.pop_back();
+    }
+    if(line.empty())
+      continue;
+
+    // remove whitespaces at the beginning of the line
+    for (std::size_t i=0; i<line.size(); ++i)
+    {
+      if (!std::isspace(line[i]))
+      {
+        if (i!=0)
+          line = line.substr(i);
+        break;
+      }
+    }
+
+    if (line.at(0) == '#' &&
+        line.find("CGAL::Mesh_complex_3_in_triangulation_3") != std::string::npos)
+    {
+      is_CGAL_mesh = true; // with CGAL meshes, domain 0 should be kept
+      continue;
+    }
+
+    // skip non-CGAL comments
+    if (line.at(0)=='#') continue;
+
+    if(line.find("Vertices") != std::string::npos)
+    {
+      is >> nv;
+      if(verbose)
+        std::cerr << "Reading "<< nv << " vertices" << std::endl;
+      for(int i=0; i<nv; ++i)
+      {
+        FT x,y,z;
+        if(!(is >> x >> y >> z >> ref))
+        {
+          if(verbose)
+            std::cerr << "Issue while reading vertices" << std::endl;
+          return false;
+        }
+        points.emplace_back(x,y,z);
+      }
+    }
+
+    if(line.find("Triangles") != std::string::npos)
+    {
+      if(read_border_facets){
+        bool has_negative_surface_patch_ids = false;
+        Surface_patch_index max_surface_patch_id{0};
+        is >> nf;
+
+        if(verbose)
+          std::cerr << "Reading "<< nf << " triangles" << std::endl;
+
+        for(int i=0; i<nf; ++i)
+        {
+          int n[3];
+          Surface_patch_index surface_patch_id;
+          if(!(is >> n[0] >> n[1] >> n[2] >> surface_patch_id))
+          {
+            if(verbose)
+              std::cerr << "Issue while reading triangles" << std::endl;
+            return false;
+          }
+          has_negative_surface_patch_ids |= (surface_patch_id < 0);
+          max_surface_patch_id = (std::max)(max_surface_patch_id, surface_patch_id);
+          Facet facet;
+          facet[0] = n[0] - 1;
+          facet[1] = n[1] - 1;
+          facet[2] = n[2] - 1;
+
+          if(verbose)
+            std::cout << "Looking at face #" << i << ": " << n[0] << " " << n[1] << " " << n[2] << std::endl;
+
+          CGAL_warning_code(
+          for(int j=0; j<3; ++j)
+            for(int k=0; k<3; ++k)
+              if(j != k)
+                CGAL_warning(n[j] != n[k]);
+          )
+
+          // find the circular permutation that puts the smallest index in the first place.
+          int n0 = (std::min)({facet[0],facet[1], facet[2]});
+          do
+          {
+            std::rotate(std::begin(facet), std::next(std::begin(facet)), std::end(facet));
+          }
+          while(facet[0] != n0);
+
+          border_facets.emplace(facet, surface_patch_id);
+        }
+        if(has_negative_surface_patch_ids)
+        {
+          if(verbose)
+            std::cerr << "Warning: negative surface patch ids" << std::endl;
+          for(auto& facet_and_patch_id  : border_facets) {
+            if(facet_and_patch_id.second < 0)
+              facet_and_patch_id.second = max_surface_patch_id - facet_and_patch_id.second;
+          }
+        }
+      }else{
+        is >> nf;
+        for(int i=0; i<nf; ++i)
+        {
+          int n[3];
+          int surface_patch_id;
+          if(!(is >> n[0] >> n[1] >> n[2] >> surface_patch_id))
+          {
+            if(verbose)
+              std::cerr << "Issue while reading triangles" << std::endl;
+            return false;
+          }
+        }
+      }
+    }
+    if(line.find("Tetrahedra") != std::string::npos)
+    {
+      is >> ntet;
+
+      if(verbose)
+        std::cerr << "Reading "<< ntet << " tetrahedra" << std::endl;
+
+      for(int i=0; i<ntet; ++i)
+      {
+        int n[4];
+        int reference;
+
+        if(!(is >> n[0] >> n[1] >> n[2] >> n[3] >> reference))
+        {
+          if(verbose)
+            std::cerr << "Issue while reading tetrahedra" << std::endl;
+          return false;
+        }
+
+        if(verbose)
+          std::cout << "Looking at tet #" << i << ": " << n[0] << " " << n[1] << " " << n[2] << " " << n[3] << std::endl;
+
+        CGAL_warning_code(
+        for(int j=0; j<4; ++j)
+          for(int k=0; k<4; ++k)
+            if(j != k)
+              CGAL_warning(n[j] != n[k]);
+        )
+
+        Tet_with_ref t;
+        t[0] = n[0] - 1;
+        t[1] = n[1] - 1;
+        t[2] = n[2] - 1;
+        t[3] = n[3] - 1;
+
+        finite_cells.push_back(t);
+        subdomains.push_back(reference);
+      }
+    }
+  }
+
+  if (verbose)
+  {
+    std::cout << points.size() << " points" << std::endl;
+    std::cout << border_facets.size() << " border facets" << std::endl;
+    std::cout << finite_cells.size() << " cells" << std::endl;
+  }
+
+  if(finite_cells.empty())
+    return false;
+
+  CGAL_assertion(finite_cells.size() == subdomains.size());
+
+  return true;
+}
+
+} // namespace internal
+
+
+
+/*!
+ * \ingroup PkgStreamSupportIoFuncsMEDIT
+ *
+ * \brief reads the content of `is` into `points` and `finite_cells`.
+ *
+ * \attention The cell soup is not cleared, and the data from the stream are appended.
+ *
+ * \tparam PointRange a model of the concept `BackInsertionSequence` whose value type is the point type
+ * \tparam CellRange a model of the concept `BackInsertionSequence` whose `value_type` is `std::array<int,4>`
+ *
+ * \param is the input stream
+ * \param points points of the soup of cells
+ * \param finite_cells Each element in it describes a cell
+ *        using the indices of the points in `points`
+ *  \param subdomains Each element in it describes the subdomain index of the corresponding cell in `finite_cells`
+ *  \param verbose if `true`, prints information about the reading process
+ *
+ * \returns `true` if the reading was successful, `false` otherwise.
+ */
+
+template<class PointRange, class CellRange>
+bool read_MEDIT(std::istream& is,
+                PointRange& points,
+                CellRange& finite_cells,
+                std::vector<int>& subdomains,
+                bool verbose = false)
+
+{
+  boost::unordered_map<std::array<int,3>,int > border_facets;
+  constexpr bool read_border_facets = false;
+  bool is_CGAL_mesh;
+  return internal::read_MEDIT(is, points, finite_cells, subdomains, border_facets,
+                              read_border_facets, verbose, is_CGAL_mesh);
+}
+
+} // namespace IO
+
+} // namespace CGAL
+
+#endif // CGAL_IO_MEDIT_H
diff --git a/Stream_support/include/CGAL/IO/PLY.h b/Stream_support/include/CGAL/IO/PLY.h
index 4e7a365c001f..ad87a6d670a7 100644
--- a/Stream_support/include/CGAL/IO/PLY.h
+++ b/Stream_support/include/CGAL/IO/PLY.h
@@ -30,7 +30,6 @@
 #include <type_traits>
 #include <utility>
 #include <vector>
-#include <type_traits>
 
 namespace CGAL {