diff --git a/SEFramework/src/lib/Psf/VariablePsf.cpp b/SEFramework/src/lib/Psf/VariablePsf.cpp
index c8367c350..6e66e98a9 100644
--- a/SEFramework/src/lib/Psf/VariablePsf.cpp
+++ b/SEFramework/src/lib/Psf/VariablePsf.cpp
@@ -21,10 +21,14 @@
  *      Author: Alejandro Álvarez Ayllón
  */
 
-#include <ElementsKernel/Exception.h>
 #include <algorithm>
+#include <ElementsKernel/Logging.h>
+#include <ElementsKernel/Exception.h>
+#include "SEUtils/IsNan.h"
+
 #include "SEFramework/Psf/VariablePsf.h"
 
+static auto stack_logger = Elements::Logging::getLogger("PSFExPsf");
 
 namespace SourceXtractor {
 
@@ -130,6 +134,16 @@ void VariablePsf::selfTest() {
     if (coeff->getWidth() != psf_width || coeff->getHeight() != psf_height) {
       throw Elements::Exception() << "Malformed variable PSF, coefficient matrices do not have the same dimensions";
     }
+    for (auto x = 0; x < psf_width; ++x){
+      for (auto y = 0; y < psf_height; ++y) {
+        if (fastmath_isnan(coeff->at(x, y))) {
+          throw Elements::Exception() << "Malformed variable PSF, coefficient matrices contains NANs";
+        }
+        else if (fastmath_isinf(coeff->at(x, y))){
+          throw Elements::Exception() << "Malformed variable PSF, coefficient matrices contains INFs";
+        }
+      }
+    }
   }
 }
 
diff --git a/SEImplementation/src/lib/Plugin/Psf/PsfPluginConfig.cpp b/SEImplementation/src/lib/Plugin/Psf/PsfPluginConfig.cpp
index 462835f93..e0c61a9da 100644
--- a/SEImplementation/src/lib/Plugin/Psf/PsfPluginConfig.cpp
+++ b/SEImplementation/src/lib/Plugin/Psf/PsfPluginConfig.cpp
@@ -124,10 +124,13 @@ static std::shared_ptr<VariablePsf> readPsfEx(std::unique_ptr<CCfits::FITS> &pFi
     }
 
     return std::make_shared<VariablePsf>(pixel_sampling, components, group_degrees, coefficients);
-  } catch (CCfits::FITS::NoSuchHDU&) {
+  } catch (CCfits::FITS::NoSuchHDU&) { // Make sure we propagate this specific exception
     throw;
   } catch (CCfits::FitsException &e) {
     throw Elements::Exception() << "Error loading PSFEx file: " << e.message();
+  } catch (...) {
+    logger.error() << "Failed loading a psf file: " << pFits->name();
+    throw;
   }
 }
 
diff --git a/SEUtils/SEUtils/IsNan.h b/SEUtils/SEUtils/IsNan.h
new file mode 100644
index 000000000..c4dea6479
--- /dev/null
+++ b/SEUtils/SEUtils/IsNan.h
@@ -0,0 +1,80 @@
+#ifndef _SEUTILS_ISNAN_H
+#define _SEUTILS_ISNAN_H
+
+#include <cstdint>
+
+namespace SourceXtractor {
+
+/*
+ * From:
+ * https://github.com/searchivarius/BlogCode/tree/2d947d8c42ab0e10f1f28a2ad036eee10389997a/2017/6/1
+ *
+ * IEEE 754 compliant simple functions to test for NANs and INFs.
+ * This functions are necessary, b/c isnan doesn't work with -Ofast -ffast-math flags
+ * see http://searchivarius.org/blog/gcc-disables-isnan-and-isinf-when-compiling-ffast-math-flag  
+ *
+ * See also test files. The file ./regular is compiled without -ffast-math flag and it
+ * checks for a large number of values that the output of my functions is the same
+ * as the output of standard functions. For single precision numbers, i.e., for floats
+ * these checks are exhaustive. That is, I go over the set of all 4B+ possible float values.
+ * For doubles, this is not possible, so I use tests where the lower 32 bits of mantissa
+ * are set to zero.
+ * 
+ * Copyright (c) 2017   Leonid Boytsov
+ *
+ * This code is released under the
+ * Apache License Version 2.0 http://www.apache.org/licenses/.
+ *
+ */
+
+// A mask to extract an exponent from the single-precision floating point number
+// 01111111100000000000000000000000 
+const unsigned FLOAT_EXP_MASK = 0x7F800000; 
+// A mask to extract a mantissa/fractional part from the single-precision floating point number
+// 00000000011111111111111111111111 
+const unsigned FLOAT_FRAC_PART_MASK = 0x7FFFFF;
+
+inline bool fastmath_isnan(float x) {
+  union {
+    uint32_t u;
+    float    f;
+  } conv;
+  conv.f = x;
+  return ((conv.u & FLOAT_EXP_MASK) == FLOAT_EXP_MASK) && ((conv.u & FLOAT_FRAC_PART_MASK) != 0);
+};
+
+inline bool fastmath_isinf(float x) {
+  union {
+    uint32_t u;
+    float    f;
+  } conv;
+  conv.f = x;
+  return ((conv.u & FLOAT_EXP_MASK) == FLOAT_EXP_MASK) && ((conv.u & FLOAT_FRAC_PART_MASK) == 0);
+};
+
+// 0111 1111 1111 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000
+const uint64_t DOUBLE_EXP_MASK = 0x7FF0000000000000ul;
+// 0000 0000 0000 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111
+const uint64_t DOUBLE_FRAC_PART_MASK = 0x000FFFFFFFFFFFFFul;
+
+inline bool fastmath_isnan(double x) {
+  union {
+    uint64_t  u;
+    double    f;
+  } conv;
+  conv.f = x;
+  return ((conv.u & DOUBLE_EXP_MASK) == DOUBLE_EXP_MASK) && ((conv.u & DOUBLE_FRAC_PART_MASK) != 0);
+};
+
+inline bool fastmath_isinf(double x) {
+  union {
+    uint64_t  u;
+    double    f;
+  } conv;
+  conv.f = x;
+  return ((conv.u & DOUBLE_EXP_MASK) == DOUBLE_EXP_MASK) && ((conv.u & DOUBLE_FRAC_PART_MASK) == 0);
+};
+
+} // namespace SourceXtractor
+
+#endif