Merge pull request #341 from dirac-institute/prefactor

Prefactoring: Rename and move interfaces

README.md

@@ -100,7 +100,7 @@ ds = FakeDataSet(512, 512, 10)
 imgs = ds.stack.get_images()
 
 # Get the timestamp of the first image.
-t0 = imgs[0].get_time()
+t0 = imgs[0].get_obstime()
 print(f"Image times start at {t0}.")
 
 # Specify an artificial object
@@ -110,7 +110,7 @@ velocity = (2, 0)
 
 # Inject object into images
 for im in imgs:
-    dt = im.get_time() - t0
+    dt = im.get_obstime() - t0
     im.add_object(position[0] + dt * velocity[0], 
                   position[1] + dt * velocity[1], 
                   flux)

notebooks/Kbmod_Reference.ipynb

@@ -167,7 +167,7 @@
    "outputs": [],
    "source": [
     "im = kb.layered_image(im_path + \"000000.fits\", p)\n",
-    "print(f\"Loaded a {im.get_width()} by {im.get_height()} image at time {im.get_time()}\")"
+    "print(f\"Loaded a {im.get_width()} by {im.get_height()} image at time {im.get_obstime()}\")"
    ]
   },
   {
@@ -204,8 +204,8 @@
    "source": [
     "print(f\"Width = {im.get_width()}\")\n",
     "print(f\"Height = {im.get_height()}\")\n",
-    "print(f\"Pixels Per Image = {im.get_ppi()}\")\n",
-    "print(f\"Time = {im.get_time()}\")"
+    "print(f\"Pixels Per Image = {im.get_npixels()}\")\n",
+    "print(f\"Time = {im.get_obstime()}\")"
    ]
   },
   {
@@ -461,7 +461,7 @@
     "# where the stack contains images of different sizes.\n",
     "print(f\"Width = {stack.get_width()}\")\n",
     "print(f\"Height = {stack.get_height()}\")\n",
-    "print(f\"Pixels Per Image = {stack.get_ppi()}\")\n",
+    "print(f\"Pixels Per Image = {stack.get_npixels()}\")\n",
     "\n",
     "# Retrieve a list of layered_images back from the stack.\n",
     "stack.get_images()\n",

notebooks/create_fake_data.ipynb

@@ -111,9 +111,9 @@
    ],
    "source": [
     "# Check the object was inserted correctly.\n",
-    "t0 = ds.stack.get_single_image(0).get_time()\n",
+    "t0 = ds.stack.get_single_image(0).get_obstime()\n",
     "for i in range(ds.stack.img_count()):\n",
-    "    ti = ds.stack.get_single_image(i).get_time()\n",
+    "    ti = ds.stack.get_single_image(i).get_obstime()\n",
     "    dt = ti - t0\n",
     "    px = int(trj.x + dt * trj.x_v + 0.5)\n",
     "    py = int(trj.y + dt * trj.y_v + 0.5)\n",

notebooks/kbmod_visualize.ipynb

@@ -62,7 +62,7 @@
    "source": [
     "p = kb.psf(1.0)\n",
     "im = kb.layered_image(im_path + \"000000.fits\", p)\n",
-    "print(f\"Loaded a {im.get_width()} by {im.get_height()} image at time {im.get_time()}\")"
+    "print(f\"Loaded a {im.get_width()} by {im.get_height()} image at time {im.get_obstime()}\")"
    ]
   },
   {

src/kbmod/analysis_utils.py

@@ -130,7 +130,7 @@ def load_images(
             if verbose:
                 print(f"Loading file: {full_file_path}")
             img = kb.layered_image(full_file_path, psf)
-            img.set_time(time_stamp)
+            img.set_obstime(time_stamp)
 
             # Save the file, time, and image information.
             img_info.append(header_info)

src/kbmod/image_info.py

@@ -140,7 +140,7 @@ def set_epoch(self, epoch):
         """
         self.epoch_ = Time(epoch)
         if self.image is not None:
-            self.image.set_time(self.epoch_.mjd)
+            self.image.set_obstime(self.epoch_.mjd)
         self.epoch_set_ = True
 
     def get_epoch(self, none_if_unset=False):

src/kbmod/search/ImageStack.cpp

@@ -50,7 +50,7 @@ void ImageStack::extractImageTimes() {
     // Load image times
     imageTimes = std::vector<float>();
     for (auto& i : images) {
-        imageTimes.push_back(float(i.getTime()));
+        imageTimes.push_back(float(i.getObstime()));
     }
 }
 
@@ -85,7 +85,7 @@ void ImageStack::convolvePSF() {
     for (auto& i : images) i.convolvePSF();
 }
 
-void ImageStack::saveGlobalMask(const std::string& path) { globalMask.saveToFile(path, false); }
+void ImageStack::saveGlobalMask(const std::string& path) { globalMask.saveToFile(path); }
 
 void ImageStack::saveImages(const std::string& path) {
     for (auto& i : images) i.saveLayers(path);
@@ -115,21 +115,21 @@ void ImageStack::growMask(int steps) {
 }
 
 void ImageStack::createGlobalMask(int flags, int threshold) {
-    int ppi = getPPI();
+    int npixels = getNPixels();
 
     // For each pixel count the number of images where it is masked.
-    std::vector<int> counts(ppi, 0);
+    std::vector<int> counts(npixels, 0);
     for (unsigned int img = 0; img < images.size(); ++img) {
         float* imgMask = images[img].getMDataRef();
         // Count the number of times a pixel has any of the flags
-        for (unsigned int pixel = 0; pixel < ppi; ++pixel) {
+        for (unsigned int pixel = 0; pixel < npixels; ++pixel) {
             if ((flags & static_cast<int>(imgMask[pixel])) != 0) counts[pixel]++;
         }
     }
 
     // Set all pixels below threshold to 0 and all above to 1
     float* globalM = globalMask.getDataRef();
-    for (unsigned int p = 0; p < ppi; ++p) {
+    for (unsigned int p = 0; p < npixels; ++p) {
         globalM[p] = counts[p] < threshold ? 0.0 : 1.0;
     }
 }

src/kbmod/search/ImageStack.h

@@ -28,7 +28,7 @@ class ImageStack {
     unsigned imgCount() const { return images.size(); }
     unsigned getWidth() const { return images.size() > 0 ? images[0].getWidth() : 0; }
     unsigned getHeight() const { return images.size() > 0 ? images[0].getHeight() : 0; }
-    unsigned getPPI() const { return images.size() > 0 ? images[0].getPPI() : 0; }
+    unsigned getNPixels() const { return images.size() > 0 ? images[0].getNPixels() : 0; }
     std::vector<LayeredImage>& getImages() { return images; }
     const std::vector<float>& getTimes() const { return imageTimes; }
     float* getTimesDataRef() { return imageTimes.data(); }

src/kbmod/search/KBMOSearch.cpp

@@ -223,8 +223,8 @@ void KBMOSearch::saveImages(const std::string& path) {
         std::string number = std::to_string(i);
         // Add leading zeros
         number = std::string(4 - number.length(), '0') + number;
-        psiImages[i].saveToFile(path + "/psi/PSI" + number + ".fits", false);
-        phiImages[i].saveToFile(path + "/phi/PHI" + number + ".fits", false);
+        psiImages[i].saveToFile(path + "/psi/PSI" + number + ".fits");
+        phiImages[i].saveToFile(path + "/phi/PHI" + number + ".fits");
     }
 }
 
@@ -262,10 +262,10 @@ void KBMOSearch::fillPsiAndPhiVects(const std::vector<RawImage>& psiImgs,
     assert(num_images > 0);
     assert(phiImgs.size() == num_images);
 
-    int num_pixels = psiImgs[0].getPPI();
+    int num_pixels = psiImgs[0].getNPixels();
     for (int i = 0; i < num_images; ++i) {
-        assert(psiImgs[i].getPPI() == num_pixels);
-        assert(phiImgs[i].getPPI() == num_pixels);
+        assert(psiImgs[i].getNPixels() == num_pixels);
+        assert(phiImgs[i].getNPixels() == num_pixels);
     }
 
     psiVect->clear();

src/kbmod/search/LayeredImage.cpp

@@ -15,14 +15,25 @@ LayeredImage::LayeredImage(std::string path, const PointSpreadFunc& psf) : psf(p
     int fBegin = path.find_last_of("/");
     int fEnd = path.find_last_of(".fits") - 4;
     fileName = path.substr(fBegin, fEnd - fBegin);
-    readHeader(path);
-    science = RawImage(width, height);
-    mask = RawImage(width, height);
-    variance = RawImage(width, height);
-    loadLayers(path);
+
+    science = RawImage();
+    science.loadFromFile(path, 1);
+    width = science.getWidth();
+    height = science.getHeight();
+
+    mask = RawImage();
+    mask.loadFromFile(path, 2);
+
+    variance = RawImage();
+    variance.loadFromFile(path, 3);
+
+    if (width != variance.getWidth() or height != variance.getHeight())
+        throw std::runtime_error("Science and Variance layers are not the same size.");
+    if (width != mask.getWidth() or height != mask.getHeight())
+        throw std::runtime_error("Science and Mask layers are not the same size.");
 }
 
-LayeredImage::LayeredImage(const RawImage& sci, const RawImage& var, const RawImage& msk, float time,
+LayeredImage::LayeredImage(const RawImage& sci, const RawImage& var, const RawImage& msk,
                            const PointSpreadFunc& psf)
         : psf(psf), psfSQ(psf) {
     // Get the dimensions of the science layer and check for consistency with
@@ -35,7 +46,6 @@ LayeredImage::LayeredImage(const RawImage& sci, const RawImage& var, const RawIm
         throw std::runtime_error("Science and Mask layers are not the same size.");
 
     // Set the remaining variables.
-    captureTime = time;
     psfSQ.squarePSF();
 
     // Copy the image layers.
@@ -54,7 +64,6 @@ LayeredImage::LayeredImage(std::string name, int w, int h, float noiseStDev, flo
     fileName = name;
     width = w;
     height = h;
-    captureTime = time;
     psfSQ.squarePSF();
 
     std::vector<float> rawSci(width * height);
@@ -65,62 +74,14 @@ LayeredImage::LayeredImage(std::string name, int w, int h, float noiseStDev, flo
     }
     std::normal_distribution<float> distrib(0.0, noiseStDev);
     for (float& p : rawSci) p = distrib(generator);
+
     science = RawImage(w, h, rawSci);
+    science.setObstime(time);
+
     mask = RawImage(w, h, std::vector<float>(w * h, 0.0));
     variance = RawImage(w, h, std::vector<float>(w * h, pixelVariance));
 }
 
-/* Read the image dimensions and capture time from header */
-void LayeredImage::readHeader(const std::string& filePath) {
-    fitsfile* fptr;
-    int status = 0;
-    int mjdStatus = 0;
-    int fileNotFound;
-
-    // Open header to read MJD
-    if (fits_open_file(&fptr, filePath.c_str(), READONLY, &status))
-        throw std::runtime_error("Could not open file");
-
-    // Read image capture time, ignore error if does not exist
-    captureTime = 0.0;
-    fits_read_key(fptr, TDOUBLE, "MJD", &captureTime, NULL, &mjdStatus);
-
-    if (fits_close_file(fptr, &status)) fits_report_error(stderr, status);
-
-    // Reopen header for first layer to get image dimensions
-    if (fits_open_file(&fptr, (filePath + "[1]").c_str(), READONLY, &status))
-        fits_report_error(stderr, status);
-
-    // Read image Dimensions
-    long dimensions[2];
-    if (fits_read_keys_lng(fptr, "NAXIS", 1, 2, dimensions, &fileNotFound, &status))
-        fits_report_error(stderr, status);
-
-    width = dimensions[0];
-    height = dimensions[1];
-
-    if (fits_close_file(fptr, &status)) fits_report_error(stderr, status);
-}
-
-void LayeredImage::loadLayers(const std::string& filePath) {
-    // Load images from file into layers' pixels
-    readFitsImg((filePath + "[1]").c_str(), science.getDataRef());
-    readFitsImg((filePath + "[2]").c_str(), mask.getDataRef());
-    readFitsImg((filePath + "[3]").c_str(), variance.getDataRef());
-}
-
-void LayeredImage::readFitsImg(const char* name, float* target) {
-    fitsfile* fptr;
-    int nullval = 0;
-    int anynull;
-    int status = 0;
-
-    if (fits_open_file(&fptr, name, READONLY, &status)) fits_report_error(stderr, status);
-    if (fits_read_img(fptr, TFLOAT, 1, getPPI(), &nullval, target, &anynull, &status))
-        fits_report_error(stderr, status);
-    if (fits_close_file(fptr, &status)) fits_report_error(stderr, status);
-}
-
 void LayeredImage::setPSF(const PointSpreadFunc& new_psf) {
     psf = new_psf;
     psfSQ = new_psf;
@@ -165,7 +126,7 @@ void LayeredImage::applyGlobalMask(const RawImage& globalM) {
 }
 
 void LayeredImage::applyMaskThreshold(float thresh) {
-    const int numPixels = getPPI();
+    const int numPixels = getNPixels();
     float* sciPix = science.getDataRef();
     float* varPix = variance.getDataRef();
     for (int i = 0; i < numPixels; ++i) {
@@ -178,13 +139,13 @@ void LayeredImage::applyMaskThreshold(float thresh) {
 
 void LayeredImage::subtractTemplate(const RawImage& subTemplate) {
     assert(getHeight() == subTemplate.getHeight() && getWidth() == subTemplate.getWidth());
-    const int numPixels = getPPI();
+    const int numPixels = getNPixels();
 
     float* sciPix = science.getDataRef();
-    const std::vector<float>& tempPix = subTemplate.getPixels();
+    const std::vector<float>& temNPixelsx = subTemplate.getPixels();
     for (unsigned i = 0; i < numPixels; ++i) {
-        if ((sciPix[i] != NO_DATA) && (tempPix[i] != NO_DATA)) {
-            sciPix[i] -= tempPix[i];
+        if ((sciPix[i] != NO_DATA) && (temNPixelsx[i] != NO_DATA)) {
+            sciPix[i] -= temNPixelsx[i];
         }
     }
 }
@@ -193,6 +154,8 @@ void LayeredImage::saveLayers(const std::string& path) {
     fitsfile* fptr;
     int status = 0;
     long naxes[2] = {0, 0};
+    double obstime = science.getObstime();
+
     fits_create_file(&fptr, (path + fileName + ".fits").c_str(), &status);
 
     // If we are unable to create the file, check if it already exists
@@ -207,25 +170,13 @@ void LayeredImage::saveLayers(const std::string& path) {
     }
 
     fits_create_img(fptr, SHORT_IMG, 0, naxes, &status);
-    fits_update_key(fptr, TDOUBLE, "MJD", &captureTime, "[d] Generated Image time", &status);
+    fits_update_key(fptr, TDOUBLE, "MJD", &obstime, "[d] Generated Image time", &status);
     fits_close_file(fptr, &status);
     fits_report_error(stderr, status);
 
-    science.saveToFile(path + fileName + ".fits", true);
-    mask.saveToFile(path + fileName + ".fits", true);
-    variance.saveToFile(path + fileName + ".fits", true);
-}
-
-void LayeredImage::saveSci(const std::string& path) {
-    science.saveToFile(path + fileName + "SCI.fits", false);
-}
-
-void LayeredImage::saveMask(const std::string& path) {
-    mask.saveToFile(path + fileName + "MASK.fits", false);
-}
-
-void LayeredImage::saveVar(const std::string& path) {
-    variance.saveToFile(path + fileName + "VAR.fits", false);
+    science.appendLayerToFile(path + fileName + ".fits");
+    mask.appendLayerToFile(path + fileName + ".fits");
+    variance.appendLayerToFile(path + fileName + ".fits");
 }
 
 void LayeredImage::setScience(RawImage& im) {
@@ -255,7 +206,7 @@ RawImage LayeredImage::generatePsiImage() {
     float* varArray = getVDataRef();
 
     // Set each of the result pixels.
-    const int num_pixels = getPPI();
+    const int num_pixels = getNPixels();
     for (int p = 0; p < num_pixels; ++p) {
         float varPix = varArray[p];
         if (varPix != NO_DATA) {
@@ -277,7 +228,7 @@ RawImage LayeredImage::generatePhiImage() {
     float* varArray = getVDataRef();
 
     // Set each of the result pixels.
-    const int num_pixels = getPPI();
+    const int num_pixels = getNPixels();
     for (int p = 0; p < num_pixels; ++p) {
         float varPix = varArray[p];
         if (varPix != NO_DATA) {