diff --git a/README.md b/README.md index 4fa9d224..33031a74 100644 --- a/README.md +++ b/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) diff --git a/notebooks/Kbmod_Reference.ipynb b/notebooks/Kbmod_Reference.ipynb index a64656a6..3edf933e 100644 --- a/notebooks/Kbmod_Reference.ipynb +++ b/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", diff --git a/notebooks/create_fake_data.ipynb b/notebooks/create_fake_data.ipynb index a11b0867..824ee1f0 100644 --- a/notebooks/create_fake_data.ipynb +++ b/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", diff --git a/notebooks/kbmod_visualize.ipynb b/notebooks/kbmod_visualize.ipynb index 3b97184a..7c688119 100644 --- a/notebooks/kbmod_visualize.ipynb +++ b/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()}\")" ] }, { diff --git a/src/kbmod/analysis_utils.py b/src/kbmod/analysis_utils.py index 32825e28..79582482 100644 --- a/src/kbmod/analysis_utils.py +++ b/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) diff --git a/src/kbmod/image_info.py b/src/kbmod/image_info.py index 65e05eaa..22f78f77 100644 --- a/src/kbmod/image_info.py +++ b/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): diff --git a/src/kbmod/search/ImageStack.cpp b/src/kbmod/search/ImageStack.cpp index b1ec75f1..e8f9d187 100644 --- a/src/kbmod/search/ImageStack.cpp +++ b/src/kbmod/search/ImageStack.cpp @@ -50,7 +50,7 @@ void ImageStack::extractImageTimes() { // Load image times imageTimes = std::vector(); 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 counts(ppi, 0); + std::vector 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(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; } } diff --git a/src/kbmod/search/ImageStack.h b/src/kbmod/search/ImageStack.h index 88b8fa93..f6fb1324 100644 --- a/src/kbmod/search/ImageStack.h +++ b/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& getImages() { return images; } const std::vector& getTimes() const { return imageTimes; } float* getTimesDataRef() { return imageTimes.data(); } diff --git a/src/kbmod/search/KBMOSearch.cpp b/src/kbmod/search/KBMOSearch.cpp index 281cd693..749b694e 100644 --- a/src/kbmod/search/KBMOSearch.cpp +++ b/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& 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(); diff --git a/src/kbmod/search/LayeredImage.cpp b/src/kbmod/search/LayeredImage.cpp index a4cf7c35..337758d0 100644 --- a/src/kbmod/search/LayeredImage.cpp +++ b/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 rawSci(width * height); @@ -65,62 +74,14 @@ LayeredImage::LayeredImage(std::string name, int w, int h, float noiseStDev, flo } std::normal_distribution 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(w * h, 0.0)); variance = RawImage(w, h, std::vector(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& tempPix = subTemplate.getPixels(); + const std::vector& 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) { diff --git a/src/kbmod/search/LayeredImage.h b/src/kbmod/search/LayeredImage.h index d004d9a5..936782d9 100644 --- a/src/kbmod/search/LayeredImage.h +++ b/src/kbmod/search/LayeredImage.h @@ -26,7 +26,7 @@ namespace search { class LayeredImage { public: explicit LayeredImage(std::string path, const PointSpreadFunc& psf); - explicit LayeredImage(const RawImage& sci, const RawImage& var, const RawImage& msk, float time, + explicit LayeredImage(const RawImage& sci, const RawImage& var, const RawImage& msk, const PointSpreadFunc& psf); explicit LayeredImage(std::string name, int w, int h, float noiseStDev, float pixelVariance, double time, const PointSpreadFunc& psf); @@ -42,11 +42,9 @@ class LayeredImage { std::string getName() const { return fileName; } unsigned getWidth() const { return width; } unsigned getHeight() const { return height; } - unsigned getPPI() const { return width * height; } - double getTime() const { return captureTime; } - - // Basic setter functions. - void setTime(double timestamp) { captureTime = timestamp; } + unsigned getNPixels() const { return width * height; } + double getObstime() const { return science.getObstime(); } + void setObstime(double obstime) { science.setObstime(obstime); } // Getter functions for the data in the individual layers. RawImage& getScience() { return science; } @@ -72,9 +70,6 @@ class LayeredImage { // Saves the data in each later to a file. void saveLayers(const std::string& path); - void saveSci(const std::string& path); - void saveMask(const std::string& path); - void saveVar(const std::string& path); // Setter functions for the individual layers. void setScience(RawImage& im); @@ -89,15 +84,11 @@ class LayeredImage { RawImage generatePhiImage(); private: - void readHeader(const std::string& filePath); - void loadLayers(const std::string& filePath); - void readFitsImg(const char* name, float* target); void checkDims(RawImage& im); std::string fileName; unsigned width; unsigned height; - double captureTime; PointSpreadFunc psf; PointSpreadFunc psfSQ; diff --git a/src/kbmod/search/RawImage.cpp b/src/kbmod/search/RawImage.cpp index b794d552..63305464 100644 --- a/src/kbmod/search/RawImage.cpp +++ b/src/kbmod/search/RawImage.cpp @@ -16,13 +16,14 @@ namespace search { int psfSize, int psfDim, int psfRadius, float psfSum); #endif -RawImage::RawImage() : width(0), height(0) { pixels = std::vector(); } +RawImage::RawImage() : width(0), height(0), obstime(-1.0) { pixels = std::vector(); } // Copy constructor RawImage::RawImage(const RawImage& old) { width = old.getWidth(); height = old.getHeight(); pixels = old.getPixels(); + obstime = old.getObstime(); } // Copy assignment @@ -30,12 +31,13 @@ RawImage& RawImage::operator=(const RawImage& source) { width = source.width; height = source.height; pixels = source.pixels; + obstime = source.obstime; return *this; } // Move constructor RawImage::RawImage(RawImage&& source) - : width(source.width), height(source.height), pixels(std::move(source.pixels)) {} + : width(source.width), height(source.height), obstime(source.obstime), pixels(std::move(source.pixels)) {} // Move assignment RawImage& RawImage::operator=(RawImage&& source) { @@ -43,18 +45,22 @@ RawImage& RawImage::operator=(RawImage&& source) { width = source.width; height = source.height; pixels = std::move(source.pixels); + obstime = source.obstime; } return *this; } -RawImage::RawImage(unsigned w, unsigned h) : height(h), width(w), pixels(w * h) {} +RawImage::RawImage(unsigned w, unsigned h) : height(h), width(w), obstime(-1.0), pixels(w * h) {} -RawImage::RawImage(unsigned w, unsigned h, const std::vector& pix) : width(w), height(h), pixels(pix) { +RawImage::RawImage(unsigned w, unsigned h, const std::vector& pix) : width(w), height(h), obstime(-1.0), pixels(pix) { assert(w * h == pix.size()); } #ifdef Py_PYTHON_H -RawImage::RawImage(pybind11::array_t arr) { setArray(arr); } +RawImage::RawImage(pybind11::array_t arr) { + obstime = -1.0; + setArray(arr); +} void RawImage::setArray(pybind11::array_t& arr) { pybind11::buffer_info info = arr.request(); @@ -65,7 +71,7 @@ void RawImage::setArray(pybind11::array_t& arr) { height = info.shape[0]; float* pix = static_cast(info.ptr); - pixels = std::vector(pix, pix + getPPI()); + pixels = std::vector(pix, pix + getNPixels()); } #endif @@ -92,15 +98,94 @@ bool RawImage::approxEqual(const RawImage& imgB, float atol) const { return true; } -void RawImage::saveToFile(const std::string& path, bool append) { +// Load the image data from a specific layer of a FITS file. +void RawImage::loadFromFile(const std::string& filePath, int layer_num) { + // Open the file's header and read in the obstime and the dimensions. + fitsfile* fptr; + int status = 0; + int mjdStatus = 0; + int fileNotFound; + int nullval = 0; + int anynull = 0; + + // Open the correct layer to extract the RawImage. + std::string layerPath = filePath + "[" + std::to_string(layer_num) + "]"; + if (fits_open_file(&fptr, layerPath.c_str(), READONLY, &status)) { + fits_report_error(stderr, status); + throw std::runtime_error("Could not open FITS file to read RawImage"); + } + + // 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]; + + // Read in the image. + pixels = std::vector(width * height); + if (fits_read_img(fptr, TFLOAT, 1, getNPixels(), &nullval, pixels.data(), &anynull, &status)) + fits_report_error(stderr, status); + + // Read image observation time, ignore error if does not exist + obstime = -1.0; + if (fits_read_key(fptr, TDOUBLE, "MJD", &obstime, NULL, &mjdStatus)) obstime = -1.0; + if (fits_close_file(fptr, &status)) fits_report_error(stderr, status); + + // If we are reading from a sublayer and did not find a time, try the overall header. + if (obstime < 0.0) { + if (fits_open_file(&fptr, filePath.c_str(), READONLY, &status)) + throw std::runtime_error("Could not open FITS file to read RawImage"); + fits_read_key(fptr, TDOUBLE, "MJD", &obstime, NULL, &mjdStatus); + if (fits_close_file(fptr, &status)) fits_report_error(stderr, status); + } +} + +void RawImage::saveToFile(const std::string& filename) { + fitsfile* fptr; + int status = 0; + long naxes[2] = {0, 0}; + + fits_create_file(&fptr, filename.c_str(), &status); + + // If we are unable to create the file, check if it already exists + // and, if so, delete it and retry the create. + if (status == 105) { + status = 0; + fits_open_file(&fptr, filename.c_str(), READWRITE, &status); + if (status == 0) { + fits_delete_file(fptr, &status); + fits_create_file(&fptr, filename.c_str(), &status); + } + } + + // Create the primary array image (32-bit float pixels) + long dimensions[2]; + dimensions[0] = width; + dimensions[1] = height; + fits_create_img(fptr, FLOAT_IMG, 2 /*naxis*/, dimensions, &status); + fits_report_error(stderr, status); + + /* Write the array of floats to the image */ + fits_write_img(fptr, TFLOAT, 1, getNPixels(), pixels.data(), &status); + fits_report_error(stderr, status); + + // Add the basic header data. + fits_update_key(fptr, TDOUBLE, "MJD", &obstime, "[d] Generated Image time", &status); + fits_report_error(stderr, status); + + fits_close_file(fptr, &status); + fits_report_error(stderr, status); +} + +void RawImage::appendLayerToFile(const std::string& filename) { int status = 0; fitsfile* f; - // Create a new file if append is false or we cannot open - // the specified file. - if (!append || fits_open_file(&f, path.c_str(), READWRITE, &status)) { - fits_create_file(&f, (path).c_str(), &status); + // Check that we can open the file. + if (fits_open_file(&f, filename.c_str(), READWRITE, &status)) { fits_report_error(stderr, status); + throw std::runtime_error("Unable to open FITS file for appending."); } // This appends a layer (extension) if the file exists) @@ -112,8 +197,13 @@ void RawImage::saveToFile(const std::string& path, bool append) { fits_report_error(stderr, status); /* Write the array of floats to the image */ - fits_write_img(f, TFLOAT, 1, getPPI(), pixels.data(), &status); + fits_write_img(f, TFLOAT, 1, getNPixels(), pixels.data(), &status); fits_report_error(stderr, status); + + // Save the image time in the header. + fits_update_key(f, TDOUBLE, "MJD", &obstime, "[d] Generated Image time", &status); + fits_report_error(stderr, status); + fits_close_file(f, &status); fits_report_error(stderr, status); } @@ -135,6 +225,8 @@ RawImage RawImage::createStamp(float x, float y, int radius, bool interpolate, b stamp.setPixel(xoff, yoff, pixVal); } } + + stamp.setObstime(obstime); return stamp; } @@ -170,8 +262,8 @@ void RawImage::convolve_cpu(const PointSpreadFunc& psf) { } // Copy the data into the pixels vector. - const int ppi = width * height; - for(int i = 0; i < ppi; ++i) { + const int npixels = getNPixels(); + for(int i = 0; i < npixels; ++i) { pixels[i] = result[i]; } } @@ -187,8 +279,8 @@ void RawImage::convolve(PointSpreadFunc psf) { void RawImage::applyMask(int flags, const std::vector& exceptions, const RawImage& mask) { const std::vector& maskPix = mask.getPixels(); - const int num_pixels = getPPI(); - assert(num_pixels == mask.getPPI()); + const int num_pixels = getNPixels(); + assert(num_pixels == mask.getNPixels()); for (unsigned int p = 0; p < num_pixels; ++p) { int pixFlags = static_cast(maskPix[p]); bool isException = false; @@ -337,7 +429,7 @@ void RawImage::setAllPix(float value) { } std::array RawImage::computeBounds() const { - const int num_pixels = getPPI(); + const int num_pixels = getNPixels(); float minVal = FLT_MAX; float maxVal = -FLT_MAX; for (unsigned p = 0; p < num_pixels; ++p) { diff --git a/src/kbmod/search/RawImage.h b/src/kbmod/search/RawImage.h index b06e6660..2cd4c745 100644 --- a/src/kbmod/search/RawImage.h +++ b/src/kbmod/search/RawImage.h @@ -36,6 +36,7 @@ class RawImage { RawImage(RawImage&& source); // Move constructor explicit RawImage(unsigned w, unsigned h); explicit RawImage(unsigned w, unsigned h, const std::vector& pix); + #ifdef Py_PYTHON_H explicit RawImage(pybind11::array_t arr); void setArray(pybind11::array_t& arr); @@ -47,7 +48,7 @@ class RawImage { // Basic getter functions for image data. unsigned getWidth() const { return width; } unsigned getHeight() const { return height; } - unsigned getPPI() const { return width * height; } + unsigned getNPixels() const { return width * height; } // Inline pixel functions. float getPixel(int x, int y) const { @@ -71,6 +72,10 @@ class RawImage { // Check if two raw images are approximately equal. bool approxEqual(const RawImage& imgB, float atol) const; + // Functions for locally storing the image time. + double getObstime() const { return obstime; } + void setObstime(double new_time) { obstime = new_time; } + // Compute the min and max bounds of values in the image. std::array computeBounds() const; @@ -89,9 +94,13 @@ class RawImage { // Grow the area of masked pixels. void growMask(int steps); - // Save the RawImage to a file. Append indicates whether to append - // or create a new file. - void saveToFile(const std::string& path, bool append); + // Load the image data from a specific layer of a FITS file. + // Overwrites the current image data. + void loadFromFile(const std::string& filePath, int layer_num); + + // Save the RawImage to a file (single layer) or append the layer to an existing file. + void saveToFile(const std::string& filename); + void appendLayerToFile(const std::string& filename); // Convolve the image with a point spread function. void convolve(PointSpreadFunc psf); @@ -119,6 +128,7 @@ class RawImage { unsigned width; unsigned height; std::vector pixels; + double obstime; }; // Helper functions for creating composite images. diff --git a/src/kbmod/search/bindings.cpp b/src/kbmod/search/bindings.cpp index 0a80d01b..bf620a4a 100644 --- a/src/kbmod/search/bindings.cpp +++ b/src/kbmod/search/bindings.cpp @@ -88,9 +88,11 @@ PYBIND11_MODULE(search, m) { .def(py::init>()) .def("get_height", &ri::getHeight, "Returns the image's height in pixels.") .def("get_width", &ri::getWidth, "Returns the image's width in pixels.") - .def("get_ppi", &ri::getPPI, "Returns the image's total number of pixels.") + .def("get_npixels", &ri::getNPixels, "Returns the image's total number of pixels.") .def("get_all_pixels", &ri::getPixels, "Returns a list of the images pixels.") .def("set_array", &ri::setArray, "Sets all image pixels given an array of values.") + .def("get_obstime", &ri::getObstime, "Get the observation time of the image.") + .def("set_obstime", &ri::setObstime, "Set the observation time of the image.") .def("approx_equal", &ri::approxEqual, "Checks if two images are approximately equal.") .def("compute_bounds", &ri::computeBounds, "Returns min and max pixel values.") .def("find_peak", &ri::findPeak, "Returns the pixel coordinates of the maximum value.") @@ -107,13 +109,15 @@ PYBIND11_MODULE(search, m) { .def("get_pixel_interp", &ri::getPixelInterp, "Get the interoplated value of a pixel.") .def("convolve", &ri::convolve, "Convolve the image with a PSF.") .def("convolve_cpu", &ri::convolve_cpu, "Convolve the image with a PSF.") - .def("save_fits", &ri::saveToFile, "Save the image to a FITS file."); + .def("load_fits", &ri::loadFromFile, "Load the image data from a FITS file.") + .def("save_fits", &ri::saveToFile, "Save the image to a FITS file.") + .def("append_fits_layer", &ri::appendLayerToFile, "Append the image as a layer in a FITS file."); m.def("create_median_image", &search::createMedianImage); m.def("create_summed_image", &search::createSummedImage); m.def("create_mean_image", &search::createMeanImage); py::class_
  • (m, "layered_image") .def(py::init()) - .def(py::init(), R"pbdoc( + .def(py::init(), R"pbdoc( Creates a layered_image out of individual `raw_image` layers. Parameters @@ -124,8 +128,6 @@ PYBIND11_MODULE(search, m) { The `raw_image` for the cariance layer. msk : `raw_image` The `raw_image` for the mask layer. - time : `float` - The image time stamp. p : `psf` The PSF for the image. @@ -142,9 +144,6 @@ PYBIND11_MODULE(search, m) { .def("apply_mask_threshold", &li::applyMaskThreshold) .def("sub_template", &li::subtractTemplate) .def("save_layers", &li::saveLayers) - .def("save_sci", &li::saveSci) - .def("save_mask", &li::saveMask) - .def("save_var", &li::saveVar) .def("get_science", &li::getScience, "Returns the science layer raw_image.") .def("get_mask", &li::getMask, "Returns the mask layer raw_image.") .def("get_variance", &li::getVariance, "Returns the variance layer raw_image.") @@ -157,9 +156,9 @@ PYBIND11_MODULE(search, m) { .def("get_name", &li::getName, "Returns the name of the layered image.") .def("get_width", &li::getWidth, "Returns the image's width in pixels.") .def("get_height", &li::getHeight, "Returns the image's height in pixels.") - .def("get_ppi", &li::getPPI, "Returns the image's total number of pixels.") - .def("get_time", &li::getTime) - .def("set_time", &li::setTime) + .def("get_npixels", &li::getNPixels, "Returns the image's total number of pixels.") + .def("get_obstime", &li::getObstime, "Get the image's observation time.") + .def("set_obstime", &li::setObstime, "Set the image's observation time.") .def("generate_psi_image", &li::generatePsiImage) .def("generate_phi_image", &li::generatePhiImage); py::class_(m, "image_stack") @@ -181,7 +180,7 @@ PYBIND11_MODULE(search, m) { .def("convolve_psf", &is::convolvePSF) .def("get_width", &is::getWidth) .def("get_height", &is::getHeight) - .def("get_ppi", &is::getPPI); + .def("get_npixels", &is::getNPixels); py::class_(m, "stack_search") .def(py::init()) .def("save_psi_phi", &ks::savePsiPhi) diff --git a/tests/test_analysis_utils.py b/tests/test_analysis_utils.py index f66f6f75..fe1d8a37 100644 --- a/tests/test_analysis_utils.py +++ b/tests/test_analysis_utils.py @@ -387,7 +387,7 @@ def test_file_load_basic(self): img = stack.get_single_image(i) self.assertEqual(img.get_width(), 64) self.assertEqual(img.get_height(), 64) - self.assertAlmostEqual(img.get_time(), true_times[i], delta=0.005) + self.assertAlmostEqual(img.get_obstime(), true_times[i], delta=0.005) self.assertAlmostEqual(1.0, img.get_psf().get_stdev()) # Check that visit IDs and times were extracted for each file in img_info. @@ -420,7 +420,7 @@ def test_file_load_extra(self): img = stack.get_single_image(i) self.assertEqual(img.get_width(), 64) self.assertEqual(img.get_height(), 64) - self.assertAlmostEqual(img.get_time(), true_times[i], delta=0.005) + self.assertAlmostEqual(img.get_obstime(), true_times[i], delta=0.005) self.assertAlmostEqual(psfs_std[i], img.get_psf().get_stdev()) # Check that visit IDs and times were extracted for each file in img_info. diff --git a/tests/test_bilinear_interp.py b/tests/test_bilinear_interp.py index 43f70abd..ed0b071b 100644 --- a/tests/test_bilinear_interp.py +++ b/tests/test_bilinear_interp.py @@ -37,7 +37,7 @@ def test_pixel_interp(self): im = kb.raw_image(pixels) self.assertEqual(im.get_width(), 3) self.assertEqual(im.get_height(), 2) - self.assertEqual(im.get_ppi(), 6) + self.assertEqual(im.get_npixels(), 6) # The middle of a pixel should interp to the pixel's value. self.assertAlmostEqual(im.get_pixel_interp(0.5, 0.5), 0.0, delta=0.001) diff --git a/tests/test_fake_data_creator.py b/tests/test_fake_data_creator.py index b6dbacdc..319a4af0 100644 --- a/tests/test_fake_data_creator.py +++ b/tests/test_fake_data_creator.py @@ -17,7 +17,7 @@ def test_create(self): self.assertEqual(layered.get_width(), 128) self.assertEqual(layered.get_height(), 256) - t = layered.get_time() + t = layered.get_obstime() self.assertGreater(t, last_time) last_time = t @@ -31,9 +31,9 @@ def test_insert_object(self): self.assertEqual(len(ds.trajectories), 1) # Check the object was inserted correctly. - t0 = ds.stack.get_single_image(0).get_time() + t0 = ds.stack.get_single_image(0).get_obstime() for i in range(ds.stack.img_count()): - dt = ds.stack.get_single_image(i).get_time() - t0 + dt = ds.stack.get_single_image(i).get_obstime() - t0 px = int(trj.x + dt * trj.x_v + 0.5) py = int(trj.y + dt * trj.y_v + 0.5) diff --git a/tests/test_image_stack.py b/tests/test_image_stack.py index b157ddcf..e2d707b9 100644 --- a/tests/test_image_stack.py +++ b/tests/test_image_stack.py @@ -33,12 +33,12 @@ def test_create(self): self.assertEqual(self.num_images, self.im_stack.img_count()) self.assertEqual(self.im_stack.get_height(), 60) self.assertEqual(self.im_stack.get_width(), 80) - self.assertEqual(self.im_stack.get_ppi(), 60 * 80) + self.assertEqual(self.im_stack.get_npixels(), 60 * 80) def test_access(self): # Test we can access an individual image. img = self.im_stack.get_single_image(1) - self.assertEqual(img.get_time(), 2.0) + self.assertEqual(img.get_obstime(), 2.0) self.assertEqual(img.get_name(), "layered_test_1") # Test an out of bounds access. diff --git a/tests/test_layered_image.py b/tests/test_layered_image.py index 85207e1b..539c8a89 100644 --- a/tests/test_layered_image.py +++ b/tests/test_layered_image.py @@ -25,8 +25,8 @@ def test_create(self): self.assertIsNotNone(self.image) self.assertEqual(self.image.get_width(), 80) self.assertEqual(self.image.get_height(), 60) - self.assertEqual(self.image.get_ppi(), 80 * 60) - self.assertEqual(self.image.get_time(), 10.0) + self.assertEqual(self.image.get_npixels(), 80 * 60) + self.assertEqual(self.image.get_obstime(), 10.0) self.assertEqual(self.image.get_name(), "layered_test") # Create a fake layered_image. @@ -56,11 +56,11 @@ def test_create_from_layers(self): mask.set_all(0.0) # Create the layered image. - img2 = layered_image(sci, var, mask, 0.0, psf(2.0)) + img2 = layered_image(sci, var, mask, psf(2.0)) self.assertEqual(img2.get_width(), 30) self.assertEqual(img2.get_height(), 40) - self.assertEqual(img2.get_ppi(), 30.0 * 40.0) - self.assertEqual(img2.get_time(), 0.0) + self.assertEqual(img2.get_npixels(), 30.0 * 40.0) + self.assertEqual(img2.get_obstime(), -1.0) # No time given # Check the layers. science = img2.get_science() @@ -72,13 +72,6 @@ def test_create_from_layers(self): self.assertEqual(variance.get_pixel(x, y), 1.0) self.assertAlmostEqual(science.get_pixel(x, y), x + 30.0 * y) - def test_set_time(self): - self.assertIsNotNone(self.image) - self.assertEqual(self.image.get_time(), 10.0) - - self.image.set_time(20.0) - self.assertEqual(self.image.get_time(), 20.0) - def test_add_object(self): science = self.image.get_science() science_50_50 = science.get_pixel(50, 50) @@ -323,12 +316,15 @@ def test_read_write_files(self): im2 = layered_image(full_path, self.p) self.assertEqual(im1.get_height(), im2.get_height()) self.assertEqual(im1.get_width(), im2.get_width()) - self.assertEqual(im1.get_ppi(), im2.get_ppi()) + self.assertEqual(im1.get_npixels(), im2.get_npixels()) + self.assertEqual(im1.get_obstime(), im2.get_obstime()) sci1 = im1.get_science() + sci2 = im2.get_science() + self.assertEqual(sci1.get_obstime(), sci2.get_obstime()) + var1 = im1.get_variance() mask1 = im1.get_mask() - sci2 = im2.get_science() var2 = im2.get_variance() mask2 = im2.get_mask() for x in range(im1.get_width()): diff --git a/tests/test_raw_image.py b/tests/test_raw_image.py index cbf0b5e9..3b2ce688 100644 --- a/tests/test_raw_image.py +++ b/tests/test_raw_image.py @@ -14,15 +14,17 @@ def setUp(self): for x in range(self.width): for y in range(self.height): self.img.set_pixel(x, y, float(x + y * self.width)) + self.img.set_obstime(10.0) def test_create(self): self.assertEqual(self.img.get_width(), self.width) self.assertEqual(self.img.get_height(), self.height) - self.assertEqual(self.img.get_ppi(), self.width * self.height) + self.assertEqual(self.img.get_npixels(), self.width * self.height) for x in range(self.width): for y in range(self.height): self.assertTrue(self.img.pixel_has_data(x, y)) self.assertEqual(self.img.get_pixel(x, y), float(x + y * self.width)) + self.assertEqual(self.img.get_obstime(), 10.0) # Pixels outside the image have no data. self.assertFalse(self.img.pixel_has_data(-1, 5)) @@ -64,17 +66,20 @@ def test_copy(self): img2 = raw_image(self.img) self.assertEqual(img2.get_width(), self.width) self.assertEqual(img2.get_height(), self.height) - self.assertEqual(img2.get_ppi(), self.width * self.height) + self.assertEqual(img2.get_npixels(), self.width * self.height) self.assertTrue(self.img.approx_equal(img2, 0.0001)) + self.assertEqual(img2.get_obstime(), 10.0) - # Set the old image to all zeros. + # Set the old image to all zeros and change the time. self.img.set_all(0.0) + self.img.set_obstime(1.0) # Check the new image is still set correctly. for x in range(self.width): for y in range(self.height): self.assertTrue(img2.pixel_has_data(x, y)) self.assertEqual(img2.get_pixel(x, y), float(x + y * self.width)) + self.assertEqual(img2.get_obstime(), 10.0) def test_set_all(self): self.img.set_all(15.0) @@ -395,6 +400,49 @@ def test_make_stamp(self): stamp.get_pixel(2 + x, 2 + y), float((x + 2) + (y + 2) * self.width), delta=0.001 ) + # Check that the stamp has the same obstime. + self.assertEqual(stamp.get_obstime(), 10.0) + + def test_read_write_file(self): + with tempfile.TemporaryDirectory() as dir_name: + file_name = "tmp_raw_image" + full_path = "%s/%s.fits" % (dir_name, file_name) + + self.img.save_fits(full_path) + + # Reload the file. + img2 = raw_image(0, 0) + img2.load_fits(full_path, 0) + self.assertEqual(img2.get_width(), self.width) + self.assertEqual(img2.get_height(), self.height) + self.assertEqual(img2.get_npixels(), self.width * self.height) + self.assertEqual(img2.get_obstime(), 10.0) + self.assertTrue(self.img.approx_equal(img2, 1e-5)) + + def test_stack_file(self): + with tempfile.TemporaryDirectory() as dir_name: + file_name = "tmp_raw_image" + full_path = "%s/%s.fits" % (dir_name, file_name) + + # Save the image and create a file. + self.img.save_fits(full_path) + + # Add 4 more layers at different times. + for i in range(1, 5): + self.img.set_obstime(10.0 + 2.0 * i) + self.img.append_fits_layer(full_path) + + # Check that we get 5 layers with the correct times. + img2 = raw_image(0, 0) + for i in range(5): + img2.load_fits(full_path, i) + + self.assertEqual(img2.get_width(), self.width) + self.assertEqual(img2.get_height(), self.height) + self.assertEqual(img2.get_npixels(), self.width * self.height) + self.assertEqual(img2.get_obstime(), 10.0 + 2.0 * i) + self.assertTrue(self.img.approx_equal(img2, 1e-5)) + def test_create_median_image(self): img1 = raw_image(np.array([[0.0, -1.0], [2.0, 1.0], [0.7, 3.1]])) img2 = raw_image(np.array([[1.0, 0.0], [1.0, 3.5], [4.0, 3.0]])) diff --git a/tests/test_readme_example.py b/tests/test_readme_example.py index 19830ff4..c4da5849 100644 --- a/tests/test_readme_example.py +++ b/tests/test_readme_example.py @@ -16,7 +16,7 @@ def test_make_and_copy(self): imgs = ds.stack.get_images() # Get the timestamp of the first image. - t0 = imgs[0].get_time() + t0 = imgs[0].get_obstime() # Specify an artificial object flux = 275.0 @@ -25,7 +25,7 @@ def test_make_and_copy(self): # 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) # Create a new image stack with the inserted object.