Set rest frequency in moment generator

CHANGELOG.md

@@ -7,11 +7,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
-* Added support for PV image generation along a polyline region ([#1341](https://github.com/CARTAvis/carta-backend/issues/1341)).
-* Added support for loading remote FITS files from the hips2fits server ([#1379](https://github.com/CARTAvis/carta-backend/issues/1379)).
+* Add support for PV image generation along a polyline region ([#1341](https://github.com/CARTAvis/carta-backend/issues/1341)).
+* Add support for loading remote FITS files from the hips2fits server ([#1379](https://github.com/CARTAvis/carta-backend/issues/1379)).
+* Add support for setting rest frequency for moment image generation ([#1385](https://github.com/CARTAvis/carta-backend/issues/1385)).
 
 ### Fixed
-* Fixed crash when loading non-image HDU by URL ([#1365](https://github.com/CARTAvis/carta-backend/issues/1365)).
+* Fix crash when loading non-image HDU by URL ([#1365](https://github.com/CARTAvis/carta-backend/issues/1365)).
 * Fix crash when parsing FITS header long value ([#1366](https://github.com/CARTAvis/carta-backend/issues/1366)).
 * Fix incorrect parsing of SPECSYS value for ATCA FITS header ([#1375](https://github.com/CARTAvis/carta-backend/issues/1375)).
 * Fix hdf5 image distortion after animation stops ([#1368](https://github.com/CARTAvis/carta-backend/issues/1368)).

src/ImageGenerators/MomentGenerator.cc

@@ -28,20 +28,20 @@ bool MomentGenerator::CalculateMoments(int file_id, const casacore::ImageRegion&
     _success = false;
     _cancel = false;
 
-    // Set moment axis
+    // Save request settings
     SetMomentAxis(moment_request);
-
-    // Set pixel range
     SetPixelRange(moment_request);
-
-    // Set moment types
     SetMomentTypes(moment_request);
+    SetRestFrequency(moment_request);
 
-    // Reset an ImageMoments
-    ResetImageMoments(image_region);
+    // Save image rest frequency for restore
+    double image_rest_freq = _image->coordinates().spectralCoordinate().restFrequency(); // Hz
 
     // Calculate moments
     try {
+        // Reset an ImageMoments rest frequency and subimage
+        ResetImageMoments(image_region);
+
         // Start the timer
         _start_time = std::chrono::high_resolution_clock::now();
 
@@ -90,9 +90,12 @@ bool MomentGenerator::CalculateMoments(int file_id, const casacore::ImageRegion&
             }
         }
     } catch (AipsError& error) {
-        _error_msg = error.getLastMessage();
+        _error_msg = error.getMesg();
     }
 
+    // Restore original rest frequency
+    SetImageRestFrequency(image_rest_freq);
+
     // Set is the moment calculation successful or not
     moment_response.set_success(IsSuccess());
 
@@ -181,7 +184,14 @@ void MomentGenerator::SetPixelRange(const CARTA::MomentRequest& moment_request)
     }
 }
 
+void MomentGenerator::SetRestFrequency(const CARTA::MomentRequest& moment_request) {
+    _rest_frequency = moment_request.rest_freq(); // Hz
+}
+
 void MomentGenerator::ResetImageMoments(const casacore::ImageRegion& image_region) {
+    // Set the requested rest frequency in the image spectral coordinate
+    SetImageRestFrequency(_rest_frequency);
+
     // Reset the sub-image
     _sub_image.reset(new casacore::SubImage<casacore::Float>(*_image, image_region));
 
@@ -192,6 +202,19 @@ void MomentGenerator::ResetImageMoments(const casacore::ImageRegion& image_regio
     _image_moments.reset(new IM(casacore::SubImage<casacore::Float>(*_sub_image), os, this, true));
 }
 
+void MomentGenerator::SetImageRestFrequency(double rest_frequency) {
+    auto csys = _image->coordinates();
+    if (rest_frequency != csys.spectralCoordinate().restFrequency()) {
+        casacore::String error;
+        casacore::Quantity new_rest_freq(rest_frequency, "Hz");
+        if (csys.setRestFrequency(error, new_rest_freq)) {
+            _image->setCoordinateInfo(csys);
+        } else {
+            throw(casacore::AipsError(error));
+        }
+    }
+}
+
 int MomentGenerator::GetMomentMode(CARTA::Moment moment) {
     if (_moment_map.count(moment)) {
         return _moment_map[moment];

src/ImageGenerators/MomentGenerator.h

@@ -51,7 +51,9 @@ class MomentGenerator : public casa::ImageMomentsProgressMonitor {
     void SetMomentAxis(const CARTA::MomentRequest& moment_request);
     void SetMomentTypes(const CARTA::MomentRequest& moment_request);
     void SetPixelRange(const CARTA::MomentRequest& moment_request);
+    void SetRestFrequency(const CARTA::MomentRequest& moment_request);
     void ResetImageMoments(const casacore::ImageRegion& image_region);
+    void SetImageRestFrequency(double rest_frequency);
     int GetMomentMode(CARTA::Moment moment);
     casacore::String GetMomentSuffix(casacore::Int moment);
     casacore::String GetInputFileName();
@@ -72,6 +74,7 @@ class MomentGenerator : public casa::ImageMomentsProgressMonitor {
     int _axis;                                // Moment axis
     casacore::Vector<float> _include_pix;
     casacore::Vector<float> _exclude_pix;
+    double _rest_frequency; // Hz
     casacore::String _error_msg;
     bool _success;
     bool _cancel;