Skip to content

Commit

Permalink
Use distance axis for PV image from polyline
Browse files Browse the repository at this point in the history
  • Loading branch information
@pford
pford committed Jun 7, 2024
1 parent 5c7337f commit 392d41b
Show file tree
Hide file tree
Showing 3 changed files with 63 additions and 29 deletions.
57 changes: 38 additions & 19 deletions src/ImageGenerators/PvGenerator.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,8 @@ void PvGenerator::SetFileName(int index, const std::string& filename, bool is_pr
}

bool PvGenerator::GetPvImage(const std::shared_ptr<Frame>& frame, const casacore::Matrix<float>& pv_data, casacore::IPosition& pv_shape,
const casacore::Quantity& offset_increment, int start_chan, int stokes, bool reverse, GeneratedImage& pv_image, std::string& message) {
PositionAxisType axis_type, const casacore::Quantity& position_increment, int start_chan, int stokes, bool reverse,
GeneratedImage& pv_image, std::string& message) {
// Create PV image with input data.
auto input_csys = frame->CoordinateSystem();
if (!input_csys->hasSpectralAxis()) {
Expand All @@ -32,12 +33,26 @@ bool PvGenerator::GetPvImage(const std::shared_ptr<Frame>& frame, const casacore
return false;
}

// Convert spectral reference pixel value to world coordinates
// Position axis
int linear_axis = (reverse ? 1 : 0);
std::string position_name;
double position_refval(0.0), position_refpix(0.0);
if (axis_type == OFFSET) {
// refval 0 at center pixel
position_name = "Offset";
position_refpix = (pv_shape[linear_axis] - 1) / 2;
} else {
// refval 0 at first pixel
position_name = "Distance";
}

// Spectral axis: Convert spectral reference pixel value to world coordinates
double spectral_refval, spectral_pixval(start_chan);
input_csys->spectralCoordinate().toWorld(spectral_refval, spectral_pixval);

// Create casacore::TempImage coordinate system and other image info
auto image = SetupPvImage(frame->GetImage(), input_csys, pv_shape, stokes, offset_increment, spectral_refval, reverse, message);
auto image = SetupPvImage(frame->GetImage(), input_csys, pv_shape, position_name, position_increment, position_refval, position_refpix,
spectral_refval, stokes, reverse, message);
if (!image) {
message = "PV image setup failed.";
return false;
Expand Down Expand Up @@ -73,10 +88,12 @@ void PvGenerator::SetPvImageName(const std::string& filename, int index, bool is

std::shared_ptr<casacore::ImageInterface<casacore::Float>> PvGenerator::SetupPvImage(
const std::shared_ptr<casacore::ImageInterface<float>>& input_image, const std::shared_ptr<casacore::CoordinateSystem>& input_csys,
casacore::IPosition& pv_shape, int stokes, const casacore::Quantity& offset_increment, double spectral_refval, bool reverse,
std::string& message) {
casacore::IPosition& pv_shape, const std::string& position_name, const casacore::Quantity& position_increment, double position_refval,
double position_refpix, double spectral_refval, int stokes, bool reverse, std::string& message) {
// Create temporary image (no data) using input image. Return casacore::TempImage.
casacore::CoordinateSystem pv_csys = GetPvCoordinateSystem(input_csys, pv_shape, stokes, offset_increment, spectral_refval, reverse);
// Position axis
casacore::CoordinateSystem pv_csys = GetPvCoordinateSystem(
input_csys, pv_shape, position_name, position_increment, position_refval, position_refpix, spectral_refval, stokes, reverse);

std::shared_ptr<casacore::ImageInterface<casacore::Float>> image(
new casacore::TempImage<casacore::Float>(casacore::TiledShape(pv_shape), pv_csys));
Expand All @@ -96,35 +113,37 @@ std::shared_ptr<casacore::ImageInterface<casacore::Float>> PvGenerator::SetupPvI
}

casacore::CoordinateSystem PvGenerator::GetPvCoordinateSystem(const std::shared_ptr<casacore::CoordinateSystem>& input_csys,
casacore::IPosition& pv_shape, int stokes, const casacore::Quantity& offset_increment, double spectral_refval, bool reverse) {
casacore::IPosition& pv_shape, const std::string& position_name, const casacore::Quantity& position_increment, double position_refval,
double position_refpix, double spectral_refval, int stokes, bool reverse) {
// Set PV coordinate system with LinearCoordinate and input coordinates for spectral and stokes
casacore::CoordinateSystem pv_csys;
int offset_axis = reverse ? 1 : 0;

// Add linear coordinate (offset); needs to have 2 axes or pc matrix will fail in wcslib.
// Add linear coordinate (offset or distance); needs to have 2 axes or pc matrix will fail in wcslib.
// Will remove degenerate linear axis below
casacore::Vector<casacore::String> name(2, "Offset");
casacore::Vector<casacore::String> unit(2, offset_increment.getUnit());
casacore::Vector<casacore::Double> crval(2, 0.0); // center offset is 0
casacore::Vector<casacore::Double> inc(2, offset_increment.getValue());
casacore::Vector<casacore::String> name(2, position_name);
casacore::Vector<casacore::String> unit(2, position_increment.getUnit());
casacore::Vector<casacore::Double> refval(2, position_refval);
casacore::Vector<casacore::Double> inc(2, position_increment.getValue());
casacore::Matrix<casacore::Double> pc(2, 2, 1);
pc(0, 1) = 0.0;
pc(1, 0) = 0.0;
casacore::Vector<casacore::Double> crpix(2, (pv_shape[offset_axis] - 1) / 2);
casacore::LinearCoordinate linear_coord(name, unit, crval, inc, pc, crpix);
casacore::Vector<casacore::Double> refpix(2, position_refpix);
casacore::LinearCoordinate linear_coord(name, unit, refval, inc, pc, refpix);

// Set spectral coordinate
auto& input_spectral_coord = input_csys->spectralCoordinate();
auto freq_type = input_spectral_coord.frequencySystem();
auto freq_inc = input_spectral_coord.increment()(0);
auto rest_freq = input_spectral_coord.restFrequency();
casacore::Double refpix(0.0);
casacore::SpectralCoordinate spectral_coord(freq_type, spectral_refval, freq_inc, refpix, rest_freq);
casacore::Double spectral_refpix(0.0);
casacore::SpectralCoordinate spectral_coord(freq_type, spectral_refval, freq_inc, spectral_refpix, rest_freq);

// Add offset and spectral coordinates
// Add offset and spectral coordinates in axis order
int linear_axis(0);
if (reverse) {
pv_csys.addCoordinate(spectral_coord);
pv_csys.addCoordinate(linear_coord);
linear_axis = 1;
} else {
pv_csys.addCoordinate(linear_coord);
pv_csys.addCoordinate(spectral_coord);
Expand All @@ -140,7 +159,7 @@ casacore::CoordinateSystem PvGenerator::GetPvCoordinateSystem(const std::shared_
}

// Remove second linear axis
pv_csys.removeWorldAxis(offset_axis + 1, 0.0);
pv_csys.removeWorldAxis(linear_axis + 1, 0.0);

pv_csys.setObsInfo(input_csys->obsInfo());

Expand Down
17 changes: 11 additions & 6 deletions src/ImageGenerators/PvGenerator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,27 +20,32 @@ namespace carta {

class PvGenerator {
public:
enum PositionAxisType {OFFSET, DISTANCE};

PvGenerator();
~PvGenerator() = default;

// For PV generator (not preview)
void SetFileName(int index, const std::string& filename, bool is_preview = false);

// Create generated PV image from input data. If reverse, [spectral, offset] instead of normal [offset, spectral].
// Create generated PV image from input data.
// cut_region_type determines position axis as offset from center (line) or distance from beginning (polyline).
// reverse determines order of position and spectral axes.
// Returns generated image and success, with message if failure.
bool GetPvImage(const std::shared_ptr<Frame>& frame, const casacore::Matrix<float>& pv_data, casacore::IPosition& pv_shape,
const casacore::Quantity& offset_increment, int start_chan, int stokes, bool reverse, GeneratedImage& pv_image,
std::string& message);
PositionAxisType axis_type, const casacore::Quantity& position_increment, int start_chan, int stokes, bool reverse,
GeneratedImage& pv_image, std::string& message);

private:
void SetPvImageName(const std::string& filename, int index, bool is_preview);

std::shared_ptr<casacore::ImageInterface<casacore::Float>> SetupPvImage(
const std::shared_ptr<casacore::ImageInterface<float>>& input_image, const std::shared_ptr<casacore::CoordinateSystem>& input_csys,
casacore::IPosition& pv_shape, int stokes, const casacore::Quantity& offset_increment, double spectral_refval, bool reverse,
std::string& message);
casacore::IPosition& pv_shape, const std::string& position_name, const casacore::Quantity& position_increment,
double position_refval, double position_refpix, double spectral_refval, int stokes, bool reverse, std::string& message);
casacore::CoordinateSystem GetPvCoordinateSystem(const std::shared_ptr<casacore::CoordinateSystem>& input_csys,
casacore::IPosition& pv_shape, int stokes, const casacore::Quantity& offset_increment, double spectral_refval, bool reverse);
casacore::IPosition& pv_shape, const std::string& position_name, const casacore::Quantity& position_increment,
double position_refval, double position_refpix, double spectral_refval, int stokes, bool reverse);

// GeneratedImage parameters
int _file_id;
Expand Down
18 changes: 14 additions & 4 deletions src/Region/RegionHandler.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1292,13 +1292,15 @@ bool RegionHandler::CalculatePvPreviewImage(int frame_id, int preview_id, bool q
RemoveRegion(preview_cut_id);

// Use PvGenerator to set PV image for headers only
PvGenerator::PositionAxisType pos_axis_type = (preview_cut_state.type == CARTA::LINE ? PvGenerator::OFFSET : PvGenerator::DISTANCE);
casacore::Matrix<float> no_preview_data; // do not copy actual preview data into image
int start_channel(0); // spectral range applied in preview image
int stokes(preview_cube->GetStokes());
PvGenerator pv_generator;
pv_generator.SetFileName(preview_id, preview_cube->GetSourceFileName(), true);

if (pv_generator.GetPvImage(preview_frame, no_preview_data, data_shape, increment, start_channel, stokes, reverse, pv_image, error)) {
if (pv_generator.GetPvImage(
preview_frame, no_preview_data, data_shape, pos_axis_type, increment, start_channel, stokes, reverse, pv_image, error)) {
int width = data_shape(0);
int height = data_shape(1);

Expand Down Expand Up @@ -1341,13 +1343,20 @@ bool RegionHandler::CalculatePvPreviewImage(int frame_id, int preview_id, bool q

bool RegionHandler::CalculatePvImage(int file_id, int region_id, int width, AxisRange& spectral_range, bool reverse, bool keep,
std::shared_ptr<Frame>& frame, GeneratorProgressCallback progress_callback, CARTA::PvResponse& pv_response, GeneratedImage& pv_image) {
// Generate PV image by approximating line as box regions and getting spectral profile for each.
// Generate PV image by approximating line/polyline as box regions and getting spectral profile for each.
// Sends updates via progress callback.
// Return parameters: PvResponse, GeneratedImage, and preview data if is preview.
// Returns whether PV image was generated.
pv_response.set_success(false);
pv_response.set_cancel(false);

auto region = GetRegion(region_id);
if (!region) {
pv_response.set_message("PV cut region not set");
return false;
}
auto cut_region_type = region->GetRegionState().type;

// Reset stop flag
_stop_pv[file_id] = false;

Expand All @@ -1364,6 +1373,7 @@ bool RegionHandler::CalculatePvImage(int file_id, int region_id, int width, Axis
if (GetLineProfiles(file_id, region_id, width, spectral_range, per_z, stokes_index, "", progress_callback, pv_data, offset_increment,
cancelled, message, reverse)) {
auto pv_shape = pv_data.shape();
PvGenerator::PositionAxisType pos_axis_type = (cut_region_type == CARTA::LINE ? PvGenerator::OFFSET : PvGenerator::DISTANCE);
int start_chan(spectral_range.from); // Used for reference value in returned PV image

// Set PV index suffix (_pv1, _pv2, etc) to keep previously opened PV image
Expand All @@ -1379,8 +1389,8 @@ bool RegionHandler::CalculatePvImage(int file_id, int region_id, int width, Axis
// Create GeneratedImage in PvGenerator
PvGenerator pv_generator;
pv_generator.SetFileName(name_index, source_filename);
pv_success =
pv_generator.GetPvImage(frame, pv_data, pv_shape, offset_increment, start_chan, stokes_index, reverse, pv_image, message);
pv_success = pv_generator.GetPvImage(
frame, pv_data, pv_shape, pos_axis_type, offset_increment, start_chan, stokes_index, reverse, pv_image, message);
cancelled &= _stop_pv[file_id];

// Cleanup
Expand Down

0 comments on commit 392d41b

Please sign in to comment.