diff --git a/src/common/PointCloud.h b/src/common/PointCloud.h
index ed966d8..1b11a29 100644
--- a/src/common/PointCloud.h
+++ b/src/common/PointCloud.h
@@ -42,19 +42,19 @@ namespace pubgeo {
         int yOff;
         int zOff;
 
-        inline float x(int i) {
+        inline float x(int i) const {
             if (pv != nullptr)
                 return (float) (pv->getFieldAs<double>(pdal::Dimension::Id::X, i) - xOff);
             throw "Point set has not be initialized.";
         }
 
-        inline float y(int i) {
+        inline float y(int i) const {
             if (pv != nullptr)
                 return (float) (pv->getFieldAs<double>(pdal::Dimension::Id::Y, i) - yOff);
             throw "Point set has not be initialized.";
         }
 
-        inline float z(int i) {
+        inline float z(int i) const {
             if (pv != nullptr)
                 return (float) (pv->getFieldAs<double>(pdal::Dimension::Id::Z, i) - zOff);
             throw "Point set has not be initialized.";
diff --git a/src/common/orthoimage.h b/src/common/orthoimage.h
index ca0610e..4b30a53 100644
--- a/src/common/orthoimage.h
+++ b/src/common/orthoimage.h
@@ -286,8 +286,11 @@ namespace pubgeo {
             // Read a PSET file (e.g., BPF or LAS).
             PointCloud pset;
             bool ok = pset.Read(fileName);
-            if (!ok) return false;
+            return ok && readFromPointCloud(pset, gsdMeters, mode);
+        }
 
+        // Create image from point cloud.
+        bool readFromPointCloud(const PointCloud& pset, float gsdMeters, MIN_MAX_TYPE mode = MIN_VALUE) {
             // Calculate scale and offset for conversion to TYPE.
             float minVal = pset.bounds.zMin - 1;    // Reserve zero for noData value
             float maxVal = pset.bounds.zMax + 1;
diff --git a/src/shr3d/main.cpp b/src/shr3d/main.cpp
index f5e7475..5b562c1 100644
--- a/src/shr3d/main.cpp
+++ b/src/shr3d/main.cpp
@@ -103,7 +103,12 @@ int main(int argc, char **argv) {
 #endif
     } else if ((strcmp(ext, "las") == 0) || (strcmp(ext, "bpf") == 0)) {
         // First get the max Z values for the DSM.
-        bool ok = dsmImage.readFromPointCloud(inputFileName, (float) dh_meters, shr3d::MAX_VALUE);
+        // Read a PSET file (e.g., BPF or LAS).
+        shr3d::PointCloud pset;
+        bool ok = pset.Read(inputFileName);
+        if (!ok) return -1;
+
+        ok = dsmImage.readFromPointCloud(pset, (float) dh_meters, shr3d::MAX_VALUE);
         if (!ok) return -1;
 
         // Median filter, replacing only points differing by more than the AGL threshold.
@@ -118,7 +123,7 @@ int main(int argc, char **argv) {
         dsmImage.write(dsmOutFileName, true);
 
         // Now get the minimum Z values for the DTM.
-        ok = minImage.readFromPointCloud(inputFileName, (float) dh_meters, shr3d::MIN_VALUE);
+        ok = minImage.readFromPointCloud(pset, (float) dh_meters, shr3d::MIN_VALUE);
         if (!ok) return -1;
 
         // Median filter, replacing only points differing by more than the AGL threshold.