Skip to content

Commit

Permalink
Remove the creation of a beam deconvolved image
Browse files Browse the repository at this point in the history
  • Loading branch information
markccchiang committed May 27, 2024
1 parent afdcae8 commit afb7067
Show file tree
Hide file tree
Showing 11 changed files with 72 additions and 121 deletions.
2 changes: 1 addition & 1 deletion carta-protobuf
7 changes: 3 additions & 4 deletions src/Frame/Frame.cc
Original file line number Diff line number Diff line change
Expand Up @@ -1805,8 +1805,7 @@ void Frame::StopMomentCalc() {
}

bool Frame::FitImage(const CARTA::FittingRequest& fitting_request, CARTA::FittingResponse& fitting_response, GeneratedImage& model_image,
GeneratedImage& deconvolved_model_image, GeneratedImage& residual_image, GeneratorProgressCallback progress_callback,
StokesRegion* stokes_region) {
GeneratedImage& residual_image, GeneratorProgressCallback progress_callback, StokesRegion* stokes_region) {
if (!_image_fitter) {
_image_fitter = std::make_unique<ImageFitter>();
}
Expand Down Expand Up @@ -1892,8 +1891,8 @@ bool Frame::FitImage(const CARTA::FittingRequest& fitting_request, CARTA::Fittin
GetImageRegion(file_id, AxisRange(CurrentZ()), CurrentStokes(), output_stokes_region);
}
casa::SPIIF in_image(_loader->GetStokesImage(output_stokes_region.stokes_source));
success = _image_fitter->GetGeneratedImages(in_image, output_stokes_region.image_region, GetFileName(), model_image,
deconvolved_model_image, residual_image, fitting_response);
success = _image_fitter->GetGeneratedImages(
in_image, output_stokes_region.image_region, GetFileName(), model_image, residual_image, fitting_response);
}
}

Expand Down
3 changes: 1 addition & 2 deletions src/Frame/Frame.h
Original file line number Diff line number Diff line change
Expand Up @@ -196,8 +196,7 @@ class Frame {

// Image fitting
bool FitImage(const CARTA::FittingRequest& fitting_request, CARTA::FittingResponse& fitting_response, GeneratedImage& model_image,
GeneratedImage& deconvolved_model_image, GeneratedImage& residual_image, GeneratorProgressCallback progress_callback,
StokesRegion* stokes_region = nullptr);
GeneratedImage& residual_image, GeneratorProgressCallback progress_callback, StokesRegion* stokes_region = nullptr);
void StopFitting();

// Save as a new file or export sub-image to CASA/FITS format
Expand Down
71 changes: 36 additions & 35 deletions src/ImageFitter/Deconvolver.cc
Original file line number Diff line number Diff line change
Expand Up @@ -25,44 +25,44 @@ void Deconvolver::GetDeconvolutionResults(CARTA::FittingResponse& fitting_respon
const CARTA::GaussianComponent& in_gauss = in_gauss_vec[i];
const CARTA::GaussianComponent& in_gauss_error = in_gauss_error_vec[i];
DeconvolutionResult world_result;
if (DoDeconvolution(in_gauss, world_result)) {
if (world_result.amplitude > 0) {
DeconvolutionResult pixel_result;
if (GetWorldWidthToPixel(world_result, pixel_result)) {
// Add a deconvolved gaussian values in pixel unit
auto decon_values = fitting_response.add_result_decon_values();
auto decon_values_center = decon_values->mutable_center();
decon_values_center->set_x(in_gauss.center().x());
decon_values_center->set_y(in_gauss.center().y());
auto decon_values_fwhm = decon_values->mutable_fwhm();
decon_values_fwhm->set_x(pixel_result.major.getValue());
decon_values_fwhm->set_y(pixel_result.minor.getValue());
decon_values->set_amp(in_gauss.amp());
double pa = world_result.pa.getValue() < 0 ? -world_result.pa.getValue() + 90 : world_result.pa.getValue();
decon_values->set_pa(pa);

// Add a deconvolved gaussian errors in pixel unit
auto decon_errors = fitting_response.add_result_decon_errors();
auto decon_errors_center = decon_errors->mutable_center();
decon_errors_center->set_x(in_gauss_error.center().x());
decon_errors_center->set_y(in_gauss_error.center().y());
auto decon_errors_fwhm = decon_errors->mutable_fwhm();
decon_errors_fwhm->set_x(pixel_result.major_err.getValue());
decon_errors_fwhm->set_y(pixel_result.minor_err.getValue());
decon_errors->set_amp(in_gauss_error.amp());
decon_errors->set_pa(world_result.pa_err.getValue());

pixel_results.emplace_back(pixel_result);
}
} else {
auto fit_log = fitting_response.mutable_log();
fit_log->append("\nCan not deconvolved from beam: component is a point source\n");
bool is_point_source;
if (DoDeconvolution(in_gauss, world_result, is_point_source)) {
DeconvolutionResult pixel_result;
if (GetWorldWidthToPixel(world_result, pixel_result)) {
// Add a deconvolved gaussian values in pixel unit
auto decon_values = fitting_response.add_result_decon_values();
auto decon_values_center = decon_values->mutable_center();
decon_values_center->set_x(in_gauss.center().x());
decon_values_center->set_y(in_gauss.center().y());
auto decon_values_fwhm = decon_values->mutable_fwhm();
decon_values_fwhm->set_x(pixel_result.major.getValue());
decon_values_fwhm->set_y(pixel_result.minor.getValue());
decon_values->set_amp(in_gauss.amp());
double pa = world_result.pa.getValue() < 0 ? -world_result.pa.getValue() + 90 : world_result.pa.getValue();
decon_values->set_pa(pa);

// Add a deconvolved gaussian errors in pixel unit
auto decon_errors = fitting_response.add_result_decon_errors();
auto decon_errors_center = decon_errors->mutable_center();
decon_errors_center->set_x(in_gauss_error.center().x());
decon_errors_center->set_y(in_gauss_error.center().y());
auto decon_errors_fwhm = decon_errors->mutable_fwhm();
decon_errors_fwhm->set_x(pixel_result.major_err.getValue());
decon_errors_fwhm->set_y(pixel_result.minor_err.getValue());
decon_errors->set_amp(in_gauss_error.amp());
decon_errors->set_pa(world_result.pa_err.getValue());

pixel_results.emplace_back(pixel_result);
}
}
if (is_point_source) {
auto fit_log = fitting_response.mutable_log();
fit_log->append(fmt::format("\nComponent #{} is a point source.\n", i));
}
}
}

bool Deconvolver::DoDeconvolution(const CARTA::GaussianComponent& in_gauss, DeconvolutionResult& result) {
bool Deconvolver::DoDeconvolution(const CARTA::GaussianComponent& in_gauss, DeconvolutionResult& result, bool& is_point_source) {
casacore::Double center_x = in_gauss.center().x();
casacore::Double center_y = in_gauss.center().y();
casacore::Double fwhm_x = in_gauss.fwhm().x();
Expand All @@ -78,12 +78,13 @@ bool Deconvolver::DoDeconvolution(const CARTA::GaussianComponent& in_gauss, Deco
bool success(false);
casacore::GaussianBeam best_sol(ori_major, ori_minor, ori_pa);
casacore::GaussianBeam best_decon_sol;
casacore::Bool is_point_source(true);
try {
is_point_source = Deconvolve(best_decon_sol, best_sol, _beam);
success = true;
} catch (const casacore::AipsError& x) {
is_point_source = true;
is_point_source = false;
spdlog::error("Fail to deconvolve from fitting results: {}", x.getMesg());
return success;
}

// Calculate errors for deconvolved gaussian from original fit results
Expand Down
2 changes: 1 addition & 1 deletion src/ImageFitter/Deconvolver.h
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@ class Deconvolver {
~Deconvolver() = default;

void GetDeconvolutionResults(CARTA::FittingResponse& fitting_response, std::vector<DeconvolutionResult>& pixel_results);
bool DoDeconvolution(const CARTA::GaussianComponent& in_gauss, DeconvolutionResult& result);
bool DoDeconvolution(const CARTA::GaussianComponent& in_gauss, DeconvolutionResult& result, bool& is_point_source);
bool GetWorldWidthToPixel(const DeconvolutionResult& world_coords, DeconvolutionResult& pixel_coords);
bool WorldWidthToPixel(casacore::Vector<casacore::Double>& pixel_params, const casacore::Vector<casacore::Quantity>& world_params);
GaussianShape PixelToWorld(
Expand Down
35 changes: 3 additions & 32 deletions src/ImageFitter/ImageFitter.cc
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,6 @@ bool ImageFitter::FitImage(size_t width, size_t height, float* image, double bea
_beam_size = beam_size;
_unit = unit;
_create_model_data = create_model_image;
_create_deconvolved_model_data = false;
_create_residual_data = create_residual_image;
_progress_callback = progress_callback;

Expand Down Expand Up @@ -117,14 +116,9 @@ bool ImageFitter::FitImage(size_t width, size_t height, float* image, double bea
}

bool ImageFitter::GetGeneratedImages(casa::SPIIF image, const casacore::ImageRegion& image_region, const std::string& filename,
GeneratedImage& model_image, GeneratedImage& deconvolved_model_image, GeneratedImage& residual_image,
CARTA::FittingResponse& fitting_response) {
GeneratedImage& model_image, GeneratedImage& residual_image, CARTA::FittingResponse& fitting_response) {
if (_create_model_data) {
model_image = GeneratedImage(GetFilename(filename, "model"), GetImageData(image, image_region, _model_data));
if (_create_deconvolved_model_data) {
deconvolved_model_image = GeneratedImage(
GetFilename(filename, "beam_deconvolved_model"), GetImageData(image, image_region, _deconvolved_model_data, true));
}
}
if (_create_residual_data) {
residual_image = GeneratedImage(GetFilename(filename, "residual"), GetImageData(image, image_region, _residual_data));
Expand All @@ -138,7 +132,7 @@ void ImageFitter::StopFitting() {

double ImageFitter::GetResidualRms() {
double residual_rms = std::numeric_limits<double>::quiet_NaN();
if (_create_residual_data) {
if (!_residual_data.empty()) {
gsl_rstat_workspace* rstat_p = gsl_rstat_alloc();
for (int i = 0; i < _residual_data.size(); ++i) {
if (!std::isnan(_residual_data[i])) {
Expand All @@ -157,29 +151,6 @@ bool ImageFitter::GetDeconvolvedResults(casacore::ImageInterface<float>* image,
carta::Deconvolver deconvolver(image->coordinates(), image->imageInfo().restoringBeam(channel, stokes), GetResidualRms());
std::vector<DeconvolutionResult> pixel_results;
deconvolver.GetDeconvolutionResults(fitting_response, pixel_results);

// Create model data for deconvolved fit results
if (_create_model_data) {
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 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 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 *
gsl_ran_gaussian_pdf(xp, sigma_x) * gsl_ran_gaussian_pdf(yp, sigma_y);
}
}
_create_deconvolved_model_data = !pixel_results.empty();
}
return true;
}
return false;
Expand Down Expand Up @@ -404,7 +375,7 @@ void ImageFitter::CalculateImageData(const gsl_vector* residual) {
if (_create_model_data) {
_model_data[i] = _fit_data.data[i] - gsl_vector_get(residual, i);
}
if (_create_residual_data) {
if (_create_residual_data || _beam_size > 0) {
_residual_data[i] = isnan(_fit_data.data[i]) ? _fit_data.data[i] : gsl_vector_get(residual, i);
}
}
Expand Down
7 changes: 1 addition & 6 deletions src/ImageFitter/ImageFitter.h
Original file line number Diff line number Diff line change
Expand Up @@ -98,8 +98,7 @@ class ImageFitter {
* @return Whether the images are successfully generated
*/
bool GetGeneratedImages(std::shared_ptr<casacore::ImageInterface<float>> image, const casacore::ImageRegion& image_region,
const std::string& filename, GeneratedImage& model_image, GeneratedImage& deconvolved_model_image, GeneratedImage& residual_image,
CARTA::FittingResponse& fitting_response);
const std::string& filename, GeneratedImage& model_image, GeneratedImage& residual_image, CARTA::FittingResponse& fitting_response);
/** @brief Stop the ongoing fitting process. */
void StopFitting();
/** @brief Get RMS of residual data from the fit. */
Expand Down Expand Up @@ -133,14 +132,10 @@ class ImageFitter {
FitStatus _fit_status;
/** @brief Whether to create a model image. */
bool _create_model_data;
/** @brief Whether to create a deconvolved model image. */
bool _create_deconvolved_model_data;
/** @brief Whether to create a residual image. */
bool _create_residual_data;
/** @brief Model image data. */
std::vector<float> _model_data;
/** @brief Model data for the deconvolved beam image. */
std::vector<float> _deconvolved_model_data;
/** @brief Residual image data. */
std::vector<float> _residual_data;
/** @brief Maximum number of fitting iterations. */
Expand Down
6 changes: 2 additions & 4 deletions src/Region/RegionHandler.cc
Original file line number Diff line number Diff line change
Expand Up @@ -1513,7 +1513,7 @@ void RegionHandler::ClosePvPreview(int preview_id) {
}

bool RegionHandler::FitImage(const CARTA::FittingRequest& fitting_request, CARTA::FittingResponse& fitting_response,
std::shared_ptr<Frame> frame, GeneratedImage& model_image, GeneratedImage& deconvolved_model_image, GeneratedImage& residual_image,
std::shared_ptr<Frame> frame, GeneratedImage& model_image, GeneratedImage& residual_image,
GeneratorProgressCallback progress_callback) {
int file_id(fitting_request.file_id());
int region_id(fitting_request.region_id());
Expand Down Expand Up @@ -1552,9 +1552,7 @@ bool RegionHandler::FitImage(const CARTA::FittingRequest& fitting_request, CARTA
return false;
}

bool success = false;
success = frame->FitImage(
fitting_request, fitting_response, model_image, deconvolved_model_image, residual_image, progress_callback, &stokes_region);
bool success = frame->FitImage(fitting_request, fitting_response, model_image, residual_image, progress_callback, &stokes_region);

if (region_id == TEMP_FOV_REGION_ID) {
RemoveRegion(region_id);
Expand Down
3 changes: 1 addition & 2 deletions src/Region/RegionHandler.h
Original file line number Diff line number Diff line change
Expand Up @@ -93,8 +93,7 @@ class RegionHandler {

// Image fitting
bool FitImage(const CARTA::FittingRequest& fitting_request, CARTA::FittingResponse& fitting_response, std::shared_ptr<Frame> frame,
GeneratedImage& model_image, GeneratedImage& deconvolved_model_image, GeneratedImage& residual_image,
GeneratorProgressCallback progress_callback);
GeneratedImage& model_image, GeneratedImage& residual_image, GeneratorProgressCallback progress_callback);

private:
// Get unique region id: max id (from 0) + 1
Expand Down
13 changes: 3 additions & 10 deletions src/Session/Session.cc
Original file line number Diff line number Diff line change
Expand Up @@ -1427,7 +1427,6 @@ void Session::OnFittingRequest(const CARTA::FittingRequest& fitting_request, uin
bool success(false);
int region_id(fitting_request.region_id());
GeneratedImage model_image;
GeneratedImage deconvolved_model_image;
GeneratedImage residual_image;

// Set fitting progress callback function
Expand All @@ -1440,23 +1439,17 @@ void Session::OnFittingRequest(const CARTA::FittingRequest& fitting_request, uin
if (!_region_handler) {
_region_handler = std::unique_ptr<RegionHandler>(new RegionHandler());
}
success = _region_handler->FitImage(fitting_request, fitting_response, _frames.at(file_id), model_image,
deconvolved_model_image, residual_image, progress_callback);
success = _region_handler->FitImage(
fitting_request, fitting_response, _frames.at(file_id), model_image, residual_image, progress_callback);
} else {
success = _frames.at(file_id)->FitImage(
fitting_request, fitting_response, model_image, deconvolved_model_image, residual_image, progress_callback);
success = _frames.at(file_id)->FitImage(fitting_request, fitting_response, model_image, residual_image, progress_callback);
}

if (success) {
int next_file_id = GetNextFileId();
if (fitting_request.create_model_image()) {
auto* model_image_open_file_ack = fitting_response.mutable_model_image();
OnOpenFile(next_file_id, model_image.name, model_image.image, model_image_open_file_ack);
if (deconvolved_model_image.image) {
auto* deconvolved_model_image_open_file_ack = fitting_response.mutable_deconvolved_model_image();
OnOpenFile(
++next_file_id, deconvolved_model_image.name, deconvolved_model_image.image, deconvolved_model_image_open_file_ack);
}
}
if (fitting_request.create_residual_image()) {
auto* residual_image_open_file_ack = fitting_response.mutable_residual_image();
Expand Down
Loading

0 comments on commit afb7067

Please sign in to comment.