Skip to content

Commit

Permalink
Increase the precision of beam deconvolved image by five times
Browse files Browse the repository at this point in the history
  • Loading branch information
markccchiang committed May 24, 2024
1 parent df6936e commit d08d239
Show file tree
Hide file tree
Showing 2 changed files with 31 additions and 10 deletions.
35 changes: 27 additions & 8 deletions src/ImageFitter/ImageFitter.cc
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@

using namespace carta;

ImageFitter::ImageFitter() {
ImageFitter::ImageFitter() : _precision(5) {
_fdf.f = FuncF;
_fdf.df = nullptr; // internally computed using finite difference approximations of f when set to NULL
_fdf.fvv = nullptr;
Expand Down Expand Up @@ -159,18 +159,20 @@ bool ImageFitter::GetDeconvolvedResults(casacore::ImageInterface<float>* image,

// Create model data for deconvolved fit results
if (_create_deconvolved_model_data && !pixel_results.empty()) {
width *= _precision;
height *= _precision;
size_t data_size = width * height;
_deconvolved_model_data.resize(data_size, 0);
for (auto gauss : pixel_results) {
double sigma_x = gauss.major.getValue() / 2.355;
double sigma_y = gauss.minor.getValue() / 2.355;
double sigma_x = _precision * gauss.major.getValue() / 2.355;
double sigma_y = _precision * gauss.minor.getValue() / 2.355;
double theta = gauss.pa.getValue();

for (int i = 0; i < data_size; ++i) {
size_t row = i % width;
size_t col = i / width;
double x = offset_x + (double)row - gauss.center_x.getValue();
double y = offset_y + (double)col - gauss.center_y.getValue();
double x = _precision * offset_x + (double)row - _precision * gauss.center_x.getValue();
double y = _precision * offset_y + (double)col - _precision * gauss.center_y.getValue();
double xp = x * std::cos(theta) + y * std::sin(theta);
double yp = -x * std::sin(theta) + y * std::cos(theta);
_deconvolved_model_data[i] += 2 * casacore::C::pi * sigma_x * sigma_y * gauss.amplitude *
Expand Down Expand Up @@ -441,17 +443,34 @@ std::string ImageFitter::GetLog() {
}

casa::SPIIF ImageFitter::GetImageData(
casa::SPIIF image, const casacore::ImageRegion& image_region, std::vector<float> image_data, bool remove_beam_info) {
casa::SPIIF image, const casacore::ImageRegion& image_region, std::vector<float> image_data, bool is_beam_deconvolved) {
casa::SPIIF sub_image(new casacore::SubImage<casacore::Float>(*image, image_region));
casacore::CoordinateSystem csys = sub_image->coordinates();
casacore::IPosition shape = sub_image->shape();
if (csys.hasDirectionCoordinate() && is_beam_deconvolved) {
if (csys.showType(0) == "Direction") {
auto dir_coord = csys.directionCoordinate();
auto pix_deltas = dir_coord.increment();
pix_deltas /= (double)_precision;
dir_coord.setIncrement(pix_deltas);
auto pix_refs = dir_coord.referencePixel();
pix_refs = pix_refs * (double)_precision + (double)_precision / 2;
dir_coord.setReferencePixel(pix_refs);
csys.replaceCoordinate(dir_coord, 0);
}
if (shape.size() > 1) {
shape[0] *= _precision;
shape[1] *= _precision;
}
}

casa::SPIIF output_image(new casacore::TempImage<casacore::Float>(casacore::TiledShape(shape), csys));
output_image->setUnits(sub_image->units());
output_image->setMiscInfo(sub_image->miscInfo());
output_image->appendLog(sub_image->logger());

auto image_info = sub_image->imageInfo();
if (remove_beam_info) {
if (is_beam_deconvolved) {
image_info.removeRestoringBeam();
} else if (image_info.hasMultipleBeams()) {
// Use first beam, per imageanalysis ImageCollapser
Expand All @@ -464,7 +483,7 @@ casa::SPIIF ImageFitter::GetImageData(
casacore::Array<float> data_array(shape, image_data.data());
output_image->put(data_array);

if (sub_image->isMasked()) {
if (sub_image->isMasked() && !is_beam_deconvolved) {
output_image->makeMask("mask0", true, true);
casacore::Array<casacore::Bool> mask_array = sub_image->getMask();
casacore::Lattice<casacore::Bool>& mask_out = output_image->pixelMask();
Expand Down
6 changes: 4 additions & 2 deletions src/ImageFitter/ImageFitter.h
Original file line number Diff line number Diff line change
Expand Up @@ -148,6 +148,8 @@ class ImageFitter {
const size_t _max_iter = 200;
/** @brief Callback function for updating fitting progress. */
GeneratorProgressCallback _progress_callback;
/** @brief A factor used to change the delta value of direction axes from the original image header **/
size_t _precision;

/**
* @brief Calculate the number of NaN values and standard deviation of the image data.
Expand Down Expand Up @@ -184,11 +186,11 @@ class ImageFitter {
* @param image The casacore ImageInterface object of the entire image
* @param image_region The fitting region
* @param image_data The image data for the generated casacore ImageInterface object
* @param remove_beam_info Whether to remove the beam information from image header
* @param is_beam_deconvolved Whether to remove the beam information from image header
* @return A casacore ImageInterface object
*/
casa::SPIIF GetImageData(
casa::SPIIF image, const casacore::ImageRegion& image_region, std::vector<float> image_data, bool remove_beam_info = false);
casa::SPIIF image, const casacore::ImageRegion& image_region, std::vector<float> image_data, bool is_beam_deconvolved = false);
/**
* @brief Generate filenames by adding a suffix.
* @param filename Name of the fitting image file
Expand Down

0 comments on commit d08d239

Please sign in to comment.