From 79205ace2377d5bee887eb12922c754738f471c6 Mon Sep 17 00:00:00 2001 From: Martin Falk Date: Tue, 10 Feb 2026 14:36:01 +0100 Subject: [PATCH 01/16] Astrophysics: Astrophysics module --- modules/astrophysics/CMakeLists.txt | 40 ++ modules/astrophysics/depends.cmake | 21 + .../astrophysics/astrophysicsmodule.h | 47 ++ .../astrophysics/astrophysicsmoduledefine.h | 22 + .../processors/fitsvolumeraycaster.h | 67 +++ .../fitsnonlineardepthcomponent.h | 75 +++ modules/astrophysics/python/astroviscommon.py | 186 +++++++ .../python/processors/AxisSliceAnnotation.py | 85 +++ .../processors/FitsPolarizationSource.py | 122 +++++ .../python/processors/FitsVolumeSource.py | 496 ++++++++++++++++++ modules/astrophysics/readme.md | 3 + .../astrophysics/src/astrophysicsmodule.cpp | 52 ++ .../src/processors/fitsvolumeraycaster.cpp | 72 +++ .../fitsnonlineardepthcomponent.cpp | 303 +++++++++++ .../unittests/astrophysics-unittest-main.cpp | 57 ++ 15 files changed, 1648 insertions(+) create mode 100644 modules/astrophysics/CMakeLists.txt create mode 100644 modules/astrophysics/depends.cmake create mode 100644 modules/astrophysics/include/infravis/astrophysics/astrophysicsmodule.h create mode 100644 modules/astrophysics/include/infravis/astrophysics/astrophysicsmoduledefine.h create mode 100644 modules/astrophysics/include/infravis/astrophysics/processors/fitsvolumeraycaster.h create mode 100644 modules/astrophysics/include/infravis/astrophysics/shadercomponents/fitsnonlineardepthcomponent.h create mode 100644 modules/astrophysics/python/astroviscommon.py create mode 100644 modules/astrophysics/python/processors/AxisSliceAnnotation.py create mode 100644 modules/astrophysics/python/processors/FitsPolarizationSource.py create mode 100644 modules/astrophysics/python/processors/FitsVolumeSource.py create mode 100644 modules/astrophysics/readme.md create mode 100644 modules/astrophysics/src/astrophysicsmodule.cpp create mode 100644 modules/astrophysics/src/processors/fitsvolumeraycaster.cpp create mode 100644 modules/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp create mode 100644 modules/astrophysics/tests/unittests/astrophysics-unittest-main.cpp diff --git a/modules/astrophysics/CMakeLists.txt b/modules/astrophysics/CMakeLists.txt new file mode 100644 index 000000000..9cc54cf32 --- /dev/null +++ b/modules/astrophysics/CMakeLists.txt @@ -0,0 +1,40 @@ +ivw_module(AstroPhysics) + +set(HEADER_FILES + include/infravis/astrophysics/astrophysicsmodule.h + include/infravis/astrophysics/astrophysicsmoduledefine.h + include/infravis/astrophysics/processors/fitsvolumeraycaster.h + include/infravis/astrophysics/shadercomponents/fitsnonlineardepthcomponent.h +) +ivw_group("Header Files" ${HEADER_FILES}) + +set(SOURCE_FILES + src/astrophysicsmodule.cpp + src/processors/fitsvolumeraycaster.cpp + src/shadercomponents/fitsnonlineardepthcomponent.cpp +) +ivw_group("Source Files" ${SOURCE_FILES}) + +set(PYTHON_FILES + python/processors/FitsPolarizationSource.py + python/processors/FitsVolumeSource.py + python/astroviscommon.py +) +ivw_group("Python Files" ${PYTHON_FILES}) + +set(SHADER_FILES + # Add shaders +) +ivw_group("Shader Files" ${SHADER_FILES}) + +set(TEST_FILES + tests/unittests/astrophysics-unittest-main.cpp +) +ivw_add_unittest(${TEST_FILES}) + +ivw_create_module(${SOURCE_FILES} ${HEADER_FILES} ${SHADER_FILES} ${PYTHON_FILES}) + +set_property(GLOBAL APPEND PROPERTY IVW_PYTHON_ADDITIONAL_MODULES astropy) + +# Add shader directory to install package +#ivw_add_to_module_pack(${CMAKE_CURRENT_SOURCE_DIR}/glsl) diff --git a/modules/astrophysics/depends.cmake b/modules/astrophysics/depends.cmake new file mode 100644 index 000000000..5c984fd76 --- /dev/null +++ b/modules/astrophysics/depends.cmake @@ -0,0 +1,21 @@ +# Inviwo module dependencies for current module +# List modules on the format "InviwoModule" +set(dependencies + InviwoOpenGLModule + InviwoBaseGLModule + InviwoPython3Module +) + +# Add an alias for this module. Several modules can share an alias. +# Useful if several modules implement the same functionality. +# A depending module can depend on the alias, and a drop down of +# available implementations will be presented in CMake +#set(aliases ModuleAlias) + +# Mark the module as protected to prevent it from being reloaded +# when using runtime module reloading. +#set(protected ON) + +# By calling set(EnableByDefault ON) the module will be set to enabled +# when initially being added to CMake. Default OFF. +set(EnableByDefault ON) \ No newline at end of file diff --git a/modules/astrophysics/include/infravis/astrophysics/astrophysicsmodule.h b/modules/astrophysics/include/infravis/astrophysics/astrophysicsmodule.h new file mode 100644 index 000000000..7d9813fed --- /dev/null +++ b/modules/astrophysics/include/infravis/astrophysics/astrophysicsmodule.h @@ -0,0 +1,47 @@ +/********************************************************************************* + * + * Inviwo - Interactive Visualization Workshop + * + * Copyright (c) 2026 Inviwo Foundation + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + *********************************************************************************/ +#pragma once + +#include +#include +#include +#include + +namespace inviwo { + +class IVW_MODULE_ASTROPHYSICS_API AstroPhysicsModule : public InviwoModule { +public: + explicit AstroPhysicsModule(InviwoApplication* app); + virtual ~AstroPhysicsModule() = default; + + pyutil::ModulePath scripts_; + PythonProcessorFolderObserver pythonFolderObserver_; +}; + +} // namespace inviwo diff --git a/modules/astrophysics/include/infravis/astrophysics/astrophysicsmoduledefine.h b/modules/astrophysics/include/infravis/astrophysics/astrophysicsmoduledefine.h new file mode 100644 index 000000000..1c54b3532 --- /dev/null +++ b/modules/astrophysics/include/infravis/astrophysics/astrophysicsmoduledefine.h @@ -0,0 +1,22 @@ +#pragma once + +// clang-format off +#ifdef INVIWO_ALL_DYN_LINK //DYNAMIC + // If we are building DLL files we must declare dllexport/dllimport + #ifdef IVW_MODULE_ASTROPHYSICS_EXPORTS + #ifdef _WIN32 + #define IVW_MODULE_ASTROPHYSICS_API __declspec(dllexport) + #else //UNIX (GCC) + #define IVW_MODULE_ASTROPHYSICS_API __attribute__ ((visibility ("default"))) + #endif + #else + #ifdef _WIN32 + #define IVW_MODULE_ASTROPHYSICS_API __declspec(dllimport) + #else + #define IVW_MODULE_ASTROPHYSICS_API + #endif + #endif +#else //STATIC + #define IVW_MODULE_ASTROPHYSICS_API +#endif +// clang-format on diff --git a/modules/astrophysics/include/infravis/astrophysics/processors/fitsvolumeraycaster.h b/modules/astrophysics/include/infravis/astrophysics/processors/fitsvolumeraycaster.h new file mode 100644 index 000000000..9349aaef0 --- /dev/null +++ b/modules/astrophysics/include/infravis/astrophysics/processors/fitsvolumeraycaster.h @@ -0,0 +1,67 @@ +/********************************************************************************* + * + * Inviwo - Interactive Visualization Workshop + * + * Copyright (c) 2026 Inviwo Foundation + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + *********************************************************************************/ + +#pragma once + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace inviwo { + +class IVW_MODULE_ASTROPHYSICS_API FitsVolumeRaycaster : public VolumeRaycasterBase { +public: + explicit FitsVolumeRaycaster(std::string_view identifier = {}, std::string_view displayName = {}); + + virtual void process() override; + + virtual const ProcessorInfo& getProcessorInfo() const override; + static const ProcessorInfo processorInfo_; + +private: + FitsNonlinearDepthComponent fits_; + //VolumeComponent volume_; + EntryExitComponent entryExit_; + BackgroundComponent background_; + IsoTFComponent<1> isoTF_; + RaycastingComponent raycasting_; + CameraComponent camera_; + LightComponent light_; +}; + +} // namespace inviwo diff --git a/modules/astrophysics/include/infravis/astrophysics/shadercomponents/fitsnonlineardepthcomponent.h b/modules/astrophysics/include/infravis/astrophysics/shadercomponents/fitsnonlineardepthcomponent.h new file mode 100644 index 000000000..7628dceab --- /dev/null +++ b/modules/astrophysics/include/infravis/astrophysics/shadercomponents/fitsnonlineardepthcomponent.h @@ -0,0 +1,75 @@ +/********************************************************************************* + * + * Inviwo - Interactive Visualization Workshop + * + * Copyright (c) 2026 Inviwo Foundation + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + *********************************************************************************/ +#pragma once + +#include + +#include +#include +#include +#include +#include + +namespace inviwo { + +class Inport; +class Property; +class Shader; +class TextureUnitContainer; + +class IVW_MODULE_ASTROPHYSICS_API FitsNonlinearDepthComponent : public ShaderComponent { +public: + explicit FitsNonlinearDepthComponent(std::string_view volume, Document help = {}); + + virtual std::string_view getName() const override; + virtual void process(Shader& shader, TextureUnitContainer&) override; + virtual void initializeResources(Shader& shader) override; + virtual std::vector> getInports() override; + virtual std::vector getProperties() override; + + virtual std::vector getSegments() override; + + std::optional channelsForVolume() const; + + VolumeInport volumePort_; + +private: + CompositeProperty fitsParameters_; + FloatProperty assumedDistance_; + FloatProperty velocityInf_; + FloatProperty velocityStar_; + FloatProperty crpix3_; + FloatProperty crval3_; + FloatProperty cdelt3_; + FloatProperty restFrequency_; + FloatVec3Property dataExtent_; + FloatVec3Property dataOffset_; +}; + +} // namespace inviwo diff --git a/modules/astrophysics/python/astroviscommon.py b/modules/astrophysics/python/astroviscommon.py new file mode 100644 index 000000000..6bd91f0d1 --- /dev/null +++ b/modules/astrophysics/python/astroviscommon.py @@ -0,0 +1,186 @@ + +import numpy as np +from pathlib import Path +from collections.abc import Callable +from enum import Enum +from dataclasses import dataclass +from functools import lru_cache +from astropy.io import fits + + +class LatLong(Enum): + AbsoluteDeg = 0 + RelativeDeg = 1 + RelativeArcsec = 2 + + +@dataclass +class FitsParams: + naxis: tuple[int, ...] + ctype: tuple[str, str, str] + cunit: tuple[str, str, str] + crpix: tuple[float, float, float] + crval: tuple[float, float, float] + cdelt: tuple[float, float, float] + btype: str + bunit: str + rest_frequency: float # [Hz] + + +def extractFITSParams(header: fits.header.Header) -> FitsParams: + rest_frequency: float = header['restfrq'] # equals 345_339_760_000 + # mean difference 0.0080740896 + # rest_frequency = 345_339_769_300 # [hz] + + return FitsParams( + naxis=tuple(header[f'naxis{i}'] for i in range(1, header['naxis'] + 1)), + ctype=tuple(header[f'ctype{i}'].lower() for i in range(1, 4)), + cunit=tuple(header[f'cunit{i}'] for i in range(1, 4)), + crpix=tuple(header[f'crpix{i}'] for i in range(1, 4)), + crval=tuple(header[f'crval{i}'] for i in range(1, 4)), + cdelt=tuple(header[f'cdelt{i}'] for i in range(1, 4)), + btype=header['btype'], + bunit=header['bunit'], + rest_frequency=rest_frequency) + + +@lru_cache(maxsize=2) +def readFITSData(filepath: Path, + dataset_index: int | None = None, + **kwargs) -> tuple[FitsParams, np.ndarray]: + if not dataset_index: + dataset_index = 0 + with fits.open(filepath, **kwargs) as hdul: + fits_parameters: FitsParams = extractFITSParams(hdul[dataset_index].header) + data: np.ndarray = hdul[dataset_index].data # TODO: copy data? + return (fits_parameters, data) + + +def pixelCoordinatesFunc(params: FitsParams) -> Callable[[np.ndarray], [int, int]]: + def func(i: int, j: int) -> np.ndarray: + # FIXME: use i or (i+1) ? + alpha = params.crval[0] + params.cdelt[0] * (i - params.crpix[0]) + beta = params.crval[1] + params.cdelt[1] * (j - params.crpix[1]) + return np.array([alpha, beta]) + + if not params.ctype[0] == 'ra---sin': + raise ValueError( + f'Expected unit type "ctype1" to be "RA---SIN". Found "{params.ctype[0]}".') + if not params.ctype[1] == 'dec--sin': + raise ValueError( + f'Expected unit type "ctype2" to be "DEC--SIN". Found "{params.ctype[1]}".') + return func + + +def estimateDepth(params: FitsParams, + velocity: float, + assumed_distance: float = 123.0, + velocity_inf: float = 14.5, + velocity_star: float = -26.5) -> np.ndarray: + """ + Conversion of velocity axis to 3rd coordinate: + - In this case via the assumption of a radial velocity law of the form: + |V - V*| = Vinf * (1-((Rw-r)/Rw)^2.5) for rRw + Vinf = 14.5 km/s + V* = -26.5 km/s + Rw = 54 au + (assumed distance 123 pc, so 1" = 123 au) + + Since in this particular case, we're at scales larger than Rw we're in the + simplest case where the velocity is constant V=Vinf=14.5 km/s. + + So the line of sight velocity Vlos = V - V* = Vinf * (z/r) where r=sqrt(x^2+y^2+z^2) + + taking rho = sqrt(x^2+y^2) and mu = Vlos/Vinf we have + Z = +- rho * mu/sqrt(1-mu^2) + + The positive and negative sign are for the blue-shifted (in front) and red-shifted + (behind the plane of the sky) parts of the envelope. + We should in this case disregard all velocity channels with |V-V*|=|Vlos|>Vinf. + And remember for the coordinate transformation from (relative to the center of + the emission) arcseconds to au is using 1"=123 au). + + Parameters + ---------- + params : FitsParams + Parameters extracted from a fits.header.Header + velocity : float + in km/s. + assumed_distance : float, optional + in parsecs. The default is 123.0. + velocity_inf: float, optional + in km/s. The default is 14.5. + velocity_star: float, optional + in km/s. The default is -26.5. + + + Returns + ------- + 2D array of depth values for the given velocity + + """ + v_los: float = velocity - velocity_star + if np.abs(v_los) > velocity_inf: + raise ValueError(f'v_los={v_los} exceeds v_inf={velocity_inf}') + + mu: float = v_los / velocity_inf + + coords: Callable[[np.ndarray], [int, int]] = pixelCoordinatesFunc(params) + # coordinates are given in degree, convert to arcseconds then parsecs (au) + delta: np.ndarray = (coords(1, 1) - coords(0, 0)) * \ + 3600.0 * assumed_distance + dims = np.array(params.naxis[:2]) + center = (dims - 1) / 2.0 + + depths = np.empty(dims, dtype=float) + with np.nditer(depths, flags=['multi_index'], op_flags=['writeonly']) as it: + for value in it: + d = (it.multi_index - center) * delta + rho: float = np.sqrt(d[0]**2 + d[1]**2) + value[...] = rho * mu / np.sqrt(1.0 - mu**2) + return depths + + +def getLatLongBasis(params: FitsParams, + lat_long: LatLong = LatLong.RelativeArcsec + ) -> tuple[np.ndarray, np.ndarray]: + """ + Calculate the lat/long extent and offset given a fits header for pixel centers + + Parameters + ---------- + params : FitsParams + Parameters extracted from a fits.header.Header + lat_long : LatLong, optional + Determines the coordinate system (absolute/relative, degree/arcsec, ...). + The default is LatLong.RelativeArcsec. + + Returns + ------- + tuple[np.ndarray, np.ndarray] + lat/long extent and offset + + """ + coords: Callable[[np.ndarray], [int, int]] = pixelCoordinatesFunc(params) + origin: np.ndarray = coords(0, 0) + delta: np.ndarray = coords(1, 1) - origin + extent: np.ndarray = coords(params.naxis[0], params.naxis[1]) - origin + + for i in range(2): + if extent[i] < 0.0: + origin[i] += extent[i] + extent[i] = -extent[i] + delta[i] = -delta[i] + + if lat_long == LatLong.RelativeArcsec: + extent *= 3600 # convert from degree to arcseconds + + if lat_long == LatLong.AbsoluteDeg: + offset = origin - delta * 0.5 + elif lat_long in [LatLong.RelativeDeg, LatLong.RelativeArcsec]: + offset = extent * 0.5 + else: + raise ValueError(f'Invalid Lat/Long value {lat_long}') + + return (extent, offset) diff --git a/modules/astrophysics/python/processors/AxisSliceAnnotation.py b/modules/astrophysics/python/processors/AxisSliceAnnotation.py new file mode 100644 index 000000000..02be38c3c --- /dev/null +++ b/modules/astrophysics/python/processors/AxisSliceAnnotation.py @@ -0,0 +1,85 @@ +# Name: AxisSliceAnnotation + +import inviwopy as ivw + +customLabelHelp: str = r'''Custom formatted label. The following placeholders are defined: ++ `{i}` slice index ++ `{n}` axis name ++ `{v}` axis value of the current slice index ++ `{u}` axis unit''' + + +class AxisSliceAnnotation(ivw.Processor): + """ + Composes a formatted string for a given slice index and axis using the axis information + of the input volume. + """ + + def __init__(self, id, name): + ivw.Processor.__init__(self, id, name) + self.inport = ivw.data.VolumeInport("volume") + self.addInport(self.inport) + + cb = ivw.properties.ConstraintBehavior + self.sliceAxis = \ + ivw.properties.OptionPropertyInt( + "sliceAxis", "Slice Axis", + ivw.md2doc("Defines the volume axis of the slice"), + [ivw.properties.IntOption("x", "X axis", 0), + ivw.properties.IntOption("y", "Y axis", 1), + ivw.properties.IntOption("z", "Z axis", 2)]) + self.slice = ivw.properties.IntProperty("sliceNumber", "Slice", + ivw.md2doc("Index of the slice"), + 128, min=(1, cb.Immutable), max=(256, cb.Ignore)) + self.labelType = \ + ivw.properties.OptionPropertyString( + "labelType", "Type of Slice Label", + [ivw.properties.StringOption("slice", "Slice Index", "{i}"), + ivw.properties.StringOption("axisValueUnit", "Axis Value + Unit", "{v} [{u}]"), + ivw.properties.StringOption("axisValue", "Axis Value only", "{v}"), + ivw.properties.StringOption("custom", + "Custom Format (example 'Slice {i}, {n} {v} [{u}]')", + "custom")]) + self.customLabel = ivw.properties.StringProperty( + "customLabel", "Custom Label", ivw.md2doc(customLabelHelp), "Slice {i}, {n} {v} [{u}]") + self.output = ivw.properties.StringProperty( + "output", "Output", + ivw.md2doc("Resulting label for the given slice and axis."), "") + self.output.readOnly = True + + self.addProperties([self.sliceAxis, self.slice, self.labelType, + self.customLabel, self.output]) + + @staticmethod + def processorInfo(): + return ivw.ProcessorInfo( + classIdentifier="org.inviwo.AxisSliceAnnotation", + displayName="Axis Slice Annotation", + category="Python", + codeState=ivw.CodeState.Stable, + tags=ivw.Tags([ivw.Tags.PY, ivw.Tag("Volume")]), + help=ivw.unindentMd2doc(AxisSliceAnnotation.__doc__) + ) + + def getProcessorInfo(self): + return AxisSliceAnnotation.processorInfo() + + def initializeResources(self): + pass + + def process(self): + volume = self.inport.getData() + + axis = self.sliceAxis.value + value: float = (volume.basis[axis] * (self.slice.value - 1) + / volume.dimensions[axis] + volume.offset)[axis] + + formatString = self.labelType.selectedValue \ + if self.labelType.selectedValue != "custom" else self.customLabel.value + + self.output.value = formatString.format( + i=self.slice.value, + v=value, + n=volume.axes[axis].name, + u=volume.axes[axis].unit + ) diff --git a/modules/astrophysics/python/processors/FitsPolarizationSource.py b/modules/astrophysics/python/processors/FitsPolarizationSource.py new file mode 100644 index 000000000..d2652035c --- /dev/null +++ b/modules/astrophysics/python/processors/FitsPolarizationSource.py @@ -0,0 +1,122 @@ +# Name: FitsPolarizationSource + +import inviwopy as ivw +import ivwbase + +import os +import numpy as np +from astropy.io import fits + + +class FitsPolarizationSource(ivw.Processor): + """ + Reads two FITS files containing the angle and intensity of the polarization and + outputs the resulting 2D vector field + """ + + def __init__(self, id, name): + ivw.Processor.__init__(self, id, name) + self.deserializing = False + + self.outport = ivw.data.VolumeOutport("polarization") + self.addOutport(self.outport) + + self.filenameAngle = ivw.properties.FileProperty( + "filenameAngle", "Angle", "", "fits") + self.filenameAngle.addNameFilter(ivw.properties.FileExtension("fits", "FITS file format")) + self.filenameIntensity = ivw.properties.FileProperty( + "filenameIntensity", "Intensity", "", "fits") + self.filenameIntensity.addNameFilter( + ivw.properties.FileExtension("fits", "FITS file format")) + self.basis = ivwbase.properties.BasisProperty("basis", "Basis and offset") + self.information = ivwbase.properties.VolumeInformationProperty( + "information", "Data information") + self.addProperties([self.filenameAngle, self.filenameIntensity, + self.basis, self.information]) + + @staticmethod + def processorInfo(): + return ivw.ProcessorInfo( + classIdentifier="org.inviwo.FitsPolarizationSource", + displayName="FitsPolarizationSource", + category="Python", + codeState=ivw.CodeState.Stable, + tags=ivw.Tags.PY, + help=ivw.unindentMd2doc(FitsPolarizationSource.__doc__) + ) + + def getProcessorInfo(self): + return FitsPolarizationSource.processorInfo() + + def initializeResources(self): + pass + + @staticmethod + def createVectorField(angle: np.array, length: np.array, + use_degree: bool = True) -> ivw.data.Volume: + + # angle = np.ma.fix_invalid(angle, fill_value=0.).data + # length = np.ma.fix_invalid(length, fill_value=1.0).data + + if use_degree: + angle = angle * np.pi / 180.0 + + v = np.stack([np.cos(angle + np.pi / 2.0) * length, + np.sin(angle + np.pi / 2.0) * length], axis=3) + + volumerep = ivw.data.VolumePy(np.copy(v.astype(dtype=np.float32), order='C')) + + volume = ivw.data.Volume(volumerep) + + # TODO: extract extent from hdul[0].header + volume.basis = ivw.glm.mat3(1) + volume.offset = ivw.glm.vec3(-0.5, -0.5, -0.5) + + max_len: float = np.nanmax(length) + datarange = ivw.glm.dvec2(-max_len, max_len) + volume.dataMap.dataRange = datarange + volume.dataMap.valueRange = datarange + return volume + + def process(self): + if not self.filenameAngle.value or not self.filenameIntensity: + return + + filenameAngle = ivwbase.io.downloadAndCacheIfUrl(self.filenameAngle.valueAsString()) + if not os.path.isfile(filenameAngle): + return + filenameIntensity = ivwbase.io.downloadAndCacheIfUrl(self.filenameIntensity.valueAsString()) + if not os.path.isfile(filenameIntensity): + return + + with fits.open(filenameAngle) as hdul_angle, fits.open(filenameIntensity) as hdul_intensity: + if hdul_angle[0].data.shape != hdul_intensity[0].data.shape: + raise ValueError(f'Shape differs between polarization angles ' + f'({hdul_angle[0].data.shape}) and intensity ' + f'({hdul_intensity[0].data.shape})') + use_degree: bool = hdul_angle[0].header['bunit'].lower() == "deg" + + volume: ivw.data.Volume = FitsPolarizationSource.createVectorField( + hdul_angle[0].data, + hdul_intensity[0].data, + use_degree=use_degree) + + volume.dataMap.valueAxis.name = "Polarity" + # hdul_intensity[0].header['bunit'] is 'Jy/beam ', which is not a valid unit + volume.dataMap.valueAxis.unit = ivw.data.Unit("Jy") + volume.swizzlemask = ivw.data.SwizzleMask.defaultData(2) + + from inviwopy.properties import OverwriteState as ows + + self.basis.updateForNewEntity(volume, self.deserializing) + self.information.updateForNewVolume(volume, ows.Yes if self.deserializing else ows.No) + self.deserializing = False + + self.basis.updateEntity(volume) + self.information.updateVolume(volume) + + self.outport.setData(volume) + + def deserialize(self, deserializer: ivw.Deserializer): + ivw.Processor.deserialize(self, deserializer) + self.deserializing = True diff --git a/modules/astrophysics/python/processors/FitsVolumeSource.py b/modules/astrophysics/python/processors/FitsVolumeSource.py new file mode 100644 index 000000000..b06558262 --- /dev/null +++ b/modules/astrophysics/python/processors/FitsVolumeSource.py @@ -0,0 +1,496 @@ +# Name: FitsVolumeSource + +import inviwopy as ivw +from inviwopy import glm +import ivwbase + +import numpy as np +from pathlib import Path +from collections.abc import Callable +from enum import Enum +from dataclasses import dataclass +from contextlib import suppress +from astropy.io import fits + + +class LatLong(Enum): + AbsoluteDeg = 0 + RelativeDeg = 1 + RelativeArcsec = 2 + + +@dataclass +class FitsParams: + naxis: tuple[int, ...] + ctype: tuple[str, str, str] + cunit: tuple[str, str, str] + crpix: tuple[float, float, float] + crval: tuple[float, float, float] + cdelt: tuple[float, float, float] + btype: str + bunit: str + rest_frequency: float # [Hz] + + +def extractFitsParams(header: fits.header.Header) -> FitsParams: + rest_frequency: float = header['restfrq'] # equals 345_339_760_000 + # mean difference 0.0080740896 + # rest_frequency = 345_339_769_300 # [hz] + + return FitsParams( + naxis=tuple(header[f'naxis{i}'] for i in range(1, header['naxis'] + 1)), + ctype=tuple(header[f'ctype{i}'].lower() for i in range(1, 4)), + cunit=tuple(header[f'cunit{i}'] for i in range(1, 4)), + crpix=tuple(header[f'crpix{i}'] for i in range(1, 4)), + crval=tuple(header[f'crval{i}'] for i in range(1, 4)), + cdelt=tuple(header[f'cdelt{i}'] for i in range(1, 4)), + btype=header['btype'], + bunit=header['bunit'], + rest_frequency=rest_frequency) + + +def pixelCoordinatesFunc(params: FitsParams) -> Callable[[np.ndarray], [int, int]]: + def func(i: int, j: int) -> np.ndarray: + # FIXME: use i or (i+1) ? + alpha = params.crval[0] + params.cdelt[0] * (i - params.crpix[0]) + beta = params.crval[1] + params.cdelt[1] * (j - params.crpix[1]) + return np.array([alpha, beta]) + + if not params.ctype[0] == 'ra---sin': + raise ValueError( + f'Expected unit type "ctype1" to be "RA---SIN". Found "{params.ctype[0]}".') + if not params.ctype[1] == 'dec--sin': + raise ValueError( + f'Expected unit type "ctype2" to be "DEC--SIN". Found "{params.ctype[1]}".') + return func + + +""" +Conversion of velocity axis to 3rd coordinate: +- In this case via the assumption of a radial velocity law of the form: + |V - V*| = Vinf * (1-((Rw-r)/Rw)^2.5) for rRw +Vinf = 14.5 km/s +V* = -26.5 km/s +Rw = 54 au +(assumed distance 123 pc, so 1" = 123 au) + +Since in this particular case, we're at scales larger than Rw we're in the +simplest case where the velocity is constant V=Vinf=14.5 km/s. + +So the line of sight velocity Vlos = V - V* = Vinf * (z/r) where r=sqrt(x^2+y^2+z^2) + +taking rho = sqrt(x^2+y^2) and mu = Vlos/Vinf we have +Z = +- rho * mu/sqrt(1-mu^2) + +The positive and negative sign are for the blue-shifted (in front) and red-shifted +(behind the plane of the sky) parts of the envelope. +We should in this case disregard all velocity channels with |V-V*|=|Vlos|>Vinf. +And remember for the coordinate transformation from (relative to the center of +the emission) arcseconds to au is using 1"=123 au). +""" + + +def estimateDepth(params: FitsParams, + velocity: float, + assumed_distance: float = 123.0, + velocity_inf: float = 14.5, + velocity_star: float = -26.5) -> np.ndarray: + """ + + Parameters + ---------- + params : FitsParams + Parameters extracted from a fits.header.Header + velocity : float + in km/s. + assumed_distance : float, optional + in parsecs. The default is 123.0. + velocity_inf: float, optional + in km/s. The default is 14.5. + velocity_star: float, optional + in km/s. The default is -26.5. + + + Returns + ------- + 2D array of depth values for the given velocity + + """ + v_los: float = velocity - velocity_star + if np.abs(v_los) > velocity_inf: + raise ValueError(f'v_los={v_los} exceeds v_inf={velocity_inf}') + + mu: float = v_los / velocity_inf + + coords: Callable[[np.ndarray], [int, int]] = pixelCoordinatesFunc(params) + # coordinates are given in degree, convert to arcseconds then parsecs (au) + delta: np.ndarray = (coords(1, 1) - coords(0, 0)) * \ + 3600.0 * assumed_distance + dims = np.array(params.naxis[:2]) + center = (dims - 1) / 2.0 + + depths = np.empty(dims, dtype=float) + with np.nditer(depths, flags=['multi_index'], op_flags=['writeonly']) as it: + for value in it: + d = (it.multi_index - center) * delta + rho: float = np.sqrt(d[0]**2 + d[1]**2) + value[...] = rho * mu / np.sqrt(1.0 - mu**2) + return depths + + +def getLatLongBasis(params: FitsParams, + lat_long: LatLong = LatLong.RelativeArcsec + ) -> tuple[np.ndarray, np.ndarray]: + """ + Calculate the lat/long extent and offset given a fits header for pixel centers + + Parameters + ---------- + params : FitsParams + Parameters extracted from a fits.header.Header + lat_long : LatLong, optional + Determines the coordinate system (absolute/relative, degree/arcsec, ...). + The default is LatLong.RelativeArcsec. + + Returns + ------- + tuple[np.ndarray, np.ndarray] + lat/long extent and offset + + """ + coords: Callable[[np.ndarray], [int, int]] = pixelCoordinatesFunc(params) + origin: np.ndarray = coords(0, 0) + delta: np.ndarray = coords(1, 1) - origin + extent: np.ndarray = coords(params.naxis[0], params.naxis[1]) - origin + + for i in range(2): + if extent[i] < 0.0: + origin[i] += extent[i] + extent[i] = -extent[i] + delta[i] = -delta[i] + + if lat_long == LatLong.RelativeArcsec: + extent *= 3600 # convert from degree to arcseconds + + if lat_long == LatLong.AbsoluteDeg: + offset = origin - delta * 0.5 + elif lat_long in [LatLong.RelativeDeg, LatLong.RelativeArcsec]: + offset = extent * 0.5 + else: + raise ValueError(f'Invalid Lat/Long value {lat_long}') + + return (extent, offset) + + +def createDepthMesh(params: FitsParams, + velocity_func: Callable[[float], int], + tf: ivw.data.TransferFunction, + lat_long: LatLong, + assumed_distance: float = 123.0, + velocity_inf: float = 14.5, + velocity_star: float = -26.5 + ) -> (ivw.data.Mesh, np.ndarray): + dims = np.array(params.naxis) + + positions: list[np.ndarray] = [] + + mesh_count: int = 0 + mesh = ivw.data.Mesh(dt=ivw.data.DrawType.Triangles, ct=ivw.data.ConnectivityType.Unconnected) + for i in range(0, dims[2]): + with suppress(ValueError): + with np.nditer(estimateDepth(params, velocity_func(i)), flags=['multi_index']) as it: + for value in it: + p = np.array([it.multi_index[0], it.multi_index[1], value]) + positions.append(p) + + offset: int = mesh_count * dims[0] * dims[1] + for y in range(dims[1] - 1): + indices: list[int] = [] + for x in range(dims[0]): + indices.append(y * dims[0] + x + offset) + indices.append((y + 1) * dims[0] + x + offset) + mesh.addIndices(ivw.data.MeshInfo(dt=ivw.data.DrawType.Triangles, + ct=ivw.data.ConnectivityType.Strip), + ivw.data.IndexBufferUINT32(np.array(indices, dtype=np.uint32))) + mesh_count += 1 + + if len(positions) == 0: + return (mesh, np.array([0, 0])) + + positions_np = np.array(positions) + min_max = [np.min(positions_np[:, 2]), np.max(positions_np[:, 2])] + ivw.logInfo(f'min/max: {min_max}') + depth_scaling = np.max(np.abs(min_max)) + # scaling = np.array([1.0 / dims[0], 1.0 / dims[1], 1.0 / depth_scaling]) + + colors = [tf.sample(value / depth_scaling * 0.5 + 0.5) for value in positions_np[:, 2]] + + extent, offset = getLatLongBasis(params, lat_long) + ivw.logInfo(f'{extent=}\n{offset=}') + positions_np = (positions_np * np.array([extent[0] / (dims[0] - 1), + extent[1] / (dims[1] - 1), 1.0]) + - np.array([offset[0], offset[1], 0.0])) + + mesh.addBuffer(ivw.data.BufferType.PositionAttrib, + ivw.data.Buffer(np.array(positions_np, dtype=np.float32))) + mesh.addBuffer(ivw.data.BufferType.ColorAttrib, + ivw.data.Buffer(np.array(colors, dtype=np.float32))) + + # m = ivw.glm.mat4([[1.0 / dims[0], 0.0, 0.0, 0.0], + # [0.0, 1.0 / dims[1], 0.0, 0.0], + # [0.0, 0.0, 1.0 / depth_scaling, 0.0], + # [-0.5, -0.5, 0.0, 1.0]]) + + m = ivw.glm.mat4(1.0) + + w = ivw.glm.mat4(1.0) + w[2][2] = 1.0 / depth_scaling * 20.0 + # w[3] = [-0.5, -0.5, 0.0, 1.0] + + mesh.modelMatrix = m + mesh.worldMatrix = w + + return (mesh, np.array(min_max)) + + +def createFitsCompositeProperty(identifier: str, + displayName: str) -> ivw.properties.CompositeProperty: + cb = ivw.properties.ConstraintBehavior + inv = ivw.properties.InvalidationLevel + semantics = ivw.properties.PropertySemantics + + assumedDistance = ivw.properties.FloatProperty("assumedDistance", "Assumed Distance [pc]", + ivw.md2doc("Assumed distance of the star in parsec [pc]."), # noqa e501 + 123.0, increment=0.001, + min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), + semantics=semantics.Text) + velocityInf = ivw.properties.FloatProperty("velocityInf", "v_infinity [km/s]", + ivw.md2doc( + "The hyperbolic excess speed v_inf [km/s]."), + 14.5, increment=0.001, + min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), + semantics=semantics.Text) + velocityStar = ivw.properties.FloatProperty("velocityStar", "v_star [km/s]", + ivw.md2doc("v_star [km/s]"), + -26.5, increment=0.001, + min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), + semantics=semantics.Text) + crpix3 = ivw.properties.FloatProperty("crpix3", "crpix3", + ivw.md2doc("Reference pixel in axis3"), + 0.0, increment=0.001, + min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), + semantics=semantics.Text, + invalidationLevel=inv.Valid) + crval3 = ivw.properties.FloatProperty("crval3", "crval3", + ivw.md2doc("Physical value of the reference pixel"), + 0.0, increment=0.001, + min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), + semantics=semantics.Text, + invalidationLevel=inv.Valid) + cdelt3 = ivw.properties.FloatProperty("cdelt3", "cdelt3", + ivw.md2doc(""), + 0.0, increment=0.001, + min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), + semantics=semantics.Text, + invalidationLevel=inv.Valid) + restFrequency = ivw.properties.FloatProperty("restFrequency", "Rest Frequency [Hz", + ivw.md2doc("restfreq [Hz]"), + 0.0, increment=0.001, + min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), + semantics=semantics.Text, + invalidationLevel=inv.Valid) + extent = ivw.properties.FloatVec3Property("dataExtent", "Extent [au]", + ivw.md2doc("Extent of the dataset in au."), + glm.vec3(0.0), increment=glm.vec3(0.001), + min=(glm.vec3(-100.0), cb.Ignore), + max=(glm.vec3(100.0), cb.Ignore), + semantics=semantics.Text, + invalidationLevel=inv.Valid) + offset = ivw.properties.FloatVec3Property("dataOffset", "Offset [au]", + ivw.md2doc("Coordinate of the lower left corner of the dataset in au."), # noqa e501 + glm.vec3(0.0), increment=glm.vec3(0.001), + min=(glm.vec3(-100.0), cb.Ignore), + max=(glm.vec3(100.0), cb.Ignore), + semantics=semantics.Text, + invalidationLevel=inv.Valid) + dataRange = ivw.properties.DoubleMinMaxProperty("dataRange", + "Data Range", + 0.0, 1.0, -1.70e308, 1.79e308, + increment=0.0001, + semantics=semantics.Text, + invalidationLevel=inv.Valid) + + for p in [crpix3, crval3, cdelt3, restFrequency, extent, offset, dataRange]: + p.readOnly = True + + prop = ivw.properties.CompositeProperty(identifier, displayName) + prop.addProperties([assumedDistance, velocityInf, velocityStar, + crpix3, crval3, cdelt3, restFrequency, + extent, offset, dataRange]) + return prop + + +class FitsVolumeSource(ivw.Processor): + """ + Reads a FITS file and outputs the first array as a Volume + """ + + def __init__(self, id, name): + ivw.Processor.__init__(self, id, name) + self.deserializing = False + + self.outport = ivw.data.VolumeOutport("Intensity") + self.addOutport(self.outport) + + self.depth = ivw.data.MeshOutport("depth") + self.addOutport(self.depth) + + self.filename = ivw.properties.FileProperty( + "filename", "Filename", "", "fits") + self.filename.addNameFilter(ivw.properties.FileExtension("fits", "FITS file format")) + + options = [ivw.properties.IntOption("absolute", "Absolute (degree)", LatLong.AbsoluteDeg.value), # noqa e501 + ivw.properties.IntOption("relativeDeg", "Relative (degree)", LatLong.RelativeDeg.value), # noqa e501 + ivw.properties.IntOption("relativeArcsec", "Relative (arcsec)", LatLong.RelativeArcsec.value)] # noqa e501 + self.latLonCoords = ivw.properties.OptionPropertyInt( + "latLonCoords", "Lat/Long Coordinates", + help=ivw.unindentMd2doc(r'''Defines how Latitude and Longitude coordinates are shown. + +* _Absolute_ full coordinate in decimal degrees +* _Relative (degree)_ coordinates in arcseconds relative to the star/center +* _Relative (arcsec)_ coordinates in arcseconds relative to the star/center'''), + options=options, selectedIndex=2) + + self.fitsParameters = createFitsCompositeProperty("fitsParameters", "Fits Parameters") + + self.basis = ivwbase.properties.BasisProperty("basis", "Basis and offset") + self.information = ivwbase.properties.VolumeInformationProperty( + "information", "Data information") + + self.tf = ivw.properties.TransferFunctionProperty("colormap", "Mesh Colormap", + ivw.doc.Document("Colormap used for coloring the frequency slices according to depth")) # noqa e501 + self.addProperties([self.filename, self.fitsParameters, self.latLonCoords, self.tf, + self.basis, self.information]) + + @staticmethod + def processorInfo(): + return ivw.ProcessorInfo( + classIdentifier="org.inviwo.FitsVolumeSource", + displayName="FitsVolumeSource", + category="Python", + codeState=ivw.CodeState.Stable, + tags=ivw.Tags.PY, + help=ivw.unindentMd2doc(FitsVolumeSource.__doc__) + ) + + def getProcessorInfo(self): + return FitsVolumeSource.processorInfo() + + def initializeResources(self): + pass + + @staticmethod + def createVolume(data: np.array) -> ivw.data.Volume: + # data = np.nan_to_num(data) + volumerep = ivw.data.VolumePy(np.copy(data.astype(dtype=np.float32), order='C')) + + volume = ivw.data.Volume(volumerep) + volume.swizzlemask = ivw.data.SwizzleMask.defaultData(1) + + # TODO: extract extent from hdul[0].header + volume.basis = ivw.glm.mat3(1) + volume.offset = ivw.glm.vec3(-0.5, -0.5, -0.5) + + datarange = ivw.glm.dvec2(np.nanmin(data), np.nanmax(data)) + volume.dataMap.dataRange = datarange + volume.dataMap.valueRange = datarange + return volume + + @staticmethod + def frequencyToVelocityFunc(params: FitsParams) -> Callable[[float], int]: + def func(channel: int) -> float: + freq: int = (params.rest_frequency + - ((channel + 1 - params.crpix[2]) * params.cdelt[2] + params.crval[2])) + return vc * freq / params.rest_frequency / 1000 + + vc = 299_792_458 # speed of light + return func + + def process(self): + if not self.filename.value: + return + + filename: Path = ivwbase.io.downloadAndCacheIfUrl(self.filename.valueAsString()) + if not filename.exists(): + return + + with fits.open(filename) as hdul: + dim = hdul[0].data.shape + + fits_parameters: FitsParams = extractFitsParams(hdul[0].header) + + func: Callable[[float], int] = FitsVolumeSource.frequencyToVelocityFunc(fits_parameters) + velocity: list = [func(channel) for channel in range(dim[0])] + print(f'Velocity matching frequency: {velocity}') + + depth_mesh, min_max = \ + createDepthMesh(fits_parameters, func, self.tf.value, + LatLong(self.latLonCoords.value), + assumed_distance=self.fitsParameters.assumedDistance.value, + velocity_inf=self.fitsParameters.velocityInf.value, + velocity_star=self.fitsParameters.velocityStar.value) + + volume_intensity = FitsVolumeSource.createVolume(hdul[0].data) + volume_intensity.dataMap.valueAxis.name = 'Intensity' + # hdul[0].header['bunit'] is 'Jy/beam ', which is not a valid unit in Inviwo + volume_intensity.dataMap.valueAxis.unit = ivw.data.Unit( + fits_parameters.bunit.replace('/beam', '')) + + extent, offset = getLatLongBasis(fits_parameters, LatLong(self.latLonCoords.value)) + + m = ivw.glm.mat4(1.0) + m[0][0] = extent[0] + m[1][1] = extent[1] + m[2][2] = np.abs(velocity[-1] - velocity[0]) + m[3] = ivw.glm.vec4(offset[0], offset[1], velocity[0], 1.0) + + axes = [("x", "\""), ("y", "\""), ("Velocity", "km/s")] + for i, (label, unit) in enumerate(axes): + volume_intensity.axes[i].name = label + volume_intensity.axes[i].unit = ivw.data.Unit(unit) + + w = ivw.glm.mat4(1.0) + w[2][2] = 2.0 / m[2][2] + # w[3] = [-0.5, -0.5, 0.0, 1.0] + volume_intensity.modelMatrix = m + volume_intensity.worldMatrix = w + + self.fitsParameters.crpix3.value = fits_parameters.crpix[2] + self.fitsParameters.crval3.value = fits_parameters.crval[2] + self.fitsParameters.cdelt3.value = fits_parameters.cdelt[2] + self.fitsParameters.restFrequency.value = fits_parameters.rest_frequency + self.fitsParameters.dataRange.value = volume_intensity.dataMap.dataRange + + extent, offset = getLatLongBasis(fits_parameters, LatLong.RelativeArcsec) + extent *= self.fitsParameters.assumedDistance.value + offset *= self.fitsParameters.assumedDistance.value + self.fitsParameters.dataExtent.value = glm.vec3(extent[0], extent[1], + min_max[1] - min_max[0]) + self.fitsParameters.dataOffset.value = glm.vec3(offset[0], offset[1], min_max[0]) + + from inviwopy.properties import OverwriteState as ows + + self.basis.updateForNewEntity(volume_intensity, self.deserializing) + self.information.updateForNewVolume( + volume_intensity, ows.Yes if self.deserializing else ows.No) + self.basis.updateEntity(volume_intensity) + self.information.updateVolume(volume_intensity) + self.deserializing = False + + self.outport.setData(volume_intensity) + self.depth.setData(depth_mesh) + + def deserialize(self, deserializer: ivw.Deserializer): + ivw.Processor.deserialize(self, deserializer) + self.deserializing = True diff --git a/modules/astrophysics/readme.md b/modules/astrophysics/readme.md new file mode 100644 index 000000000..5cd47b350 --- /dev/null +++ b/modules/astrophysics/readme.md @@ -0,0 +1,3 @@ +# AstroPhysics Module + +Description of the AstroPhysics module diff --git a/modules/astrophysics/src/astrophysicsmodule.cpp b/modules/astrophysics/src/astrophysicsmodule.cpp new file mode 100644 index 000000000..9f405084b --- /dev/null +++ b/modules/astrophysics/src/astrophysicsmodule.cpp @@ -0,0 +1,52 @@ +/********************************************************************************* + * + * Inviwo - Interactive Visualization Workshop + * + * Copyright (c) 2026 Inviwo Foundation + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + *********************************************************************************/ + +#include +#include + +#include + +#include +#include + +namespace inviwo { + +AstroPhysicsModule::AstroPhysicsModule(InviwoApplication* app) try + : InviwoModule{app, "AstroPhysics"} + , scripts_{InviwoModule::getPath() / "python"} + , pythonFolderObserver_{app, InviwoModule::getPath() / "python/processors", *this} { + + // ShaderManager::getPtr()->addShaderSearchPath(getPath(ModulePath::GLSL)); + + registerProcessor(); +} catch (const std::exception& e) { + throw ModuleInitException(e.what()); +} + +} // namespace inviwo diff --git a/modules/astrophysics/src/processors/fitsvolumeraycaster.cpp b/modules/astrophysics/src/processors/fitsvolumeraycaster.cpp new file mode 100644 index 000000000..3b5cb464f --- /dev/null +++ b/modules/astrophysics/src/processors/fitsvolumeraycaster.cpp @@ -0,0 +1,72 @@ +/********************************************************************************* + * + * Inviwo - Interactive Visualization Workshop + * + * Copyright (c) 2026 Inviwo Foundation + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + *********************************************************************************/ + +#include + +#include +#include + +namespace inviwo { + +// The Class Identifier has to be globally unique. Use a reverse DNS naming scheme +const ProcessorInfo FitsVolumeRaycaster::processorInfo_{ + "org.infravis.FitsVolumeRaycaster", // Class identifier + "Fits Volume Raycaster", // Display name + "Volume Rendering", // Category + CodeState::Experimental, // Code state + Tags::GL | Tag{"Volume"} | Tag{"Raycaster"} | Tag{"Fits"}, // Tags + R"(Visualizing volumetric astrophysical observation data by also taking non-linear + depth into consideration.)"_unindentHelp, +}; + +const ProcessorInfo& FitsVolumeRaycaster::getProcessorInfo() const { return processorInfo_; } + +FitsVolumeRaycaster::FitsVolumeRaycaster(std::string_view identifier, std::string_view displayName) + : VolumeRaycasterBase{identifier, displayName} + //, volume_{"volume", VolumeComponent::Gradients::Single, + // "input volume (Only one channel will be rendered)"_help} + , fits_{"volume"} + , entryExit_{} + , background_{*this} + , isoTF_{fits_.volumePort_} + , raycasting_{fits_.getName(), isoTF_.isotfs[0]} + , camera_{"camera", util::boundingBox(fits_.volumePort_)} + , light_{&camera_.camera} +{ + + registerComponents(fits_, entryExit_, background_, raycasting_, isoTF_, camera_, + light_); +} + +void FitsVolumeRaycaster::process() { + util::checkValidChannel(raycasting_.selectedChannel(), fits_.channelsForVolume().value_or(0)); + VolumeRaycasterBase::process(); +} + +} // namespace inviwo diff --git a/modules/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp b/modules/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp new file mode 100644 index 000000000..edc86c0b2 --- /dev/null +++ b/modules/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp @@ -0,0 +1,303 @@ +/********************************************************************************* + * + * Inviwo - Interactive Visualization Workshop + * + * Copyright (c) 2026 Inviwo Foundation + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + *********************************************************************************/ + +#include + +#include +#include +#include + +#include +#include +#include + +namespace inviwo { + +FitsNonlinearDepthComponent::FitsNonlinearDepthComponent(std::string_view name, Document help) + : ShaderComponent{} + , volumePort_{name, std::move(help)} + , fitsParameters_{"fitsParameters", "FitsParameters", + "Various parameters required for coordinate transformations and other calculations."_help} + , assumedDistance_{"assumedDistance", "Assumed Distance [pc]", + util::ordinalSymmetricVector(0.0f) + .setInc(0.001f) + .set(PropertySemantics::Text) + .set("Assumed distance of the star in parsec [pc]."_help)} + , velocityInf_{"velocityInf", "v_infinity [km/s]", + util::ordinalSymmetricVector(0.0f) + .setInc(0.001f) + .set(PropertySemantics::Text) + .set("The hyperbolic excess speed v_inf [km/s]."_help)} + , velocityStar_{"velocityStar", "v_star [km/s]", + util::ordinalSymmetricVector(0.0f) + .setInc(0.001f) + .set(PropertySemantics::Text) + .set("v_star [km/s]"_help)} + , crpix3_{"crpix3", "crpix3", + util::ordinalSymmetricVector(0.0f) + .setInc(0.001f) + .set(PropertySemantics::Text) + .set("Reference pixel in axis3"_help)} + , crval3_{"crval3", "crval3", + util::ordinalSymmetricVector(0.0f) + .setInc(0.001f) + .set(PropertySemantics::Text) + .set("Physical value of the reference pixel"_help)} + , cdelt3_{"cdelt3", "cdelt3", + util::ordinalSymmetricVector(0.0f) + .setInc(0.001f) + .set(PropertySemantics::Text) + .set(""_help)} + , restFrequency_{"restFrequency", "Rest Frequency [Hz]", + util::ordinalSymmetricVector(0.0f) + .setInc(0.001f) + .set(PropertySemantics::Text) + .set("restfreq [Hz]"_help)} + , dataExtent_{"dataExtent", "Extent [au]", + util::ordinalSymmetricVector(vec3{0.0f}) + .setInc(vec3{0.001f}) + .set(PropertySemantics::Text) + .set("Extent of the dataset in au."_help)} + , dataOffset_{"dataOffset", "Offset [au]", + util::ordinalSymmetricVector(vec3{0.0f}) + .setInc(vec3{0.001f}) + .set(PropertySemantics::Text) + .set("Coordinate of the lower left corner of the dataset in au"_help)} { + + fitsParameters_.addProperties(assumedDistance_, velocityInf_, velocityStar_, crpix3_, crval3_, + cdelt3_, restFrequency_, dataExtent_, dataOffset_); + fitsParameters_.setCollapsed(true); +} + +std::string_view FitsNonlinearDepthComponent::getName() const { + return volumePort_.getIdentifier(); +} + +void FitsNonlinearDepthComponent::process(Shader& shader, TextureUnitContainer& cont) { + utilgl::bindAndSetUniforms(shader, cont, volumePort_); + + const auto data = volumePort_.getData(); + const vec3 delta = dataExtent_.get() / vec3{data->getDimensions()}; + mat4 m = glm::scale(delta * vec3{data->getDimensions() + size3_t{1u}}); + m[3] = vec4{dataOffset_.get() - delta * 0.5f, 1.0f}; + + StrBuffer buf{"{}FitsParameters", getName()}; + std::string prefix{buf}; + shader.setUniform(buf.replace("{}.textureToModel", prefix), m); + shader.setUniform(buf.replace("{}.assumedDistance", prefix), assumedDistance_); + shader.setUniform(buf.replace("{}.velocityInf", prefix), velocityInf_); + shader.setUniform(buf.replace("{}.velocityStar", prefix), velocityStar_); + shader.setUniform(buf.replace("{}.crpix3", prefix), crpix3_); + shader.setUniform(buf.replace("{}.crval3", prefix), crval3_); + shader.setUniform(buf.replace("{}.cdelt3", prefix), cdelt3_); + shader.setUniform(buf.replace("{}.restFrequency", prefix), restFrequency_); +} + +std::vector> FitsNonlinearDepthComponent::getInports() { + return {{&volumePort_, std::string{"volumes"}}}; +} + +void FitsNonlinearDepthComponent::initializeResources(Shader&) { + // shader definitions... +} + +std::vector FitsNonlinearDepthComponent::getProperties() { return {&fitsParameters_}; } + +namespace { + +constexpr std::string_view includes = util::trim(R"( +struct FitsParameters { + mat4 textureToModel; // transform from [0 1] to relative distances [au] + + float assumedDistance; // [parsec] + float velocityInf; // [km/s] + float velocityStar; // [km/s] + + float crpix3; + float crval3; + float cdelt3; + float restFrequency; // [Hz] +}; + +const float invalidValue = -1.e20; +)"); + +constexpr std::string_view uniforms = util::trim(R"( +uniform VolumeParameters {0}Parameters; +uniform sampler3D {0}; + +uniform FitsParameters {0}FitsParameters; +)"); + +constexpr std::string_view utilityFuncs = util::trim(R"( +// Convert a velocity [km/s] to frequency channel [-]. +float velocityToChannel(in FitsParameters params, in float velocity) {{ + float vc = 299792458; // speed of light [m/s] + return (((params.restFrequency * + (1.0 - velocity * 1000 / vc) - params.crval3) / params.cdelt3 + params.crpix3 - 1)); +}} + +// Convert a relative spatial coordinate [au] to velocity [km/s] +float spatialCoordinateToVelocity(in FitsParameters params, in vec3 coord) {{ + // convert xy coord from arcseconds to au + //coord.xy *= params.assumedDistance; + return params.velocityInf * coord.z / length(coord) + params.velocityStar; +}} + +// TODO: function from normalized texture coords to normalized, non-linear frequency domain +// using velocityToChannel(spatialCoordinateToVelocity(...)) +// requires extent in au + +// Sample the volume texture by transforming normalized texture coordinates @p texCoord into +// the non-linear frequency domain using velocityToChannel(spatialCoordinateToVelocity(...)). +// +// @param texCoord normalized texture coordinate [0, 1] +// @return volume sample at re-normalized texture coordinate, if v_los > v_star. +// Otherwise vec4(invalidValue); +vec4 getNormalizedVoxel(in FitsParameters params, in vec3 texCoord) {{ + vec3 coord = (params.textureToModel * vec4(texCoord, 1.0)).xyz; + float velocity = spatialCoordinateToVelocity(params, coord); + + float v_los = velocity - params.velocityStar; + if (abs(v_los) > params.velocityInf) {{ + return vec4(invalidValue); + }} + texCoord = vec3(texCoord.xy, velocityToChannel(params, velocity) / textureSize(volume, 0).z); + return getNormalizedVoxel({0}, {0}Parameters, texCoord); +}} +)"); + +// Initialize the VoxelPrev value to the same as the first voxel value. This value is important +// mainly for the isosurface rendering. Setting it to the same voxel value prevents isosurfaces +// being rendered at the volume boundaries. +constexpr std::string_view voxelFirst = util::trim(R"( +vec4 {0}Voxel = getNormalizedVoxel({0}FitsParameters, samplePosition); +vec4 {0}VoxelPrev = {0}Voxel; +)"); + +constexpr std::string_view voxel = util::trim(R"( +{0}VoxelPrev = {0}Voxel; +{0}Voxel = getNormalizedVoxel({0}FitsParameters, samplePosition); +)"); + +constexpr std::string_view gradientFirst = util::trim(R"( +vec3 {0}GradientPrev = vec3(0); +vec3 {0}Gradient = vec3(0); +#if defined(GRADIENTS_ENABLED) +{0}Gradient = useSurfaceNormals ? -texture(surfaceNormal, texCoords).xyz : + normalize(COMPUTE_GRADIENT_FOR_CHANNEL({0}Voxel, {0}, {0}Parameters, + samplePosition, channel)); +if (!useSurfaceNormals) {{ + {0}Gradient *= sign({0}Voxel[channel] / {0}Parameters.texToNormalized.scale + {0}Parameters.texToNormalized.offset); +}} +#endif +)"); + +constexpr std::string_view gradient = util::trim(R"( +#if defined(GRADIENTS_ENABLED) +{0}GradientPrev = {0}Gradient; +{0}Gradient = normalize(COMPUTE_GRADIENT_FOR_CHANNEL({0}Voxel, {0}, {0}Parameters, + samplePosition, channel)); +{0}Gradient *= sign({0}Voxel[channel] / {0}Parameters.texToNormalized.scale + {0}Parameters.texToNormalized.offset); +#endif +)"); + +constexpr std::string_view allGradientsFirst = util::trim(R"( +mat4x3 {0}AllGradientsPrev = mat4x3(0); +mat4x3 {0}AllGradients = mat4x3(0); +#if defined(GRADIENTS_ENABLED) +vec3 surfaceNormal = useSurfaceNormals ? -texture(surfaceNormal, texCoords).xyz : vec3(0); +{0}AllGradients = useSurfaceNormals ? + mat4x3(surfaceNormal, surfaceNormal, surfaceNormal, surfaceNormal) : + COMPUTE_ALL_GRADIENTS({0}Voxel, {0}, {0}Parameters, samplePosition); +{0}AllGradients[0] = normalize({0}AllGradients[0]); +{0}AllGradients[1] = normalize({0}AllGradients[1]); +{0}AllGradients[2] = normalize({0}AllGradients[2]); +{0}AllGradients[3] = normalize({0}AllGradients[3]); +#endif +)"); + +constexpr std::string_view allGradients = util::trim(R"( +#if defined(GRADIENTS_ENABLED) +{0}AllGradientsPrev = {0}AllGradients; +{0}AllGradients = COMPUTE_ALL_GRADIENTS({0}Voxel, {0}, {0}Parameters, samplePosition); +{0}AllGradients[0] = normalize({0}AllGradients[0]); +{0}AllGradients[1] = normalize({0}AllGradients[1]); +{0}AllGradients[2] = normalize({0}AllGradients[2]); +{0}AllGradients[3] = normalize({0}AllGradients[3]); +#endif +)"); + +} // namespace + +auto FitsNonlinearDepthComponent::getSegments() -> std::vector { + using fmt::literals::operator""_a; + + std::vector segments{ + {.snippet = std::string{includes}, .placeholder = placeholder::include, .priority = 400}, + {.snippet = fmt::format(uniforms, getName()), + .placeholder = placeholder::uniform, + .priority = 400}, + {.snippet = fmt::format(utilityFuncs, getName()), + .placeholder = placeholder::uniform, + .priority = 800}, + {.snippet = fmt::format(voxelFirst, getName()), + .placeholder = placeholder::first, + .priority = 400}, + {.snippet = fmt::format(voxel, getName()), + .placeholder = placeholder::loop, + .priority = 400}, + }; + + // if (gradients != Gradients::None) { + // segments.emplace_back(std::string{R"(#include "utils/gradients.glsl")"}, + // placeholder::include, 400); + // } + // if (gradients == Gradients::Single) { + // segments.emplace_back(fmt::format(gradientFirst, getName()), placeholder::first, 410); + // segments.emplace_back(fmt::format(gradient, getName()), placeholder::loop, 410); + // } + + // if (gradients == Gradients::All) { + // segments.emplace_back(fmt::format(allGradientsFirst, getName()), placeholder::first, + // 410); segments.emplace_back(fmt::format(allGradients, getName()), placeholder::loop, + // 410); + // } + + return segments; +} + +std::optional FitsNonlinearDepthComponent::channelsForVolume() const { + if (auto data = volumePort_.getData()) { + return data->getDataFormat()->getComponents(); + } + return std::nullopt; +} + +} // namespace inviwo diff --git a/modules/astrophysics/tests/unittests/astrophysics-unittest-main.cpp b/modules/astrophysics/tests/unittests/astrophysics-unittest-main.cpp new file mode 100644 index 000000000..acb90bbf0 --- /dev/null +++ b/modules/astrophysics/tests/unittests/astrophysics-unittest-main.cpp @@ -0,0 +1,57 @@ +/********************************************************************************* + * + * Inviwo - Interactive Visualization Workshop + * + * Copyright (c) 2026 Inviwo Foundation + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + *********************************************************************************/ + +#ifdef _MSC_VER +#pragma comment(linker, "/SUBSYSTEM:CONSOLE") +#endif + +#include +#include +#include + +#include +#include +#include +#include + +int main(int argc, char** argv) { + using namespace inviwo; + LogCentral::init(); + auto logger = std::make_shared(); + LogCentral::getPtr()->setVerbosity(LogVerbosity::Error); + LogCentral::getPtr()->registerLogger(logger); + + int ret = -1; + { + ::testing::InitGoogleTest(&argc, argv); + inviwo::ConfigurableGTestEventListener::setup(); + ret = RUN_ALL_TESTS(); + } + return ret; +} From 3d0a0bdd999d47c68f520be66a5afb8bbcace2b9 Mon Sep 17 00:00:00 2001 From: Martin Falk Date: Fri, 27 Feb 2026 16:24:06 +0100 Subject: [PATCH 02/16] Astrophysics: Python refactoring --- modules/astrophysics/python/astroviscommon.py | 79 ++- modules/astrophysics/python/ivwastrovis.py | 174 ++++++ .../python/processors/AngleToVector2D.py | 67 ++ .../python/processors/AxisSliceAnnotation.py | 2 +- .../python/processors/FitsDepthCorrection.py | 74 +++ .../python/processors/FitsVolumeSource.py | 582 +++++------------- 6 files changed, 545 insertions(+), 433 deletions(-) create mode 100644 modules/astrophysics/python/ivwastrovis.py create mode 100644 modules/astrophysics/python/processors/AngleToVector2D.py create mode 100644 modules/astrophysics/python/processors/FitsDepthCorrection.py diff --git a/modules/astrophysics/python/astroviscommon.py b/modules/astrophysics/python/astroviscommon.py index 6bd91f0d1..e80101196 100644 --- a/modules/astrophysics/python/astroviscommon.py +++ b/modules/astrophysics/python/astroviscommon.py @@ -10,12 +10,25 @@ class LatLong(Enum): AbsoluteDeg = 0 - RelativeDeg = 1 + FractionalArcsec = 1 # discard degree and arcmin RelativeArcsec = 2 +def getBasisUnit(lat_long: LatLong) -> str: + match lat_long: + case LatLong.AbsoluteDeg: + return "°" + case LatLong.FractionalArcsec: + return "\"" + case LatLong.RelativeArcsec: + return "\"" + case _: + raise ValueError(f'Invalid LatLong value {lat_long}') + + @dataclass class FitsParams: + objectName: str naxis: tuple[int, ...] ctype: tuple[str, str, str] cunit: tuple[str, str, str] @@ -27,12 +40,30 @@ class FitsParams: rest_frequency: float # [Hz] +@dataclass +class FitsData: + params: FitsParams + data: np.ndarray + + +@lru_cache(maxsize=1) +def readFITSData(filepath: Path, + dataset_index: int | None = None, + **kwargs) -> FitsData: + if not dataset_index: + dataset_index = 0 + with fits.open(filepath, **kwargs) as hdul: + fits_parameters: FitsParams = extractFITSParams(hdul[dataset_index].header) + data: np.ndarray = hdul[dataset_index].data # TODO: copy data? + return FitsData(params=fits_parameters, data=data) + + def extractFITSParams(header: fits.header.Header) -> FitsParams: rest_frequency: float = header['restfrq'] # equals 345_339_760_000 - # mean difference 0.0080740896 # rest_frequency = 345_339_769_300 # [hz] return FitsParams( + objectName=str(header['object']), naxis=tuple(header[f'naxis{i}'] for i in range(1, header['naxis'] + 1)), ctype=tuple(header[f'ctype{i}'].lower() for i in range(1, 4)), cunit=tuple(header[f'cunit{i}'] for i in range(1, 4)), @@ -44,16 +75,22 @@ def extractFITSParams(header: fits.header.Header) -> FitsParams: rest_frequency=rest_frequency) -@lru_cache(maxsize=2) -def readFITSData(filepath: Path, - dataset_index: int | None = None, - **kwargs) -> tuple[FitsParams, np.ndarray]: - if not dataset_index: - dataset_index = 0 - with fits.open(filepath, **kwargs) as hdul: - fits_parameters: FitsParams = extractFITSParams(hdul[dataset_index].header) - data: np.ndarray = hdul[dataset_index].data # TODO: copy data? - return (fits_parameters, data) +def frequencyFunc(params: FitsParams) -> Callable[[float], int]: + def func(channel: int) -> float: + freq: float = (channel + 1 - params.crpix[2]) * params.cdelt[2] + params.crval[2] + return freq + + return func + + +def frequencyToVelocityFunc(params: FitsParams) -> Callable[[float], int]: + def func(channel: int) -> float: + freq: float = (params.rest_frequency + - ((channel + 1 - params.crpix[2]) * params.cdelt[2] + params.crval[2])) + return vc * freq / params.rest_frequency / 1000 + + vc = 299_792_458 # speed of light + return func def pixelCoordinatesFunc(params: FitsParams) -> Callable[[np.ndarray], [int, int]]: @@ -144,7 +181,7 @@ def estimateDepth(params: FitsParams, def getLatLongBasis(params: FitsParams, lat_long: LatLong = LatLong.RelativeArcsec - ) -> tuple[np.ndarray, np.ndarray]: + ) -> tuple[np.ndarray, np.ndarray, str]: """ Calculate the lat/long extent and offset given a fits header for pixel centers @@ -158,8 +195,8 @@ def getLatLongBasis(params: FitsParams, Returns ------- - tuple[np.ndarray, np.ndarray] - lat/long extent and offset + tuple[np.ndarray, np.ndarray, str] + lat/long extent and offset, and unit """ coords: Callable[[np.ndarray], [int, int]] = pixelCoordinatesFunc(params) @@ -173,14 +210,18 @@ def getLatLongBasis(params: FitsParams, extent[i] = -extent[i] delta[i] = -delta[i] - if lat_long == LatLong.RelativeArcsec: + if lat_long in [LatLong.RelativeArcsec, LatLong.FractionalArcsec]: extent *= 3600 # convert from degree to arcseconds + origin *= 3600 + delta *= 3600 if lat_long == LatLong.AbsoluteDeg: offset = origin - delta * 0.5 - elif lat_long in [LatLong.RelativeDeg, LatLong.RelativeArcsec]: - offset = extent * 0.5 + elif lat_long == LatLong.FractionalArcsec: + offset = np.fmod(origin - delta * 0.5, 60) + elif lat_long == LatLong.RelativeArcsec: + offset = -extent * 0.5 else: raise ValueError(f'Invalid Lat/Long value {lat_long}') - return (extent, offset) + return (extent, offset, getBasisUnit(lat_long)) diff --git a/modules/astrophysics/python/ivwastrovis.py b/modules/astrophysics/python/ivwastrovis.py new file mode 100644 index 000000000..90ad15cff --- /dev/null +++ b/modules/astrophysics/python/ivwastrovis.py @@ -0,0 +1,174 @@ +import astroviscommon +from astroviscommon import FitsParams, LatLong + +import inviwopy as ivw +import ivwbase +from inviwopy import glm + +import numpy as np +from contextlib import suppress +from collections.abc import Callable + + +def createDepthMesh(params: FitsParams, + tf: ivw.data.TransferFunction, + lat_long: LatLong, + assumed_distance: float = 123.0, + velocity_inf: float = 14.5, + velocity_star: float = -26.5 + ) -> (ivw.data.Mesh, np.ndarray): + + velocity_func: Callable[[float], int] = astroviscommon.frequencyToVelocityFunc(params) + + dims = np.array(params.naxis) + + positions: list[np.ndarray] = [] + + mesh_count: int = 0 + mesh = ivw.data.Mesh(dt=ivw.data.DrawType.Triangles, ct=ivw.data.ConnectivityType.Unconnected) + for i in range(0, dims[2]): + with suppress(ValueError): + with np.nditer(astroviscommon.estimateDepth(params, velocity_func(i)), + flags=['multi_index']) as it: + for value in it: + p = np.array([it.multi_index[0], it.multi_index[1], value]) + positions.append(p) + + offset: int = mesh_count * dims[0] * dims[1] + for y in range(dims[1] - 1): + indices: list[int] = [] + for x in range(dims[0]): + indices.append(y * dims[0] + x + offset) + indices.append((y + 1) * dims[0] + x + offset) + mesh.addIndices(ivw.data.MeshInfo(dt=ivw.data.DrawType.Triangles, + ct=ivw.data.ConnectivityType.Strip), + ivw.data.IndexBufferUINT32(np.array(indices, dtype=np.uint32))) + mesh_count += 1 + + if len(positions) == 0: + return (mesh, np.array([0, 0])) + + positions_np = np.array(positions) + min_max = [np.min(positions_np[:, 2]), np.max(positions_np[:, 2])] + ivw.logInfo(f'min/max: {min_max}') + depth_scaling = np.max(np.abs(min_max)) + # scaling = np.array([1.0 / dims[0], 1.0 / dims[1], 1.0 / depth_scaling]) + + colors = [tf.sample(value / depth_scaling * 0.5 + 0.5) for value in positions_np[:, 2]] + + extent, offset = astroviscommon.getLatLongBasis(params, lat_long) + ivw.logInfo(f'{extent=}\n{offset=}') + positions_np = (positions_np * np.array([extent[0] / (dims[0] - 1), + extent[1] / (dims[1] - 1), 1.0]) + - np.array([offset[0], offset[1], 0.0])) + + mesh.addBuffer(ivw.data.BufferType.PositionAttrib, + ivw.data.Buffer(np.array(positions_np, dtype=np.float32))) + mesh.addBuffer(ivw.data.BufferType.ColorAttrib, + ivw.data.Buffer(np.array(colors, dtype=np.float32))) + + # m = ivw.glm.mat4([[1.0 / dims[0], 0.0, 0.0, 0.0], + # [0.0, 1.0 / dims[1], 0.0, 0.0], + # [0.0, 0.0, 1.0 / depth_scaling, 0.0], + # [-0.5, -0.5, 0.0, 1.0]]) + + m = ivw.glm.mat4(1.0) + + w = ivw.glm.mat4(1.0) + w[2][2] = 1.0 / depth_scaling * 20.0 + # w[3] = [-0.5, -0.5, 0.0, 1.0] + + mesh.modelMatrix = m + mesh.worldMatrix = w + + return (mesh, np.array(min_max)) + + +def createFitsCompositeProperty(identifier: str, + displayName: str) -> ivw.properties.CompositeProperty: + cb = ivw.properties.ConstraintBehavior + inv = ivw.properties.InvalidationLevel + semantics = ivw.properties.PropertySemantics + + objectName = ivw.properties.StringProperty( + "objectName", "Object", + ivw.md2doc( + "Contents of the 'objects' field in the FITS Header."), "", + invalidationLevel=inv.Valid) + assumedDistance = ivw.properties.FloatProperty("assumedDistance", "Assumed Distance [pc]", + ivw.md2doc("Assumed distance of the star in parsec [pc]."), # noqa e501 + 123.0, increment=0.001, + min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), + semantics=semantics.Text) + velocityInf = ivw.properties.FloatProperty("velocityInf", "v_infinity [km/s]", + ivw.md2doc( + "The hyperbolic excess speed v_inf [km/s]."), + 14.5, increment=0.001, + min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), + semantics=semantics.Text) + velocityStar = ivw.properties.FloatProperty("velocityStar", "v_star [km/s]", + ivw.md2doc("v_star [km/s]"), + -26.5, increment=0.001, + min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), + semantics=semantics.Text) + crpix3 = ivw.properties.FloatProperty("crpix3", "crpix3", + ivw.md2doc("Reference pixel in axis3"), + 0.0, increment=0.001, + min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), + semantics=semantics.Text, + invalidationLevel=inv.Valid) + crval3 = ivw.properties.FloatProperty("crval3", "crval3", + ivw.md2doc("Physical value of the reference pixel"), + 0.0, increment=0.001, + min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), + semantics=semantics.Text, + invalidationLevel=inv.Valid) + cdelt3 = ivw.properties.FloatProperty("cdelt3", "cdelt3", + ivw.md2doc(""), + 0.0, increment=0.001, + min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), + semantics=semantics.Text, + invalidationLevel=inv.Valid) + restFrequency = ivw.properties.FloatProperty("restFrequency", "Rest Frequency [Hz", + ivw.md2doc("restfreq [Hz]"), + 0.0, increment=0.001, + min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), + semantics=semantics.Text, + invalidationLevel=inv.Valid) + extent = ivw.properties.FloatVec3Property("dataExtent", "Extent [au]", + ivw.md2doc("Extent of the dataset in au."), + glm.vec3(0.0), increment=glm.vec3(0.001), + min=(glm.vec3(-100.0), cb.Ignore), + max=(glm.vec3(100.0), cb.Ignore), + semantics=semantics.Text, + invalidationLevel=inv.Valid) + offset = ivw.properties.FloatVec3Property("dataOffset", "Offset [au]", + ivw.md2doc("Coordinate of the lower left corner of the dataset in au."), # noqa e501 + glm.vec3(0.0), increment=glm.vec3(0.001), + min=(glm.vec3(-100.0), cb.Ignore), + max=(glm.vec3(100.0), cb.Ignore), + semantics=semantics.Text, + invalidationLevel=inv.Valid) + valueName = ivw.properties.StringProperty( + "valueName", "btype", + ivw.md2doc("Corresponds to the 'btype' field in the FITS Header."), "", + invalidationLevel=inv.Valid) + valueUnit = ivw.properties.StringProperty( + "valueUnit", "bunit", + ivw.md2doc("Corresponds to the 'bunit' field in the FITS Header."), "", + invalidationLevel=inv.Valid) + + dataRange = ivwbase.properties.DataRangeProperty("dataRange", "Data Range") + dataRange.help = ivw.md2doc("Min/max values of the dataset, ignoring NaN") + dataRange.getPropertyByIdentifier("valueRange", True).visible = False + dataRange.getPropertyByIdentifier("customValueRange", True).visible = False + + for p in [objectName, crpix3, crval3, cdelt3, restFrequency, extent, offset, + valueName, valueUnit]: + p.readOnly = True + + prop = ivw.properties.CompositeProperty(identifier, displayName) + prop.addProperties([objectName, assumedDistance, velocityInf, velocityStar, + crpix3, crval3, cdelt3, restFrequency, + extent, offset, valueName, valueUnit, dataRange]) + return prop diff --git a/modules/astrophysics/python/processors/AngleToVector2D.py b/modules/astrophysics/python/processors/AngleToVector2D.py new file mode 100644 index 000000000..2e22bed3d --- /dev/null +++ b/modules/astrophysics/python/processors/AngleToVector2D.py @@ -0,0 +1,67 @@ +# Name: AngleToVector2D + +import inviwopy as ivw +import numpy as np + + +class AngleToVector2D(ivw.Processor): + """ + Converts the angles of the input volume into vectors and scales these with an optional magnitude + """ + + def __init__(self, id, name): + ivw.Processor.__init__(self, id, name) + self.addInport(ivw.data.VolumeInport("angle")) + self.addInport(ivw.data.VolumeInport("magnitude")) + self.inports.magnitude.optional = True + + self.addOutport(ivw.data.VolumeOutport("outport")) + + @staticmethod + def processorInfo(): + return ivw.ProcessorInfo( + classIdentifier="org.inviwo.AngleToVector2D", + displayName="AngleToVector2D", + category="Python", + codeState=ivw.CodeState.Stable, + tags=ivw.Tags("PY, FITS, Astrophysics"), + help=ivw.unindentMd2doc(AngleToVector2D.__doc__) + ) + + def getProcessorInfo(self): + return AngleToVector2D.processorInfo() + + def initializeResources(self): + pass + + def process(self): + volume_angle = self.inports.angle.getData() + volume_magnitude = self.inports.magnitude.getData() + + use_degree: bool = True + # assume degree if anything other than radians is stated as unit + if volume_angle.dataMap.valueAxis.unit == ivw.data.Unit("rad"): + use_degree = False + + angle = volume_angle.data + magnitude = volume_magnitude.data + + if use_degree: + angle *= np.pi / 180.0 + data = np.stack([np.cos(angle + np.pi / 2.0) * magnitude, + np.sin(angle + np.pi / 2.0) * magnitude], axis=3) + volumerep = ivw.data.VolumePy(np.copy(data.astype(dtype=np.float32), order='C')) + volume = ivw.data.Volume(volumerep) + + max_len: float = np.nanmax(magnitude) + datarange = ivw.glm.dvec2(-max_len, max_len) + volume.dataMap.dataRange = datarange + volume.dataMap.valueRange = datarange + + volume.dataMap.valueAxis = volume_magnitude.dataMap.valueAxis + volume.modelMatrix = volume_magnitude.modelMatrix + volume.worldMatrix = volume_magnitude.worldMatrix + volume.axes = volume_magnitude.axes + volume.copyMetaDataFrom(volume_magnitude) + + self.outports.outport.setData(volume) diff --git a/modules/astrophysics/python/processors/AxisSliceAnnotation.py b/modules/astrophysics/python/processors/AxisSliceAnnotation.py index 02be38c3c..58aa9935f 100644 --- a/modules/astrophysics/python/processors/AxisSliceAnnotation.py +++ b/modules/astrophysics/python/processors/AxisSliceAnnotation.py @@ -72,7 +72,7 @@ def process(self): axis = self.sliceAxis.value value: float = (volume.basis[axis] * (self.slice.value - 1) - / volume.dimensions[axis] + volume.offset)[axis] + / (volume.dimensions[axis] - 1) + volume.offset)[axis] formatString = self.labelType.selectedValue \ if self.labelType.selectedValue != "custom" else self.customLabel.value diff --git a/modules/astrophysics/python/processors/FitsDepthCorrection.py b/modules/astrophysics/python/processors/FitsDepthCorrection.py new file mode 100644 index 000000000..6f518a5f3 --- /dev/null +++ b/modules/astrophysics/python/processors/FitsDepthCorrection.py @@ -0,0 +1,74 @@ +# Name: FitsDepthCorrection + +import inviwopy as ivw +from inviwopy import glm +import astroviscommon +from astroviscommon import LatLong +import ivwastrovis + +import time + + +class FitsDepthCorrection(ivw.Processor): + """ + Reconstructs the depth information given by a FITS image stack and a set of FitsParameters. + + The corresponding spatial extent in au is also made available via the extent and + offset properties. + """ + + def __init__(self, id, name): + ivw.Processor.__init__(self, id, name) + + self.addInport(ivw.PythonInport("fitsdata")) + self.addOutport(ivw.data.MeshOutport("depthmesh")) + + self.fitsParameters = ivwastrovis.createFitsCompositeProperty( + "fitsParameters", "Fits Parameters") + self.tf = ivw.properties.TransferFunctionProperty( + "colormap", "Mesh Colormap", + ivw.doc.Document("Colormap used for coloring the frequency slices according to depth")) + + self.addProperties([self.fitsParameters, self.tf]) + + @staticmethod + def processorInfo(): + return ivw.ProcessorInfo( + classIdentifier="org.inviwo.FitsDepthCorrection", + displayName="FitsDepthCorrection", + category="Python", + codeState=ivw.CodeState.Stable, + tags=ivw.Tags("PY, FITS, Astrophysics"), + help=ivw.unindentMd2doc(FitsDepthCorrection.__doc__) + ) + + def getProcessorInfo(self): + return FitsDepthCorrection.processorInfo() + + def initializeResources(self): + pass + + def process(self): + fits_data = self.inports.fitsdata.getData() + if not isinstance(fits_data, astroviscommon.FitsData): + raise ValueError(f"Invalid input data of type '{type(fits_data)}', " + f"expected '{astroviscommon.FitsData}'.") + + t = time.perf_counter_ns() + depth_mesh, min_max = ivwastrovis.createDepthMesh( + fits_data.params, self.tf.value, + LatLong.RelativeArcsec, + assumed_distance=self.fitsParameters.assumedDistance.value, + velocity_inf=self.fitsParameters.velocityInf.value, + velocity_star=self.fitsParameters.velocityStar.value) + elapsed = time.perf_counter_ns() - t + print(f'{elapsed/1000_000.0} ms') + + extent, offset = astroviscommon.getLatLongBasis(fits_data.params, LatLong.RelativeArcsec) + extent *= self.fitsParameters.assumedDistance.value + offset *= self.fitsParameters.assumedDistance.value + self.fitsParameters.dataExtent.value = glm.vec3( + extent[0], extent[1], min_max[1] - min_max[0]) + self.fitsParameters.dataOffset.value = glm.vec3(offset[0], offset[1], min_max[0]) + + self.outports.depthmesh.setData(depth_mesh) diff --git a/modules/astrophysics/python/processors/FitsVolumeSource.py b/modules/astrophysics/python/processors/FitsVolumeSource.py index b06558262..bf1cc46da 100644 --- a/modules/astrophysics/python/processors/FitsVolumeSource.py +++ b/modules/astrophysics/python/processors/FitsVolumeSource.py @@ -3,332 +3,22 @@ import inviwopy as ivw from inviwopy import glm import ivwbase +import astroviscommon +from astroviscommon import LatLong +import ivwastrovis import numpy as np from pathlib import Path from collections.abc import Callable from enum import Enum -from dataclasses import dataclass -from contextlib import suppress -from astropy.io import fits - - -class LatLong(Enum): - AbsoluteDeg = 0 - RelativeDeg = 1 - RelativeArcsec = 2 - - -@dataclass -class FitsParams: - naxis: tuple[int, ...] - ctype: tuple[str, str, str] - cunit: tuple[str, str, str] - crpix: tuple[float, float, float] - crval: tuple[float, float, float] - cdelt: tuple[float, float, float] - btype: str - bunit: str - rest_frequency: float # [Hz] - - -def extractFitsParams(header: fits.header.Header) -> FitsParams: - rest_frequency: float = header['restfrq'] # equals 345_339_760_000 - # mean difference 0.0080740896 - # rest_frequency = 345_339_769_300 # [hz] - - return FitsParams( - naxis=tuple(header[f'naxis{i}'] for i in range(1, header['naxis'] + 1)), - ctype=tuple(header[f'ctype{i}'].lower() for i in range(1, 4)), - cunit=tuple(header[f'cunit{i}'] for i in range(1, 4)), - crpix=tuple(header[f'crpix{i}'] for i in range(1, 4)), - crval=tuple(header[f'crval{i}'] for i in range(1, 4)), - cdelt=tuple(header[f'cdelt{i}'] for i in range(1, 4)), - btype=header['btype'], - bunit=header['bunit'], - rest_frequency=rest_frequency) - - -def pixelCoordinatesFunc(params: FitsParams) -> Callable[[np.ndarray], [int, int]]: - def func(i: int, j: int) -> np.ndarray: - # FIXME: use i or (i+1) ? - alpha = params.crval[0] + params.cdelt[0] * (i - params.crpix[0]) - beta = params.crval[1] + params.cdelt[1] * (j - params.crpix[1]) - return np.array([alpha, beta]) - - if not params.ctype[0] == 'ra---sin': - raise ValueError( - f'Expected unit type "ctype1" to be "RA---SIN". Found "{params.ctype[0]}".') - if not params.ctype[1] == 'dec--sin': - raise ValueError( - f'Expected unit type "ctype2" to be "DEC--SIN". Found "{params.ctype[1]}".') - return func - - -""" -Conversion of velocity axis to 3rd coordinate: -- In this case via the assumption of a radial velocity law of the form: - |V - V*| = Vinf * (1-((Rw-r)/Rw)^2.5) for rRw -Vinf = 14.5 km/s -V* = -26.5 km/s -Rw = 54 au -(assumed distance 123 pc, so 1" = 123 au) - -Since in this particular case, we're at scales larger than Rw we're in the -simplest case where the velocity is constant V=Vinf=14.5 km/s. - -So the line of sight velocity Vlos = V - V* = Vinf * (z/r) where r=sqrt(x^2+y^2+z^2) - -taking rho = sqrt(x^2+y^2) and mu = Vlos/Vinf we have -Z = +- rho * mu/sqrt(1-mu^2) - -The positive and negative sign are for the blue-shifted (in front) and red-shifted -(behind the plane of the sky) parts of the envelope. -We should in this case disregard all velocity channels with |V-V*|=|Vlos|>Vinf. -And remember for the coordinate transformation from (relative to the center of -the emission) arcseconds to au is using 1"=123 au). -""" - - -def estimateDepth(params: FitsParams, - velocity: float, - assumed_distance: float = 123.0, - velocity_inf: float = 14.5, - velocity_star: float = -26.5) -> np.ndarray: - """ - - Parameters - ---------- - params : FitsParams - Parameters extracted from a fits.header.Header - velocity : float - in km/s. - assumed_distance : float, optional - in parsecs. The default is 123.0. - velocity_inf: float, optional - in km/s. The default is 14.5. - velocity_star: float, optional - in km/s. The default is -26.5. - - - Returns - ------- - 2D array of depth values for the given velocity - """ - v_los: float = velocity - velocity_star - if np.abs(v_los) > velocity_inf: - raise ValueError(f'v_los={v_los} exceeds v_inf={velocity_inf}') - - mu: float = v_los / velocity_inf - - coords: Callable[[np.ndarray], [int, int]] = pixelCoordinatesFunc(params) - # coordinates are given in degree, convert to arcseconds then parsecs (au) - delta: np.ndarray = (coords(1, 1) - coords(0, 0)) * \ - 3600.0 * assumed_distance - dims = np.array(params.naxis[:2]) - center = (dims - 1) / 2.0 - - depths = np.empty(dims, dtype=float) - with np.nditer(depths, flags=['multi_index'], op_flags=['writeonly']) as it: - for value in it: - d = (it.multi_index - center) * delta - rho: float = np.sqrt(d[0]**2 + d[1]**2) - value[...] = rho * mu / np.sqrt(1.0 - mu**2) - return depths - - -def getLatLongBasis(params: FitsParams, - lat_long: LatLong = LatLong.RelativeArcsec - ) -> tuple[np.ndarray, np.ndarray]: - """ - Calculate the lat/long extent and offset given a fits header for pixel centers +import time - Parameters - ---------- - params : FitsParams - Parameters extracted from a fits.header.Header - lat_long : LatLong, optional - Determines the coordinate system (absolute/relative, degree/arcsec, ...). - The default is LatLong.RelativeArcsec. - Returns - ------- - tuple[np.ndarray, np.ndarray] - lat/long extent and offset - - """ - coords: Callable[[np.ndarray], [int, int]] = pixelCoordinatesFunc(params) - origin: np.ndarray = coords(0, 0) - delta: np.ndarray = coords(1, 1) - origin - extent: np.ndarray = coords(params.naxis[0], params.naxis[1]) - origin - - for i in range(2): - if extent[i] < 0.0: - origin[i] += extent[i] - extent[i] = -extent[i] - delta[i] = -delta[i] - - if lat_long == LatLong.RelativeArcsec: - extent *= 3600 # convert from degree to arcseconds - - if lat_long == LatLong.AbsoluteDeg: - offset = origin - delta * 0.5 - elif lat_long in [LatLong.RelativeDeg, LatLong.RelativeArcsec]: - offset = extent * 0.5 - else: - raise ValueError(f'Invalid Lat/Long value {lat_long}') - - return (extent, offset) - - -def createDepthMesh(params: FitsParams, - velocity_func: Callable[[float], int], - tf: ivw.data.TransferFunction, - lat_long: LatLong, - assumed_distance: float = 123.0, - velocity_inf: float = 14.5, - velocity_star: float = -26.5 - ) -> (ivw.data.Mesh, np.ndarray): - dims = np.array(params.naxis) - - positions: list[np.ndarray] = [] - - mesh_count: int = 0 - mesh = ivw.data.Mesh(dt=ivw.data.DrawType.Triangles, ct=ivw.data.ConnectivityType.Unconnected) - for i in range(0, dims[2]): - with suppress(ValueError): - with np.nditer(estimateDepth(params, velocity_func(i)), flags=['multi_index']) as it: - for value in it: - p = np.array([it.multi_index[0], it.multi_index[1], value]) - positions.append(p) - - offset: int = mesh_count * dims[0] * dims[1] - for y in range(dims[1] - 1): - indices: list[int] = [] - for x in range(dims[0]): - indices.append(y * dims[0] + x + offset) - indices.append((y + 1) * dims[0] + x + offset) - mesh.addIndices(ivw.data.MeshInfo(dt=ivw.data.DrawType.Triangles, - ct=ivw.data.ConnectivityType.Strip), - ivw.data.IndexBufferUINT32(np.array(indices, dtype=np.uint32))) - mesh_count += 1 - - if len(positions) == 0: - return (mesh, np.array([0, 0])) - - positions_np = np.array(positions) - min_max = [np.min(positions_np[:, 2]), np.max(positions_np[:, 2])] - ivw.logInfo(f'min/max: {min_max}') - depth_scaling = np.max(np.abs(min_max)) - # scaling = np.array([1.0 / dims[0], 1.0 / dims[1], 1.0 / depth_scaling]) - - colors = [tf.sample(value / depth_scaling * 0.5 + 0.5) for value in positions_np[:, 2]] - - extent, offset = getLatLongBasis(params, lat_long) - ivw.logInfo(f'{extent=}\n{offset=}') - positions_np = (positions_np * np.array([extent[0] / (dims[0] - 1), - extent[1] / (dims[1] - 1), 1.0]) - - np.array([offset[0], offset[1], 0.0])) - - mesh.addBuffer(ivw.data.BufferType.PositionAttrib, - ivw.data.Buffer(np.array(positions_np, dtype=np.float32))) - mesh.addBuffer(ivw.data.BufferType.ColorAttrib, - ivw.data.Buffer(np.array(colors, dtype=np.float32))) - - # m = ivw.glm.mat4([[1.0 / dims[0], 0.0, 0.0, 0.0], - # [0.0, 1.0 / dims[1], 0.0, 0.0], - # [0.0, 0.0, 1.0 / depth_scaling, 0.0], - # [-0.5, -0.5, 0.0, 1.0]]) - - m = ivw.glm.mat4(1.0) - - w = ivw.glm.mat4(1.0) - w[2][2] = 1.0 / depth_scaling * 20.0 - # w[3] = [-0.5, -0.5, 0.0, 1.0] - - mesh.modelMatrix = m - mesh.worldMatrix = w - - return (mesh, np.array(min_max)) - - -def createFitsCompositeProperty(identifier: str, - displayName: str) -> ivw.properties.CompositeProperty: - cb = ivw.properties.ConstraintBehavior - inv = ivw.properties.InvalidationLevel - semantics = ivw.properties.PropertySemantics - - assumedDistance = ivw.properties.FloatProperty("assumedDistance", "Assumed Distance [pc]", - ivw.md2doc("Assumed distance of the star in parsec [pc]."), # noqa e501 - 123.0, increment=0.001, - min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), - semantics=semantics.Text) - velocityInf = ivw.properties.FloatProperty("velocityInf", "v_infinity [km/s]", - ivw.md2doc( - "The hyperbolic excess speed v_inf [km/s]."), - 14.5, increment=0.001, - min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), - semantics=semantics.Text) - velocityStar = ivw.properties.FloatProperty("velocityStar", "v_star [km/s]", - ivw.md2doc("v_star [km/s]"), - -26.5, increment=0.001, - min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), - semantics=semantics.Text) - crpix3 = ivw.properties.FloatProperty("crpix3", "crpix3", - ivw.md2doc("Reference pixel in axis3"), - 0.0, increment=0.001, - min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), - semantics=semantics.Text, - invalidationLevel=inv.Valid) - crval3 = ivw.properties.FloatProperty("crval3", "crval3", - ivw.md2doc("Physical value of the reference pixel"), - 0.0, increment=0.001, - min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), - semantics=semantics.Text, - invalidationLevel=inv.Valid) - cdelt3 = ivw.properties.FloatProperty("cdelt3", "cdelt3", - ivw.md2doc(""), - 0.0, increment=0.001, - min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), - semantics=semantics.Text, - invalidationLevel=inv.Valid) - restFrequency = ivw.properties.FloatProperty("restFrequency", "Rest Frequency [Hz", - ivw.md2doc("restfreq [Hz]"), - 0.0, increment=0.001, - min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), - semantics=semantics.Text, - invalidationLevel=inv.Valid) - extent = ivw.properties.FloatVec3Property("dataExtent", "Extent [au]", - ivw.md2doc("Extent of the dataset in au."), - glm.vec3(0.0), increment=glm.vec3(0.001), - min=(glm.vec3(-100.0), cb.Ignore), - max=(glm.vec3(100.0), cb.Ignore), - semantics=semantics.Text, - invalidationLevel=inv.Valid) - offset = ivw.properties.FloatVec3Property("dataOffset", "Offset [au]", - ivw.md2doc("Coordinate of the lower left corner of the dataset in au."), # noqa e501 - glm.vec3(0.0), increment=glm.vec3(0.001), - min=(glm.vec3(-100.0), cb.Ignore), - max=(glm.vec3(100.0), cb.Ignore), - semantics=semantics.Text, - invalidationLevel=inv.Valid) - dataRange = ivw.properties.DoubleMinMaxProperty("dataRange", - "Data Range", - 0.0, 1.0, -1.70e308, 1.79e308, - increment=0.0001, - semantics=semantics.Text, - invalidationLevel=inv.Valid) - - for p in [crpix3, crval3, cdelt3, restFrequency, extent, offset, dataRange]: - p.readOnly = True - - prop = ivw.properties.CompositeProperty(identifier, displayName) - prop.addProperties([assumedDistance, velocityInf, velocityStar, - crpix3, crval3, cdelt3, restFrequency, - extent, offset, dataRange]) - return prop +class AxisType(Enum): + SliceIndex = 0 + Dataset = 1 + Velocity = 2 class FitsVolumeSource(ivw.Processor): @@ -340,37 +30,71 @@ def __init__(self, id, name): ivw.Processor.__init__(self, id, name) self.deserializing = False - self.outport = ivw.data.VolumeOutport("Intensity") - self.addOutport(self.outport) + global astroviscommon + global LatLong + from importlib import reload + astroviscommon = reload(astroviscommon) + from astroviscommon import LatLong - self.depth = ivw.data.MeshOutport("depth") - self.addOutport(self.depth) + self.addOutport(ivw.data.VolumeOutport("intensity")) + self.addOutport(ivw.PythonOutport("fitsdata")) self.filename = ivw.properties.FileProperty( "filename", "Filename", "", "fits") self.filename.addNameFilter(ivw.properties.FileExtension("fits", "FITS file format")) - options = [ivw.properties.IntOption("absolute", "Absolute (degree)", LatLong.AbsoluteDeg.value), # noqa e501 - ivw.properties.IntOption("relativeDeg", "Relative (degree)", LatLong.RelativeDeg.value), # noqa e501 - ivw.properties.IntOption("relativeArcsec", "Relative (arcsec)", LatLong.RelativeArcsec.value)] # noqa e501 + cb = ivw.properties.ConstraintBehavior + IntOption = ivw.properties.IntOption + inv = ivw.properties.InvalidationLevel + semantics = ivw.properties.PropertySemantics + + options = [IntOption("absolute", "Absolute (degree)", LatLong.AbsoluteDeg.value), + IntOption("fractionalArcsec", "Fractional (arcsec)", + LatLong.FractionalArcsec.value), + IntOption("relativeArcsec", "Relative (arcsec)", LatLong.RelativeArcsec.value)] self.latLonCoords = ivw.properties.OptionPropertyInt( "latLonCoords", "Lat/Long Coordinates", help=ivw.unindentMd2doc(r'''Defines how Latitude and Longitude coordinates are shown. * _Absolute_ full coordinate in decimal degrees -* _Relative (degree)_ coordinates in arcseconds relative to the star/center +* _Relative (arcsec)_ coordinates in arcseconds relative to the star/center * _Relative (arcsec)_ coordinates in arcseconds relative to the star/center'''), options=options, selectedIndex=2) - self.fitsParameters = createFitsCompositeProperty("fitsParameters", "Fits Parameters") + options = [IntOption("slice", "Slice Indices", AxisType.SliceIndex.value), + IntOption("dataset", "From DataSet", AxisType.Dataset.value), + IntOption("velocity", "Velocity [km/s]", AxisType.Velocity.value)] + self.axisType = ivw.properties.OptionPropertyInt( + "axisType", "Z Axis", + help=ivw.unindentMd2doc(r'''Type of the z axis, that is depth. + +* _Slice Indices_ zero-based indices of the 2D slice images +* _From Dataset__ use the `ctype3` and `cunit3` fields of the FITS dataset, e.g. frequency [Hz] +* _Velocity_ velocity derived from the frequency channels [km/s]'''), + options=options, selectedIndex=1) + self.axisScaling = ivw.properties.FloatProperty( + "axisScaling", "Z Scaling", + ivw.md2doc("Scaling factor applied to the z axis"), + 1.0, increment=0.001, min=(1.0e-8, cb.Immutable), max=(10.0, cb.Ignore)) + + self.zaxisRange = ivw.properties.DoubleVec2Property( + "zaxisRange", "Z Axis Range", + ivw.md2doc("Extent of the dataset along the z axis."), + glm.dvec2(0.0, 1.0), increment=glm.dvec2(0.001), + min=(glm.dvec2(-1.70e308), cb.Ignore), + max=(glm.dvec2(1.79e308), cb.Ignore), + semantics=semantics.Text, invalidationLevel=inv.Valid) + self.zaxisRange.readOnly = True + + self.fitsParameters = ivwastrovis.createFitsCompositeProperty( + "fitsParameters", "Fits Parameters") self.basis = ivwbase.properties.BasisProperty("basis", "Basis and offset") self.information = ivwbase.properties.VolumeInformationProperty( "information", "Data information") - self.tf = ivw.properties.TransferFunctionProperty("colormap", "Mesh Colormap", - ivw.doc.Document("Colormap used for coloring the frequency slices according to depth")) # noqa e501 - self.addProperties([self.filename, self.fitsParameters, self.latLonCoords, self.tf, + self.addProperties([self.filename, self.fitsParameters, + self.latLonCoords, self.axisType, self.zaxisRange, self.axisScaling, self.basis, self.information]) @staticmethod @@ -380,7 +104,7 @@ def processorInfo(): displayName="FitsVolumeSource", category="Python", codeState=ivw.CodeState.Stable, - tags=ivw.Tags.PY, + tags=ivw.Tags("PY, FITS, Astrophysics"), help=ivw.unindentMd2doc(FitsVolumeSource.__doc__) ) @@ -391,32 +115,21 @@ def initializeResources(self): pass @staticmethod - def createVolume(data: np.array) -> ivw.data.Volume: + def createVolume(fits_data: astroviscommon.FitsData, flip_z: bool = False) -> ivw.data.Volume: # data = np.nan_to_num(data) - volumerep = ivw.data.VolumePy(np.copy(data.astype(dtype=np.float32), order='C')) + data_float = fits_data.data.astype(dtype=np.float32) + if flip_z: + data_float = np.flip(data_float, axis=0) + volumerep = ivw.data.VolumePy(np.copy(data_float, order='C')) volume = ivw.data.Volume(volumerep) volume.swizzlemask = ivw.data.SwizzleMask.defaultData(1) - # TODO: extract extent from hdul[0].header - volume.basis = ivw.glm.mat3(1) - volume.offset = ivw.glm.vec3(-0.5, -0.5, -0.5) - - datarange = ivw.glm.dvec2(np.nanmin(data), np.nanmax(data)) + datarange = ivw.glm.dvec2(np.nanmin(fits_data.data), np.nanmax(fits_data.data)) volume.dataMap.dataRange = datarange volume.dataMap.valueRange = datarange return volume - @staticmethod - def frequencyToVelocityFunc(params: FitsParams) -> Callable[[float], int]: - def func(channel: int) -> float: - freq: int = (params.rest_frequency - - ((channel + 1 - params.crpix[2]) * params.cdelt[2] + params.crval[2])) - return vc * freq / params.rest_frequency / 1000 - - vc = 299_792_458 # speed of light - return func - def process(self): if not self.filename.value: return @@ -425,71 +138,114 @@ def process(self): if not filename.exists(): return - with fits.open(filename) as hdul: - dim = hdul[0].data.shape - - fits_parameters: FitsParams = extractFitsParams(hdul[0].header) - - func: Callable[[float], int] = FitsVolumeSource.frequencyToVelocityFunc(fits_parameters) - velocity: list = [func(channel) for channel in range(dim[0])] - print(f'Velocity matching frequency: {velocity}') - - depth_mesh, min_max = \ - createDepthMesh(fits_parameters, func, self.tf.value, - LatLong(self.latLonCoords.value), - assumed_distance=self.fitsParameters.assumedDistance.value, - velocity_inf=self.fitsParameters.velocityInf.value, - velocity_star=self.fitsParameters.velocityStar.value) - - volume_intensity = FitsVolumeSource.createVolume(hdul[0].data) - volume_intensity.dataMap.valueAxis.name = 'Intensity' - # hdul[0].header['bunit'] is 'Jy/beam ', which is not a valid unit in Inviwo - volume_intensity.dataMap.valueAxis.unit = ivw.data.Unit( - fits_parameters.bunit.replace('/beam', '')) - - extent, offset = getLatLongBasis(fits_parameters, LatLong(self.latLonCoords.value)) - - m = ivw.glm.mat4(1.0) - m[0][0] = extent[0] - m[1][1] = extent[1] - m[2][2] = np.abs(velocity[-1] - velocity[0]) - m[3] = ivw.glm.vec4(offset[0], offset[1], velocity[0], 1.0) - - axes = [("x", "\""), ("y", "\""), ("Velocity", "km/s")] - for i, (label, unit) in enumerate(axes): - volume_intensity.axes[i].name = label - volume_intensity.axes[i].unit = ivw.data.Unit(unit) - - w = ivw.glm.mat4(1.0) - w[2][2] = 2.0 / m[2][2] - # w[3] = [-0.5, -0.5, 0.0, 1.0] - volume_intensity.modelMatrix = m - volume_intensity.worldMatrix = w - - self.fitsParameters.crpix3.value = fits_parameters.crpix[2] - self.fitsParameters.crval3.value = fits_parameters.crval[2] - self.fitsParameters.cdelt3.value = fits_parameters.cdelt[2] - self.fitsParameters.restFrequency.value = fits_parameters.rest_frequency - self.fitsParameters.dataRange.value = volume_intensity.dataMap.dataRange - - extent, offset = getLatLongBasis(fits_parameters, LatLong.RelativeArcsec) - extent *= self.fitsParameters.assumedDistance.value - offset *= self.fitsParameters.assumedDistance.value - self.fitsParameters.dataExtent.value = glm.vec3(extent[0], extent[1], - min_max[1] - min_max[0]) - self.fitsParameters.dataOffset.value = glm.vec3(offset[0], offset[1], min_max[0]) - - from inviwopy.properties import OverwriteState as ows - - self.basis.updateForNewEntity(volume_intensity, self.deserializing) - self.information.updateForNewVolume( - volume_intensity, ows.Yes if self.deserializing else ows.No) - self.basis.updateEntity(volume_intensity) - self.information.updateVolume(volume_intensity) - self.deserializing = False - - self.outport.setData(volume_intensity) - self.depth.setData(depth_mesh) + astroviscommon.readFITSData.cache_clear() + + t = time.perf_counter_ns() + fits_data: astroviscommon.FitsData = astroviscommon.readFITSData(filename) + dim = fits_data.data.shape + elapsed = time.perf_counter_ns() - t + print(f'{elapsed/1000.0} us') + + print(astroviscommon.readFITSData.cache_info()) + + vel_func: Callable[[float], int] = astroviscommon.frequencyToVelocityFunc(fits_data.params) + velocity: list = [vel_func(channel) for channel in range(dim[0])] + print(f'Velocity matching frequency: {velocity}') + + extent_xy, offset_xy, unit = astroviscommon.getLatLongBasis( + fits_data.params, LatLong(self.latLonCoords.value)) + + axes = [("x", unit), ("y", unit), ("z", "")] + + extent_z: float = 0.0 + offset_z: float = 0.0 + # TODO: cannot swap basis function here as this will affect the slice order! + match AxisType(self.axisType.value): + case AxisType.SliceIndex: + extent_z = float(dim[0] - 1) + offset_z = 0 + axes[2] = ("Channel", "") + case AxisType.Dataset: + freq_func = astroviscommon.frequencyFunc(fits_data.params) + freq_range = (freq_func(0), freq_func(dim[0] - 1)) + extent_z = freq_range[1] - freq_range[0] + offset_z = freq_range[0] + axes[2] = (fits_data.params.ctype[2], fits_data.params.cunit[2]) + case AxisType.Velocity: + vel_range = (vel_func(0), vel_func(dim[0] - 1)) + extent_z = vel_range[1] - vel_range[0] + offset_z = vel_range[0] + axes[2] = ("Velocity", "km/s") + case _: + raise ValueError(f"Unexpected axis type {self.axisType.value}") + + # ensure positive axis extent in z direction + flip_z: bool = False + if extent_z < 0.0: + offset_z += extent_z + extent_z = -extent_z + flip_z = True + + volume_intensity = FitsVolumeSource.createVolume(fits_data, flip_z) + volume_intensity.dataMap.valueAxis.name = fits_data.params.btype + # hdul[0].header['bunit'] is 'Jy/beam ', which is not a valid unit in Inviwo + volume_intensity.dataMap.valueAxis.unit = ivw.data.Unit( + fits_data.params.bunit.replace('/beam', '')) + + for i, (label, unit) in enumerate(axes): + volume_intensity.axes[i].name = label + volume_intensity.axes[i].unit = ivw.data.Unit(unit) + + m = ivw.glm.dmat4(1.0) + m[0][0] = extent_xy[0] + m[1][1] = extent_xy[1] + m[2][2] = extent_z + m[3] = ivw.glm.dvec4(offset_xy[0], offset_xy[1], offset_z, 1.0) + + # re-normalize model transformations using the world matrix + w = ivw.glm.dmat4(1.0) + w[0][0] = 1.0 / extent_xy[0] + w[1][1] = 1.0 / extent_xy[1] + w[2][2] = self.axisScaling.value / extent_z + w[3] = [-0.5 - m[3][0] * w[0][0], -0.5 - m[3][1] * w[1][1], -0.5 - m[3][2] * w[2][2], 1.0] + volume_intensity.modelMatrix = ivw.glm.mat4(m) + volume_intensity.worldMatrix = ivw.glm.mat4(w) + + # print(w * m, extent_z * (self.axisScaling.value / extent_z)) + + self.zaxisRange.value = glm.dvec2(offset_z, offset_z + extent_z) + + self.fitsParameters.objectName.value = fits_data.params.objectName + self.fitsParameters.crpix3.value = fits_data.params.crpix[2] + self.fitsParameters.crval3.value = fits_data.params.crval[2] + self.fitsParameters.cdelt3.value = fits_data.params.cdelt[2] + self.fitsParameters.restFrequency.value = fits_data.params.rest_frequency + self.fitsParameters.valueName.value = volume_intensity.dataMap.valueAxis.name + self.fitsParameters.valueUnit.value = str(volume_intensity.dataMap.valueAxis.unit) + self.fitsParameters.dataRange.dataRange = volume_intensity.dataMap.dataRange + # update data range of volume in case of a custom range + volume_intensity.dataMap.dataRange = self.fitsParameters.dataRange.dataRange + volume_intensity.dataMap.valueRange = self.fitsParameters.dataRange.dataRange + + # determine lat/long extent in au + extent_xy, offset_xy, _ = astroviscommon.getLatLongBasis( + fits_data.params, LatLong.RelativeArcsec) + extent_xy *= self.fitsParameters.assumedDistance.value + offset_xy *= self.fitsParameters.assumedDistance.value + self.fitsParameters.dataExtent.value = glm.vec3(extent_xy[0], extent_xy[1], 0.0) + self.fitsParameters.dataOffset.value = glm.vec3(offset_xy[0], offset_xy[1], 0.0) + + from inviwopy.properties import OverwriteState as ows + + self.basis.updateForNewEntity(volume_intensity, self.deserializing) + self.information.updateForNewVolume( + volume_intensity, ows.Yes if self.deserializing else ows.No) + self.basis.updateEntity(volume_intensity) + self.information.updateVolume(volume_intensity) + self.deserializing = False + + self.outports.intensity.setData(volume_intensity) + self.outports.fitsdata.setData(fits_data) def deserialize(self, deserializer: ivw.Deserializer): ivw.Processor.deserialize(self, deserializer) From ef04eca335f5d47af45e001ed6634aa7a6f4872b Mon Sep 17 00:00:00 2001 From: Martin Falk Date: Wed, 25 Mar 2026 18:32:58 +0100 Subject: [PATCH 03/16] AstroPhysics: raycasting, depth mesh trafo --- .../fitsnonlineardepthcomponent.h | 6 +- modules/astrophysics/python/ivwastrovis.py | 49 ++++-- .../python/processors/FitsDepthCorrection.py | 38 ++++- .../python/processors/FitsVolumeSource.py | 6 - .../fitsnonlineardepthcomponent.cpp | 149 +++++++++--------- 5 files changed, 153 insertions(+), 95 deletions(-) diff --git a/modules/astrophysics/include/infravis/astrophysics/shadercomponents/fitsnonlineardepthcomponent.h b/modules/astrophysics/include/infravis/astrophysics/shadercomponents/fitsnonlineardepthcomponent.h index 7628dceab..76f742763 100644 --- a/modules/astrophysics/include/infravis/astrophysics/shadercomponents/fitsnonlineardepthcomponent.h +++ b/modules/astrophysics/include/infravis/astrophysics/shadercomponents/fitsnonlineardepthcomponent.h @@ -31,7 +31,7 @@ #include #include -#include +#include #include #include #include @@ -61,6 +61,7 @@ class IVW_MODULE_ASTROPHYSICS_API FitsNonlinearDepthComponent : public ShaderCom private: CompositeProperty fitsParameters_; + StringProperty objectName_; FloatProperty assumedDistance_; FloatProperty velocityInf_; FloatProperty velocityStar_; @@ -70,6 +71,9 @@ class IVW_MODULE_ASTROPHYSICS_API FitsNonlinearDepthComponent : public ShaderCom FloatProperty restFrequency_; FloatVec3Property dataExtent_; FloatVec3Property dataOffset_; + IntVec2Property validChannelRange_; + FloatVec3Property invalidColor_; + FloatProperty invalidAlpha_; }; } // namespace inviwo diff --git a/modules/astrophysics/python/ivwastrovis.py b/modules/astrophysics/python/ivwastrovis.py index 90ad15cff..44f434c33 100644 --- a/modules/astrophysics/python/ivwastrovis.py +++ b/modules/astrophysics/python/ivwastrovis.py @@ -16,7 +16,29 @@ def createDepthMesh(params: FitsParams, assumed_distance: float = 123.0, velocity_inf: float = 14.5, velocity_star: float = -26.5 - ) -> (ivw.data.Mesh, np.ndarray): + ) -> (ivw.data.Mesh, np.ndarray, (int, int)): + """ + Create a Mesh containing the individual channels with the distorted z axis based on a + depth estimation of the channel velocity + + Parameters + ---------- + params : FitsParams + Parameters extracted from a fits.header.Header + tf : ivw.data.TransferFunction + Transferfunction to color the mesh according to depth + assumed_distance : float, optional + in parsecs. The default is 123.0. + velocity_inf: float, optional + in km/s. The default is 14.5. + velocity_star: float, optional + in km/s. The default is -26.5. + + Returns + ------- + A tuple containing the Mesh, min/max depth values (np.ndarray), and + first and last valid channel index (int, int) where |velocity - velocity_star| <= velocity_inf + """ velocity_func: Callable[[float], int] = astroviscommon.frequencyToVelocityFunc(params) @@ -24,9 +46,15 @@ def createDepthMesh(params: FitsParams, positions: list[np.ndarray] = [] + v = np.array([velocity_func(i) for i in range(0, dims[2])]) + valid_channels = np.argwhere(np.abs(v - velocity_star) <= velocity_inf).squeeze() + if len(valid_channels) == 0: + raise ValueError( + 'No valid velocity channels found based on |velocity - velocity_star| <= velocity_inf') + mesh_count: int = 0 mesh = ivw.data.Mesh(dt=ivw.data.DrawType.Triangles, ct=ivw.data.ConnectivityType.Unconnected) - for i in range(0, dims[2]): + for i in valid_channels: with suppress(ValueError): with np.nditer(astroviscommon.estimateDepth(params, velocity_func(i)), flags=['multi_index']) as it: @@ -46,17 +74,16 @@ def createDepthMesh(params: FitsParams, mesh_count += 1 if len(positions) == 0: - return (mesh, np.array([0, 0])) + return (mesh, np.array([0, 0]), (valid_channels[0], valid_channels[-1])) positions_np = np.array(positions) min_max = [np.min(positions_np[:, 2]), np.max(positions_np[:, 2])] - ivw.logInfo(f'min/max: {min_max}') + # ivw.logInfo(f'min/max: {min_max}') depth_scaling = np.max(np.abs(min_max)) - # scaling = np.array([1.0 / dims[0], 1.0 / dims[1], 1.0 / depth_scaling]) colors = [tf.sample(value / depth_scaling * 0.5 + 0.5) for value in positions_np[:, 2]] - extent, offset = astroviscommon.getLatLongBasis(params, lat_long) + extent, offset, _ = astroviscommon.getLatLongBasis(params, lat_long) ivw.logInfo(f'{extent=}\n{offset=}') positions_np = (positions_np * np.array([extent[0] / (dims[0] - 1), extent[1] / (dims[1] - 1), 1.0]) @@ -75,13 +102,15 @@ def createDepthMesh(params: FitsParams, m = ivw.glm.mat4(1.0) w = ivw.glm.mat4(1.0) - w[2][2] = 1.0 / depth_scaling * 20.0 - # w[3] = [-0.5, -0.5, 0.0, 1.0] + # w[0][0] = 1.0 / extent[0] + # w[1][1] = 1.0 / extent[1] + # w[2][2] = 2.0 / (min_max[1] - min_max[0]) + # w[3] = [-0.5, -0.5, -0.5, 1.0] mesh.modelMatrix = m mesh.worldMatrix = w - return (mesh, np.array(min_max)) + return (mesh, np.array(min_max), (valid_channels[0], valid_channels[-1])) def createFitsCompositeProperty(identifier: str, @@ -129,7 +158,7 @@ def createFitsCompositeProperty(identifier: str, min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), semantics=semantics.Text, invalidationLevel=inv.Valid) - restFrequency = ivw.properties.FloatProperty("restFrequency", "Rest Frequency [Hz", + restFrequency = ivw.properties.FloatProperty("restFrequency", "Rest Frequency [Hz]", ivw.md2doc("restfreq [Hz]"), 0.0, increment=0.001, min=(-100.0, cb.Ignore), max=(100.0, cb.Ignore), diff --git a/modules/astrophysics/python/processors/FitsDepthCorrection.py b/modules/astrophysics/python/processors/FitsDepthCorrection.py index 6f518a5f3..f8be2aa41 100644 --- a/modules/astrophysics/python/processors/FitsDepthCorrection.py +++ b/modules/astrophysics/python/processors/FitsDepthCorrection.py @@ -20,16 +20,29 @@ class FitsDepthCorrection(ivw.Processor): def __init__(self, id, name): ivw.Processor.__init__(self, id, name) + cb = ivw.properties.ConstraintBehavior + semantics = ivw.properties.PropertySemantics + self.addInport(ivw.PythonInport("fitsdata")) self.addOutport(ivw.data.MeshOutport("depthmesh")) self.fitsParameters = ivwastrovis.createFitsCompositeProperty( "fitsParameters", "Fits Parameters") + self.validChannelRange = ivw.properties.IntVec2Property( + "validChannelRange", "Valid Channel Range", + ivw.doc.Document("Range of valid channels based on the line-of-sight velocity"), + glm.ivec2(0, 0), + min=(glm.ivec2(0), cb.Ignore), max=(glm.ivec2(100), cb.Ignore), + semantics=semantics.Text) + self.axisScaling = ivw.properties.FloatProperty( + "axisScaling", "Z Scaling", + ivw.md2doc("Scaling factor applied to the z axis"), + 1.0, increment=0.001, min=(1.0e-8, cb.Immutable), max=(10.0, cb.Ignore)) self.tf = ivw.properties.TransferFunctionProperty( "colormap", "Mesh Colormap", ivw.doc.Document("Colormap used for coloring the frequency slices according to depth")) - self.addProperties([self.fitsParameters, self.tf]) + self.addProperties([self.fitsParameters, self.validChannelRange, self.axisScaling, self.tf]) @staticmethod def processorInfo(): @@ -55,7 +68,7 @@ def process(self): f"expected '{astroviscommon.FitsData}'.") t = time.perf_counter_ns() - depth_mesh, min_max = ivwastrovis.createDepthMesh( + depth_mesh, min_max, valid_channel_range = ivwastrovis.createDepthMesh( fits_data.params, self.tf.value, LatLong.RelativeArcsec, assumed_distance=self.fitsParameters.assumedDistance.value, @@ -64,11 +77,30 @@ def process(self): elapsed = time.perf_counter_ns() - t print(f'{elapsed/1000_000.0} ms') - extent, offset = astroviscommon.getLatLongBasis(fits_data.params, LatLong.RelativeArcsec) + extent, offset, unit = astroviscommon.getLatLongBasis( + fits_data.params, LatLong.RelativeArcsec) + axes = [("x", unit), ("y", unit), ("z", "au")] + for i, (label, unit) in enumerate(axes): + depth_mesh.axes[i].name = label + depth_mesh.axes[i].unit = ivw.data.Unit(unit) + + # m = ivw.glm.mat4(1.0) + w = ivw.glm.mat4(1.0) + w[0][0] = 1.0 / extent[0] + w[1][1] = 1.0 / extent[1] + w[2][2] = self.axisScaling.value / (min_max[1] - min_max[0]) + w[3] = [-0.5 + offset[0] * w[0][0], -0.5 + offset[1] * w[1][1], + -0.5 + (-min_max[0]) * w[2][2], 1.0] + + # depth_mesh.modelMatrix = m + depth_mesh.worldMatrix = w + + extent, offset, _ = astroviscommon.getLatLongBasis(fits_data.params, LatLong.RelativeArcsec) extent *= self.fitsParameters.assumedDistance.value offset *= self.fitsParameters.assumedDistance.value self.fitsParameters.dataExtent.value = glm.vec3( extent[0], extent[1], min_max[1] - min_max[0]) self.fitsParameters.dataOffset.value = glm.vec3(offset[0], offset[1], min_max[0]) + self.validChannelRange.value = glm.ivec2(valid_channel_range) self.outports.depthmesh.setData(depth_mesh) diff --git a/modules/astrophysics/python/processors/FitsVolumeSource.py b/modules/astrophysics/python/processors/FitsVolumeSource.py index bf1cc46da..87484b440 100644 --- a/modules/astrophysics/python/processors/FitsVolumeSource.py +++ b/modules/astrophysics/python/processors/FitsVolumeSource.py @@ -30,12 +30,6 @@ def __init__(self, id, name): ivw.Processor.__init__(self, id, name) self.deserializing = False - global astroviscommon - global LatLong - from importlib import reload - astroviscommon = reload(astroviscommon) - from astroviscommon import LatLong - self.addOutport(ivw.data.VolumeOutport("intensity")) self.addOutport(ivw.PythonOutport("fitsdata")) diff --git a/modules/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp b/modules/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp index edc86c0b2..6aebaf923 100644 --- a/modules/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp +++ b/modules/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp @@ -44,6 +44,9 @@ FitsNonlinearDepthComponent::FitsNonlinearDepthComponent(std::string_view name, , volumePort_{name, std::move(help)} , fitsParameters_{"fitsParameters", "FitsParameters", "Various parameters required for coordinate transformations and other calculations."_help} + , objectName_{"objectName", "Object", + "Contents of the 'objects' field in the FITS Header."_help, "", + InvalidationLevel::Valid} , assumedDistance_{"assumedDistance", "Assumed Distance [pc]", util::ordinalSymmetricVector(0.0f) .setInc(0.001f) @@ -88,10 +91,31 @@ FitsNonlinearDepthComponent::FitsNonlinearDepthComponent(std::string_view name, util::ordinalSymmetricVector(vec3{0.0f}) .setInc(vec3{0.001f}) .set(PropertySemantics::Text) - .set("Coordinate of the lower left corner of the dataset in au"_help)} { - - fitsParameters_.addProperties(assumedDistance_, velocityInf_, velocityStar_, crpix3_, crval3_, - cdelt3_, restFrequency_, dataExtent_, dataOffset_); + .set("Coordinate of the lower left corner of the dataset in au"_help)} + , validChannelRange_{"validChannelRange", "Valid Channel Range", + OrdinalPropertyState{ + .value = ivec2{0}, + .min = ivec2{0}, + .minConstraint = ConstraintBehavior::Ignore, + .max = ivec2{100}, + .maxConstraint = ConstraintBehavior::Ignore, + .semantics = PropertySemantics::Text, + .help = + "Range of valid channels based on the line-of-sight velocity"_help, + }} + + , invalidColor_{"invalidColor", "Invalid Value Color", + util::ordinalColor(vec3{0.2f, 0.2f, 0.2f})} + , invalidAlpha_{"invalidAlpha", + "Invalid Value Alpha", + 0.005f, + {0.0f, ConstraintBehavior::Immutable}, + {1.0f, ConstraintBehavior::Immutable}, + 0.001f} { + + fitsParameters_.addProperties(objectName_, assumedDistance_, velocityInf_, velocityStar_, + crpix3_, crval3_, cdelt3_, restFrequency_, dataExtent_, + dataOffset_, validChannelRange_); fitsParameters_.setCollapsed(true); } @@ -104,8 +128,8 @@ void FitsNonlinearDepthComponent::process(Shader& shader, TextureUnitContainer& const auto data = volumePort_.getData(); const vec3 delta = dataExtent_.get() / vec3{data->getDimensions()}; - mat4 m = glm::scale(delta * vec3{data->getDimensions() + size3_t{1u}}); - m[3] = vec4{dataOffset_.get() - delta * 0.5f, 1.0f}; + mat4 m = glm::scale(dataExtent_.get()); + m[3] = vec4{dataOffset_.get(), 1.0f}; StrBuffer buf{"{}FitsParameters", getName()}; std::string prefix{buf}; @@ -117,6 +141,8 @@ void FitsNonlinearDepthComponent::process(Shader& shader, TextureUnitContainer& shader.setUniform(buf.replace("{}.crval3", prefix), crval3_); shader.setUniform(buf.replace("{}.cdelt3", prefix), cdelt3_); shader.setUniform(buf.replace("{}.restFrequency", prefix), restFrequency_); + shader.setUniform(buf.replace("{}.validChannelRange", prefix), validChannelRange_); + shader.setUniform("invalidColor", vec4{invalidColor_.get(), invalidAlpha_.get()}); } std::vector> FitsNonlinearDepthComponent::getInports() { @@ -127,7 +153,9 @@ void FitsNonlinearDepthComponent::initializeResources(Shader&) { // shader definitions... } -std::vector FitsNonlinearDepthComponent::getProperties() { return {&fitsParameters_}; } +std::vector FitsNonlinearDepthComponent::getProperties() { + return {&fitsParameters_, &invalidColor_, &invalidAlpha_}; +} namespace { @@ -143,6 +171,8 @@ struct FitsParameters { float crval3; float cdelt3; float restFrequency; // [Hz] + + ivec2 validChannelRange; // range of valid channels based on line-of-sight velocity }; const float invalidValue = -1.e20; @@ -153,6 +183,7 @@ uniform VolumeParameters {0}Parameters; uniform sampler3D {0}; uniform FitsParameters {0}FitsParameters; +uniform vec4 invalidColor = vec4(0.2, 0.2, 0.2, 0.005); )"); constexpr std::string_view utilityFuncs = util::trim(R"( @@ -163,6 +194,14 @@ float velocityToChannel(in FitsParameters params, in float velocity) {{ (1.0 - velocity * 1000 / vc) - params.crval3) / params.cdelt3 + params.crpix3 - 1)); }} +// convert a frequency channel index to velocity [km/s] +float channelToVelocity(in FitsParameters params, in int channel) {{ + float vc = 299792458; // speed of light [m/s] + float freq = (params.restFrequency - + ((channel + 1 - params.crpix3) * params.cdelt3 + params.crval3)); + return vc * freq / params.restFrequency / 1000.0; +}} + // Convert a relative spatial coordinate [au] to velocity [km/s] float spatialCoordinateToVelocity(in FitsParameters params, in vec3 coord) {{ // convert xy coord from arcseconds to au @@ -170,25 +209,26 @@ float spatialCoordinateToVelocity(in FitsParameters params, in vec3 coord) {{ return params.velocityInf * coord.z / length(coord) + params.velocityStar; }} -// TODO: function from normalized texture coords to normalized, non-linear frequency domain -// using velocityToChannel(spatialCoordinateToVelocity(...)) -// requires extent in au - // Sample the volume texture by transforming normalized texture coordinates @p texCoord into // the non-linear frequency domain using velocityToChannel(spatialCoordinateToVelocity(...)). // // @param texCoord normalized texture coordinate [0, 1] // @return volume sample at re-normalized texture coordinate, if v_los > v_star. // Otherwise vec4(invalidValue); -vec4 getNormalizedVoxel(in FitsParameters params, in vec3 texCoord) {{ +vec4 getNormalizedVoxel(in FitsParameters params, in vec3 texCoord, in vec2 minMaxVelocity) {{ vec3 coord = (params.textureToModel * vec4(texCoord, 1.0)).xyz; float velocity = spatialCoordinateToVelocity(params, coord); + float u = velocityToChannel(params, velocity) / textureSize(volume, 0).z; + float v_normalized = (velocity - minMaxVelocity.x) / (minMaxVelocity.y - minMaxVelocity.x); + float v_los = velocity - params.velocityStar; - if (abs(v_los) > params.velocityInf) {{ + if (abs(v_los) > params.velocityInf || u < 0.0 || u > 1.0 + || v_normalized < 0.0 || v_normalized > 1.0) {{ return vec4(invalidValue); }} - texCoord = vec3(texCoord.xy, velocityToChannel(params, velocity) / textureSize(volume, 0).z); + + texCoord = vec3(texCoord.st, u); return getNormalizedVoxel({0}, {0}Parameters, texCoord); }} )"); @@ -197,61 +237,29 @@ vec4 getNormalizedVoxel(in FitsParameters params, in vec3 texCoord) {{ // mainly for the isosurface rendering. Setting it to the same voxel value prevents isosurfaces // being rendered at the volume boundaries. constexpr std::string_view voxelFirst = util::trim(R"( -vec4 {0}Voxel = getNormalizedVoxel({0}FitsParameters, samplePosition); +vec2 minMaxVelocity = vec2(channelToVelocity({0}FitsParameters, {0}FitsParameters.validChannelRange.x), + channelToVelocity({0}FitsParameters, {0}FitsParameters.validChannelRange.y)); +if (minMaxVelocity.x > minMaxVelocity.y) {{ + minMaxVelocity.xy = minMaxVelocity.yx; +}} + +vec4 {0}Voxel = getNormalizedVoxel({0}FitsParameters, samplePosition, minMaxVelocity); vec4 {0}VoxelPrev = {0}Voxel; )"); constexpr std::string_view voxel = util::trim(R"( {0}VoxelPrev = {0}Voxel; -{0}Voxel = getNormalizedVoxel({0}FitsParameters, samplePosition); -)"); - -constexpr std::string_view gradientFirst = util::trim(R"( -vec3 {0}GradientPrev = vec3(0); -vec3 {0}Gradient = vec3(0); -#if defined(GRADIENTS_ENABLED) -{0}Gradient = useSurfaceNormals ? -texture(surfaceNormal, texCoords).xyz : - normalize(COMPUTE_GRADIENT_FOR_CHANNEL({0}Voxel, {0}, {0}Parameters, - samplePosition, channel)); -if (!useSurfaceNormals) {{ - {0}Gradient *= sign({0}Voxel[channel] / {0}Parameters.texToNormalized.scale + {0}Parameters.texToNormalized.offset); -}} -#endif -)"); - -constexpr std::string_view gradient = util::trim(R"( -#if defined(GRADIENTS_ENABLED) -{0}GradientPrev = {0}Gradient; -{0}Gradient = normalize(COMPUTE_GRADIENT_FOR_CHANNEL({0}Voxel, {0}, {0}Parameters, - samplePosition, channel)); -{0}Gradient *= sign({0}Voxel[channel] / {0}Parameters.texToNormalized.scale + {0}Parameters.texToNormalized.offset); -#endif +{0}Voxel = getNormalizedVoxel({0}FitsParameters, samplePosition, minMaxVelocity); )"); -constexpr std::string_view allGradientsFirst = util::trim(R"( -mat4x3 {0}AllGradientsPrev = mat4x3(0); -mat4x3 {0}AllGradients = mat4x3(0); -#if defined(GRADIENTS_ENABLED) -vec3 surfaceNormal = useSurfaceNormals ? -texture(surfaceNormal, texCoords).xyz : vec3(0); -{0}AllGradients = useSurfaceNormals ? - mat4x3(surfaceNormal, surfaceNormal, surfaceNormal, surfaceNormal) : - COMPUTE_ALL_GRADIENTS({0}Voxel, {0}, {0}Parameters, samplePosition); -{0}AllGradients[0] = normalize({0}AllGradients[0]); -{0}AllGradients[1] = normalize({0}AllGradients[1]); -{0}AllGradients[2] = normalize({0}AllGradients[2]); -{0}AllGradients[3] = normalize({0}AllGradients[3]); -#endif +constexpr std::string_view unsetColor = util::trim(R"( +color = vec4(0.0); )"); -constexpr std::string_view allGradients = util::trim(R"( -#if defined(GRADIENTS_ENABLED) -{0}AllGradientsPrev = {0}AllGradients; -{0}AllGradients = COMPUTE_ALL_GRADIENTS({0}Voxel, {0}, {0}Parameters, samplePosition); -{0}AllGradients[0] = normalize({0}AllGradients[0]); -{0}AllGradients[1] = normalize({0}AllGradients[1]); -{0}AllGradients[2] = normalize({0}AllGradients[2]); -{0}AllGradients[3] = normalize({0}AllGradients[3]); -#endif +constexpr std::string_view invalidValue = util::trim(R"( +if (volumeVoxel[channel] <= invalidValue) {{ + color = invalidColor; +}} )"); } // namespace @@ -270,26 +278,17 @@ auto FitsNonlinearDepthComponent::getSegments() -> std::vector { {.snippet = fmt::format(voxelFirst, getName()), .placeholder = placeholder::first, .priority = 400}, + {.snippet = fmt::format(unsetColor, getName()), + .placeholder = placeholder::first, + .priority = 701}, {.snippet = fmt::format(voxel, getName()), .placeholder = placeholder::loop, .priority = 400}, + {.snippet = fmt::format(invalidValue, getName()), + .placeholder = placeholder::loop, + .priority = 701}, }; - // if (gradients != Gradients::None) { - // segments.emplace_back(std::string{R"(#include "utils/gradients.glsl")"}, - // placeholder::include, 400); - // } - // if (gradients == Gradients::Single) { - // segments.emplace_back(fmt::format(gradientFirst, getName()), placeholder::first, 410); - // segments.emplace_back(fmt::format(gradient, getName()), placeholder::loop, 410); - // } - - // if (gradients == Gradients::All) { - // segments.emplace_back(fmt::format(allGradientsFirst, getName()), placeholder::first, - // 410); segments.emplace_back(fmt::format(allGradients, getName()), placeholder::loop, - // 410); - // } - return segments; } From 15110333348ded1bd2e50aac4f195cacd07f39df Mon Sep 17 00:00:00 2001 From: Martin Falk Date: Mon, 30 Mar 2026 16:15:54 +0200 Subject: [PATCH 04/16] AstroPhysics: depth mesh trafo --- modules/astrophysics/python/ivwastrovis.py | 2 +- .../python/processors/FitsDepthCorrection.py | 62 ++++++++++++++----- .../python/processors/FitsVolumeSource.py | 13 +--- 3 files changed, 50 insertions(+), 27 deletions(-) diff --git a/modules/astrophysics/python/ivwastrovis.py b/modules/astrophysics/python/ivwastrovis.py index 44f434c33..9cf7971c8 100644 --- a/modules/astrophysics/python/ivwastrovis.py +++ b/modules/astrophysics/python/ivwastrovis.py @@ -87,7 +87,7 @@ def createDepthMesh(params: FitsParams, ivw.logInfo(f'{extent=}\n{offset=}') positions_np = (positions_np * np.array([extent[0] / (dims[0] - 1), extent[1] / (dims[1] - 1), 1.0]) - - np.array([offset[0], offset[1], 0.0])) + + np.array([offset[0], offset[1], 0.0])) mesh.addBuffer(ivw.data.BufferType.PositionAttrib, ivw.data.Buffer(np.array(positions_np, dtype=np.float32))) diff --git a/modules/astrophysics/python/processors/FitsDepthCorrection.py b/modules/astrophysics/python/processors/FitsDepthCorrection.py index f8be2aa41..16ab12442 100644 --- a/modules/astrophysics/python/processors/FitsDepthCorrection.py +++ b/modules/astrophysics/python/processors/FitsDepthCorrection.py @@ -21,6 +21,7 @@ def __init__(self, id, name): ivw.Processor.__init__(self, id, name) cb = ivw.properties.ConstraintBehavior + IntOption = ivw.properties.IntOption semantics = ivw.properties.PropertySemantics self.addInport(ivw.PythonInport("fitsdata")) @@ -28,6 +29,20 @@ def __init__(self, id, name): self.fitsParameters = ivwastrovis.createFitsCompositeProperty( "fitsParameters", "Fits Parameters") + + options = [IntOption("absolute", "Absolute (degree)", LatLong.AbsoluteDeg.value), + IntOption("fractionalArcsec", "Fractional (arcsec)", + LatLong.FractionalArcsec.value), + IntOption("relativeArcsec", "Relative (arcsec)", LatLong.RelativeArcsec.value)] + self.latLonCoords = ivw.properties.OptionPropertyInt( + "latLonCoords", "Lat/Long Coordinates", + help=ivw.unindentMd2doc(r'''Defines how Latitude and Longitude coordinates are shown. + +* _Absolute_ full coordinate in decimal degrees +* _Relative (arcsec)_ coordinates in arcseconds relative to the star/center +* _Relative (arcsec)_ coordinates in arcseconds relative to the star/center'''), + options=options, selectedIndex=2) + self.validChannelRange = ivw.properties.IntVec2Property( "validChannelRange", "Valid Channel Range", ivw.doc.Document("Range of valid channels based on the line-of-sight velocity"), @@ -42,7 +57,8 @@ def __init__(self, id, name): "colormap", "Mesh Colormap", ivw.doc.Document("Colormap used for coloring the frequency slices according to depth")) - self.addProperties([self.fitsParameters, self.validChannelRange, self.axisScaling, self.tf]) + self.addProperties([self.fitsParameters, self.validChannelRange, + self.latLonCoords, self.axisScaling, self.tf]) @staticmethod def processorInfo(): @@ -68,39 +84,53 @@ def process(self): f"expected '{astroviscommon.FitsData}'.") t = time.perf_counter_ns() + latlong_calc: LatLong = LatLong.RelativeArcsec depth_mesh, min_max, valid_channel_range = ivwastrovis.createDepthMesh( fits_data.params, self.tf.value, - LatLong.RelativeArcsec, + latlong_calc, assumed_distance=self.fitsParameters.assumedDistance.value, velocity_inf=self.fitsParameters.velocityInf.value, velocity_star=self.fitsParameters.velocityStar.value) elapsed = time.perf_counter_ns() - t print(f'{elapsed/1000_000.0} ms') + # extent and offset in same coord system as the depth mesh + mesh_extent, mesh_offset, _ = astroviscommon.getLatLongBasis(fits_data.params, latlong_calc) + extent_au = mesh_extent * self.fitsParameters.assumedDistance.value + offset_au = mesh_offset * self.fitsParameters.assumedDistance.value + self.fitsParameters.dataExtent.value = glm.vec3( + extent_au[0], extent_au[1], min_max[1] - min_max[0]) + self.fitsParameters.dataOffset.value = glm.vec3( + offset_au[0], offset_au[1], min_max[0]) + self.validChannelRange.value = glm.ivec2(valid_channel_range) + + # extent and offset for positioning and axes extent, offset, unit = astroviscommon.getLatLongBasis( - fits_data.params, LatLong.RelativeArcsec) + fits_data.params, LatLong(self.latLonCoords.value)) axes = [("x", unit), ("y", unit), ("z", "au")] for i, (label, unit) in enumerate(axes): depth_mesh.axes[i].name = label depth_mesh.axes[i].unit = ivw.data.Unit(unit) - # m = ivw.glm.mat4(1.0) + ivw.logInfo(f'd: {extent=}\n{offset=}') + + m = ivw.glm.mat4(1.0) + m[0][0] = extent[0] / mesh_extent[0] + m[1][1] = extent[1] / mesh_extent[1] + if LatLong(self.latLonCoords.value) == LatLong.AbsoluteDeg: + # TODO: why without mesh_offset? + m[3] = [offset[0], offset[1], 1.0, 1.0] + else: + m[3] = [offset[0] - mesh_offset[0], offset[1] - mesh_offset[1], 1.0, 1.0] + w = ivw.glm.mat4(1.0) w[0][0] = 1.0 / extent[0] w[1][1] = 1.0 / extent[1] w[2][2] = self.axisScaling.value / (min_max[1] - min_max[0]) - w[3] = [-0.5 + offset[0] * w[0][0], -0.5 + offset[1] * w[1][1], - -0.5 + (-min_max[0]) * w[2][2], 1.0] - - # depth_mesh.modelMatrix = m + w[3] = [- m[3][0] * w[0][0], + - m[3][1] * w[1][1], + -0.5 - min_max[0] * w[2][2], 1.0] + depth_mesh.modelMatrix = m depth_mesh.worldMatrix = w - extent, offset, _ = astroviscommon.getLatLongBasis(fits_data.params, LatLong.RelativeArcsec) - extent *= self.fitsParameters.assumedDistance.value - offset *= self.fitsParameters.assumedDistance.value - self.fitsParameters.dataExtent.value = glm.vec3( - extent[0], extent[1], min_max[1] - min_max[0]) - self.fitsParameters.dataOffset.value = glm.vec3(offset[0], offset[1], min_max[0]) - self.validChannelRange.value = glm.ivec2(valid_channel_range) - self.outports.depthmesh.setData(depth_mesh) diff --git a/modules/astrophysics/python/processors/FitsVolumeSource.py b/modules/astrophysics/python/processors/FitsVolumeSource.py index 87484b440..ccad386ea 100644 --- a/modules/astrophysics/python/processors/FitsVolumeSource.py +++ b/modules/astrophysics/python/processors/FitsVolumeSource.py @@ -112,8 +112,8 @@ def initializeResources(self): def createVolume(fits_data: astroviscommon.FitsData, flip_z: bool = False) -> ivw.data.Volume: # data = np.nan_to_num(data) data_float = fits_data.data.astype(dtype=np.float32) - if flip_z: - data_float = np.flip(data_float, axis=0) + # if flip_z: + # data_float = np.flip(data_float, axis=0) volumerep = ivw.data.VolumePy(np.copy(data_float, order='C')) volume = ivw.data.Volume(volumerep) @@ -173,14 +173,7 @@ def process(self): case _: raise ValueError(f"Unexpected axis type {self.axisType.value}") - # ensure positive axis extent in z direction - flip_z: bool = False - if extent_z < 0.0: - offset_z += extent_z - extent_z = -extent_z - flip_z = True - - volume_intensity = FitsVolumeSource.createVolume(fits_data, flip_z) + volume_intensity = FitsVolumeSource.createVolume(fits_data) volume_intensity.dataMap.valueAxis.name = fits_data.params.btype # hdul[0].header['bunit'] is 'Jy/beam ', which is not a valid unit in Inviwo volume_intensity.dataMap.valueAxis.unit = ivw.data.Unit( From 5097a9cf26c2f2ffedaa2a09d6a6ddfc97d0294e Mon Sep 17 00:00:00 2001 From: Martin Falk Date: Tue, 31 Mar 2026 09:43:03 +0200 Subject: [PATCH 05/16] AstroPhysics: expose axis ranges --- modules/astrophysics/python/astroviscommon.py | 2 +- .../python/processors/FitsDepthCorrection.py | 27 ++++++++++++-- .../python/processors/FitsVolumeSource.py | 36 +++++++++++-------- 3 files changed, 47 insertions(+), 18 deletions(-) diff --git a/modules/astrophysics/python/astroviscommon.py b/modules/astrophysics/python/astroviscommon.py index e80101196..08d9b59a9 100644 --- a/modules/astrophysics/python/astroviscommon.py +++ b/modules/astrophysics/python/astroviscommon.py @@ -46,7 +46,7 @@ class FitsData: data: np.ndarray -@lru_cache(maxsize=1) +@lru_cache(maxsize=5) def readFITSData(filepath: Path, dataset_index: int | None = None, **kwargs) -> FitsData: diff --git a/modules/astrophysics/python/processors/FitsDepthCorrection.py b/modules/astrophysics/python/processors/FitsDepthCorrection.py index 16ab12442..adea554dd 100644 --- a/modules/astrophysics/python/processors/FitsDepthCorrection.py +++ b/modules/astrophysics/python/processors/FitsDepthCorrection.py @@ -53,12 +53,30 @@ def __init__(self, id, name): "axisScaling", "Z Scaling", ivw.md2doc("Scaling factor applied to the z axis"), 1.0, increment=0.001, min=(1.0e-8, cb.Immutable), max=(10.0, cb.Ignore)) + + def axisRangeProperty(axis_name: str) -> ivw.properties.DoubleVec2Property: + inv = ivw.properties.InvalidationLevel + prop = ivw.properties.DoubleVec2Property( + f"{axis_name.lower()}axisRange", f"{axis_name.upper()} Axis Range", + ivw.md2doc(f"Extent of the dataset along the {axis_name.lower()} axis."), + glm.dvec2(0.0, 1.0), increment=glm.dvec2(0.001), + min=(glm.dvec2(-1.79e308), cb.Ignore), + max=(glm.dvec2(1.79e308), cb.Ignore), + semantics=semantics.Text, invalidationLevel=inv.Valid) + prop.readOnly = True + return prop + + self.xaxisRange = axisRangeProperty("x") + self.yaxisRange = axisRangeProperty("y") + self.zaxisRange = axisRangeProperty("z") + self.tf = ivw.properties.TransferFunctionProperty( "colormap", "Mesh Colormap", ivw.doc.Document("Colormap used for coloring the frequency slices according to depth")) - self.addProperties([self.fitsParameters, self.validChannelRange, - self.latLonCoords, self.axisScaling, self.tf]) + self.addProperties([self.fitsParameters, self.validChannelRange, self.latLonCoords, + self.xaxisRange, self.yaxisRange, self.zaxisRange, + self.axisScaling, self.tf]) @staticmethod def processorInfo(): @@ -92,7 +110,7 @@ def process(self): velocity_inf=self.fitsParameters.velocityInf.value, velocity_star=self.fitsParameters.velocityStar.value) elapsed = time.perf_counter_ns() - t - print(f'{elapsed/1000_000.0} ms') + print(f'estimating depth: {elapsed/1_000_000.0} ms') # extent and offset in same coord system as the depth mesh mesh_extent, mesh_offset, _ = astroviscommon.getLatLongBasis(fits_data.params, latlong_calc) @@ -113,6 +131,9 @@ def process(self): depth_mesh.axes[i].unit = ivw.data.Unit(unit) ivw.logInfo(f'd: {extent=}\n{offset=}') + self.xaxisRange.value = glm.dvec2(offset[0], offset[0] + extent[0]) + self.yaxisRange.value = glm.dvec2(offset[1], offset[1] + extent[1]) + self.zaxisRange.value = glm.dvec2(min_max[0], min_max[1]) m = ivw.glm.mat4(1.0) m[0][0] = extent[0] / mesh_extent[0] diff --git a/modules/astrophysics/python/processors/FitsVolumeSource.py b/modules/astrophysics/python/processors/FitsVolumeSource.py index ccad386ea..bc39de6c7 100644 --- a/modules/astrophysics/python/processors/FitsVolumeSource.py +++ b/modules/astrophysics/python/processors/FitsVolumeSource.py @@ -71,14 +71,20 @@ def __init__(self, id, name): ivw.md2doc("Scaling factor applied to the z axis"), 1.0, increment=0.001, min=(1.0e-8, cb.Immutable), max=(10.0, cb.Ignore)) - self.zaxisRange = ivw.properties.DoubleVec2Property( - "zaxisRange", "Z Axis Range", - ivw.md2doc("Extent of the dataset along the z axis."), - glm.dvec2(0.0, 1.0), increment=glm.dvec2(0.001), - min=(glm.dvec2(-1.70e308), cb.Ignore), - max=(glm.dvec2(1.79e308), cb.Ignore), - semantics=semantics.Text, invalidationLevel=inv.Valid) - self.zaxisRange.readOnly = True + def axisRangeProperty(axis_name: str) -> ivw.properties.DoubleVec2Property: + prop = ivw.properties.DoubleVec2Property( + f"{axis_name.lower()}axisRange", f"{axis_name.upper()} Axis Range", + ivw.md2doc(f"Extent of the dataset along the {axis_name.lower()} axis."), + glm.dvec2(0.0, 1.0), increment=glm.dvec2(0.001), + min=(glm.dvec2(-1.79e308), cb.Ignore), + max=(glm.dvec2(1.79e308), cb.Ignore), + semantics=semantics.Text, invalidationLevel=inv.Valid) + prop.readOnly = True + return prop + + self.xaxisRange = axisRangeProperty("x") + self.yaxisRange = axisRangeProperty("y") + self.zaxisRange = axisRangeProperty("z") self.fitsParameters = ivwastrovis.createFitsCompositeProperty( "fitsParameters", "Fits Parameters") @@ -88,8 +94,9 @@ def __init__(self, id, name): "information", "Data information") self.addProperties([self.filename, self.fitsParameters, - self.latLonCoords, self.axisType, self.zaxisRange, self.axisScaling, - self.basis, self.information]) + self.latLonCoords, self.axisType, + self.xaxisRange, self.yaxisRange, self.zaxisRange, + self.axisScaling, self.basis, self.information]) @staticmethod def processorInfo(): @@ -132,15 +139,15 @@ def process(self): if not filename.exists(): return - astroviscommon.readFITSData.cache_clear() + # astroviscommon.readFITSData.cache_clear() t = time.perf_counter_ns() fits_data: astroviscommon.FitsData = astroviscommon.readFITSData(filename) dim = fits_data.data.shape elapsed = time.perf_counter_ns() - t - print(f'{elapsed/1000.0} us') + print(f'loading FITS data: {elapsed/1_000_000.0} ms') - print(astroviscommon.readFITSData.cache_info()) + # print(astroviscommon.readFITSData.cache_info()) vel_func: Callable[[float], int] = astroviscommon.frequencyToVelocityFunc(fits_data.params) velocity: list = [vel_func(channel) for channel in range(dim[0])] @@ -199,7 +206,8 @@ def process(self): volume_intensity.worldMatrix = ivw.glm.mat4(w) # print(w * m, extent_z * (self.axisScaling.value / extent_z)) - + self.xaxisRange.value = glm.dvec2(offset_xy[0], offset_xy[0] + extent_xy[0]) + self.yaxisRange.value = glm.dvec2(offset_xy[1], offset_xy[1] + extent_xy[1]) self.zaxisRange.value = glm.dvec2(offset_z, offset_z + extent_z) self.fitsParameters.objectName.value = fits_data.params.objectName From b72f0a64ae6b83c646275dcfee06442592614088 Mon Sep 17 00:00:00 2001 From: Martin Falk Date: Tue, 31 Mar 2026 10:01:59 +0200 Subject: [PATCH 06/16] renamed modules folder --- {modules => infravis}/astrophysics/CMakeLists.txt | 0 {modules => infravis}/astrophysics/depends.cmake | 0 .../include/infravis/astrophysics/astrophysicsmodule.h | 0 .../include/infravis/astrophysics/astrophysicsmoduledefine.h | 0 .../infravis/astrophysics/processors/fitsvolumeraycaster.h | 0 .../shadercomponents/fitsnonlineardepthcomponent.h | 0 {modules => infravis}/astrophysics/python/astroviscommon.py | 0 {modules => infravis}/astrophysics/python/ivwastrovis.py | 0 .../astrophysics/python/processors/AngleToVector2D.py | 0 .../astrophysics/python/processors/AxisSliceAnnotation.py | 0 .../astrophysics/python/processors/FitsDepthCorrection.py | 0 .../astrophysics/python/processors/FitsPolarizationSource.py | 0 .../astrophysics/python/processors/FitsVolumeSource.py | 0 infravis/astrophysics/readme.md | 3 +++ {modules => infravis}/astrophysics/src/astrophysicsmodule.cpp | 0 .../astrophysics/src/processors/fitsvolumeraycaster.cpp | 0 .../src/shadercomponents/fitsnonlineardepthcomponent.cpp | 0 .../tests/unittests/astrophysics-unittest-main.cpp | 0 modules/astrophysics/readme.md | 3 --- 19 files changed, 3 insertions(+), 3 deletions(-) rename {modules => infravis}/astrophysics/CMakeLists.txt (100%) rename {modules => infravis}/astrophysics/depends.cmake (100%) rename {modules => infravis}/astrophysics/include/infravis/astrophysics/astrophysicsmodule.h (100%) rename {modules => infravis}/astrophysics/include/infravis/astrophysics/astrophysicsmoduledefine.h (100%) rename {modules => infravis}/astrophysics/include/infravis/astrophysics/processors/fitsvolumeraycaster.h (100%) rename {modules => infravis}/astrophysics/include/infravis/astrophysics/shadercomponents/fitsnonlineardepthcomponent.h (100%) rename {modules => infravis}/astrophysics/python/astroviscommon.py (100%) rename {modules => infravis}/astrophysics/python/ivwastrovis.py (100%) rename {modules => infravis}/astrophysics/python/processors/AngleToVector2D.py (100%) rename {modules => infravis}/astrophysics/python/processors/AxisSliceAnnotation.py (100%) rename {modules => infravis}/astrophysics/python/processors/FitsDepthCorrection.py (100%) rename {modules => infravis}/astrophysics/python/processors/FitsPolarizationSource.py (100%) rename {modules => infravis}/astrophysics/python/processors/FitsVolumeSource.py (100%) create mode 100644 infravis/astrophysics/readme.md rename {modules => infravis}/astrophysics/src/astrophysicsmodule.cpp (100%) rename {modules => infravis}/astrophysics/src/processors/fitsvolumeraycaster.cpp (100%) rename {modules => infravis}/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp (100%) rename {modules => infravis}/astrophysics/tests/unittests/astrophysics-unittest-main.cpp (100%) delete mode 100644 modules/astrophysics/readme.md diff --git a/modules/astrophysics/CMakeLists.txt b/infravis/astrophysics/CMakeLists.txt similarity index 100% rename from modules/astrophysics/CMakeLists.txt rename to infravis/astrophysics/CMakeLists.txt diff --git a/modules/astrophysics/depends.cmake b/infravis/astrophysics/depends.cmake similarity index 100% rename from modules/astrophysics/depends.cmake rename to infravis/astrophysics/depends.cmake diff --git a/modules/astrophysics/include/infravis/astrophysics/astrophysicsmodule.h b/infravis/astrophysics/include/infravis/astrophysics/astrophysicsmodule.h similarity index 100% rename from modules/astrophysics/include/infravis/astrophysics/astrophysicsmodule.h rename to infravis/astrophysics/include/infravis/astrophysics/astrophysicsmodule.h diff --git a/modules/astrophysics/include/infravis/astrophysics/astrophysicsmoduledefine.h b/infravis/astrophysics/include/infravis/astrophysics/astrophysicsmoduledefine.h similarity index 100% rename from modules/astrophysics/include/infravis/astrophysics/astrophysicsmoduledefine.h rename to infravis/astrophysics/include/infravis/astrophysics/astrophysicsmoduledefine.h diff --git a/modules/astrophysics/include/infravis/astrophysics/processors/fitsvolumeraycaster.h b/infravis/astrophysics/include/infravis/astrophysics/processors/fitsvolumeraycaster.h similarity index 100% rename from modules/astrophysics/include/infravis/astrophysics/processors/fitsvolumeraycaster.h rename to infravis/astrophysics/include/infravis/astrophysics/processors/fitsvolumeraycaster.h diff --git a/modules/astrophysics/include/infravis/astrophysics/shadercomponents/fitsnonlineardepthcomponent.h b/infravis/astrophysics/include/infravis/astrophysics/shadercomponents/fitsnonlineardepthcomponent.h similarity index 100% rename from modules/astrophysics/include/infravis/astrophysics/shadercomponents/fitsnonlineardepthcomponent.h rename to infravis/astrophysics/include/infravis/astrophysics/shadercomponents/fitsnonlineardepthcomponent.h diff --git a/modules/astrophysics/python/astroviscommon.py b/infravis/astrophysics/python/astroviscommon.py similarity index 100% rename from modules/astrophysics/python/astroviscommon.py rename to infravis/astrophysics/python/astroviscommon.py diff --git a/modules/astrophysics/python/ivwastrovis.py b/infravis/astrophysics/python/ivwastrovis.py similarity index 100% rename from modules/astrophysics/python/ivwastrovis.py rename to infravis/astrophysics/python/ivwastrovis.py diff --git a/modules/astrophysics/python/processors/AngleToVector2D.py b/infravis/astrophysics/python/processors/AngleToVector2D.py similarity index 100% rename from modules/astrophysics/python/processors/AngleToVector2D.py rename to infravis/astrophysics/python/processors/AngleToVector2D.py diff --git a/modules/astrophysics/python/processors/AxisSliceAnnotation.py b/infravis/astrophysics/python/processors/AxisSliceAnnotation.py similarity index 100% rename from modules/astrophysics/python/processors/AxisSliceAnnotation.py rename to infravis/astrophysics/python/processors/AxisSliceAnnotation.py diff --git a/modules/astrophysics/python/processors/FitsDepthCorrection.py b/infravis/astrophysics/python/processors/FitsDepthCorrection.py similarity index 100% rename from modules/astrophysics/python/processors/FitsDepthCorrection.py rename to infravis/astrophysics/python/processors/FitsDepthCorrection.py diff --git a/modules/astrophysics/python/processors/FitsPolarizationSource.py b/infravis/astrophysics/python/processors/FitsPolarizationSource.py similarity index 100% rename from modules/astrophysics/python/processors/FitsPolarizationSource.py rename to infravis/astrophysics/python/processors/FitsPolarizationSource.py diff --git a/modules/astrophysics/python/processors/FitsVolumeSource.py b/infravis/astrophysics/python/processors/FitsVolumeSource.py similarity index 100% rename from modules/astrophysics/python/processors/FitsVolumeSource.py rename to infravis/astrophysics/python/processors/FitsVolumeSource.py diff --git a/infravis/astrophysics/readme.md b/infravis/astrophysics/readme.md new file mode 100644 index 000000000..38e0ad10e --- /dev/null +++ b/infravis/astrophysics/readme.md @@ -0,0 +1,3 @@ +# AstroPhysics Module + +Supports importing of FITS datasets using the `astropy` Python module. Primary focus is data stemming from the ALMA radio telescope. The `FitsVolumeRaycaster` processor visualizes volumetric data while applying a depth corrections for the positions around a star. diff --git a/modules/astrophysics/src/astrophysicsmodule.cpp b/infravis/astrophysics/src/astrophysicsmodule.cpp similarity index 100% rename from modules/astrophysics/src/astrophysicsmodule.cpp rename to infravis/astrophysics/src/astrophysicsmodule.cpp diff --git a/modules/astrophysics/src/processors/fitsvolumeraycaster.cpp b/infravis/astrophysics/src/processors/fitsvolumeraycaster.cpp similarity index 100% rename from modules/astrophysics/src/processors/fitsvolumeraycaster.cpp rename to infravis/astrophysics/src/processors/fitsvolumeraycaster.cpp diff --git a/modules/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp b/infravis/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp similarity index 100% rename from modules/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp rename to infravis/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp diff --git a/modules/astrophysics/tests/unittests/astrophysics-unittest-main.cpp b/infravis/astrophysics/tests/unittests/astrophysics-unittest-main.cpp similarity index 100% rename from modules/astrophysics/tests/unittests/astrophysics-unittest-main.cpp rename to infravis/astrophysics/tests/unittests/astrophysics-unittest-main.cpp diff --git a/modules/astrophysics/readme.md b/modules/astrophysics/readme.md deleted file mode 100644 index 5cd47b350..000000000 --- a/modules/astrophysics/readme.md +++ /dev/null @@ -1,3 +0,0 @@ -# AstroPhysics Module - -Description of the AstroPhysics module From 063c8f2310ce44dbd850fc11a99daabd0fddd342 Mon Sep 17 00:00:00 2001 From: Martin Falk Date: Tue, 31 Mar 2026 10:02:11 +0200 Subject: [PATCH 07/16] Meta: updated readme --- README.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/README.md b/README.md index ff70a5c56..2fd586ddc 100644 --- a/README.md +++ b/README.md @@ -6,6 +6,11 @@ Additional Inviwo modules grouped into different categories. To enable the modul * IntegralLineFiltering: Provides functionality to convert a `IntegralLineSet` to `DataFrame` containing various metrics of the lines that can be used together with the plotting functionality in the `Plotting` and `PlottingGL` modules in core. See example workspace on how to use it. ## infovis - Information Visualization + +## infravis - various modules + +* AstroPhysics: primarily Python-based processors dealing with ALMA data in the FITS data format + ## medvis - Medical Visualization * DICOM: functionality for loading and writing DICOM files using the Grassroots DICOM ([GDCM](https://sourceforge.net/projects/gdcm/)) From 04b288f7532932ba97567eda5d8049f5104e9a70 Mon Sep 17 00:00:00 2001 From: Martin Falk Date: Tue, 31 Mar 2026 10:27:14 +0200 Subject: [PATCH 08/16] gha: astophysics build --- .github/presets/linux.json | 30 ++++++++++++++++++++++++++---- 1 file changed, 26 insertions(+), 4 deletions(-) diff --git a/.github/presets/linux.json b/.github/presets/linux.json index 83abd6826..dc1836c31 100644 --- a/.github/presets/linux.json +++ b/.github/presets/linux.json @@ -13,10 +13,32 @@ "CMAKE_C_COMPILER_LAUNCHER": "ccache", "CMAKE_CXX_COMPILER_LAUNCHER": "ccache", - "IVW_MODULE_FFMPEG": "OFF", - "IVW_MODULE_WEBBROWSER": "OFF", - "IVW_MODULE_WEBQT": "OFF", - "IVW_TEST_BENCHMARKS": "OFF" + "IVW_EXTERNAL_MODULES" : "${sourceParentDir}/modules/misc;${sourceParentDir}/modules/medvis;${sourceParentDir}/modules/molvis;${sourceParentDir}/modules/tensorvis;${sourceParentDir}/modules/topovis;${sourceParentDir}/modules/vectorvis;${sourceParentDir}/modules/infravis", + + "IVW_TEST_BENCHMARKS": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_ASTROPHYSICS": { "type": "BOOL", "value": "ON" }, + "IVW_MODULE_FFMPEG": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_WEBBROWSER": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_WEBQT": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_VTK": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_COMPUTESHADEREXAMPLES": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_DEVTOOLS": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_GRAPHVIZ": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_OPENMESH": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_SPRINGSYSTEM": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_VASP": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_GAUSSIAN": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_NANOVGUTILS": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_TENSORVISBASE": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_TENSORVISIO": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_TENSORVISPYTHON": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_INTEGRALLINEFILTERING": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_MOLVISBASE": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_MOLVISGL": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_MOLVISPYTHON": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_DATAFRAMECLUSTERING": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_C3D": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_TTK": { "type": "BOOL", "value": "OFF" } } }, { From 597af0e23c7e369f6eb07c7a9e6add0a18d475a1 Mon Sep 17 00:00:00 2001 From: Martin Falk Date: Tue, 31 Mar 2026 14:22:46 +0200 Subject: [PATCH 09/16] AstroPhysics: packaging --- infravis/astrophysics/CMakeLists.txt | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/infravis/astrophysics/CMakeLists.txt b/infravis/astrophysics/CMakeLists.txt index 9cc54cf32..1c764ab2e 100644 --- a/infravis/astrophysics/CMakeLists.txt +++ b/infravis/astrophysics/CMakeLists.txt @@ -16,9 +16,13 @@ set(SOURCE_FILES ivw_group("Source Files" ${SOURCE_FILES}) set(PYTHON_FILES + python/astroviscommon.py + python/ivwastrovis.py + python/processors/AngleToVector2D.py + python/processors/AxisSliceAnnotation.py + python/processors/FitsDepthCorrection.py python/processors/FitsPolarizationSource.py python/processors/FitsVolumeSource.py - python/astroviscommon.py ) ivw_group("Python Files" ${PYTHON_FILES}) @@ -34,6 +38,8 @@ ivw_add_unittest(${TEST_FILES}) ivw_create_module(${SOURCE_FILES} ${HEADER_FILES} ${SHADER_FILES} ${PYTHON_FILES}) +ivw_add_to_module_pack(${CMAKE_CURRENT_SOURCE_DIR}/python) + set_property(GLOBAL APPEND PROPERTY IVW_PYTHON_ADDITIONAL_MODULES astropy) # Add shader directory to install package From 5e66e582f9fe812efa824b7971db9e565ca26a52 Mon Sep 17 00:00:00 2001 From: Martin Falk Date: Tue, 31 Mar 2026 14:58:24 +0200 Subject: [PATCH 10/16] AstroPhysics: linter --- .../astrophysics/processors/fitsvolumeraycaster.h | 4 ++-- .../astrophysics/src/processors/fitsvolumeraycaster.cpp | 6 ++---- .../src/shadercomponents/fitsnonlineardepthcomponent.cpp | 1 - .../tests/unittests/astrophysics-unittest-main.cpp | 9 ++++----- 4 files changed, 8 insertions(+), 12 deletions(-) diff --git a/infravis/astrophysics/include/infravis/astrophysics/processors/fitsvolumeraycaster.h b/infravis/astrophysics/include/infravis/astrophysics/processors/fitsvolumeraycaster.h index 9349aaef0..3d0091903 100644 --- a/infravis/astrophysics/include/infravis/astrophysics/processors/fitsvolumeraycaster.h +++ b/infravis/astrophysics/include/infravis/astrophysics/processors/fitsvolumeraycaster.h @@ -46,7 +46,8 @@ namespace inviwo { class IVW_MODULE_ASTROPHYSICS_API FitsVolumeRaycaster : public VolumeRaycasterBase { public: - explicit FitsVolumeRaycaster(std::string_view identifier = {}, std::string_view displayName = {}); + explicit FitsVolumeRaycaster(std::string_view identifier = {}, + std::string_view displayName = {}); virtual void process() override; @@ -55,7 +56,6 @@ class IVW_MODULE_ASTROPHYSICS_API FitsVolumeRaycaster : public VolumeRaycasterBa private: FitsNonlinearDepthComponent fits_; - //VolumeComponent volume_; EntryExitComponent entryExit_; BackgroundComponent background_; IsoTFComponent<1> isoTF_; diff --git a/infravis/astrophysics/src/processors/fitsvolumeraycaster.cpp b/infravis/astrophysics/src/processors/fitsvolumeraycaster.cpp index 3b5cb464f..373b7b6f3 100644 --- a/infravis/astrophysics/src/processors/fitsvolumeraycaster.cpp +++ b/infravis/astrophysics/src/processors/fitsvolumeraycaster.cpp @@ -57,11 +57,9 @@ FitsVolumeRaycaster::FitsVolumeRaycaster(std::string_view identifier, std::strin , isoTF_{fits_.volumePort_} , raycasting_{fits_.getName(), isoTF_.isotfs[0]} , camera_{"camera", util::boundingBox(fits_.volumePort_)} - , light_{&camera_.camera} -{ + , light_{&camera_.camera} { - registerComponents(fits_, entryExit_, background_, raycasting_, isoTF_, camera_, - light_); + registerComponents(fits_, entryExit_, background_, raycasting_, isoTF_, camera_, light_); } void FitsVolumeRaycaster::process() { diff --git a/infravis/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp b/infravis/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp index 6aebaf923..e94b52100 100644 --- a/infravis/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp +++ b/infravis/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp @@ -127,7 +127,6 @@ void FitsNonlinearDepthComponent::process(Shader& shader, TextureUnitContainer& utilgl::bindAndSetUniforms(shader, cont, volumePort_); const auto data = volumePort_.getData(); - const vec3 delta = dataExtent_.get() / vec3{data->getDimensions()}; mat4 m = glm::scale(dataExtent_.get()); m[3] = vec4{dataOffset_.get(), 1.0f}; diff --git a/infravis/astrophysics/tests/unittests/astrophysics-unittest-main.cpp b/infravis/astrophysics/tests/unittests/astrophysics-unittest-main.cpp index acb90bbf0..8ddf58e4e 100644 --- a/infravis/astrophysics/tests/unittests/astrophysics-unittest-main.cpp +++ b/infravis/astrophysics/tests/unittests/astrophysics-unittest-main.cpp @@ -41,11 +41,10 @@ #include int main(int argc, char** argv) { - using namespace inviwo; - LogCentral::init(); - auto logger = std::make_shared(); - LogCentral::getPtr()->setVerbosity(LogVerbosity::Error); - LogCentral::getPtr()->registerLogger(logger); + inviwo::LogCentral::init(); + auto logger = std::make_shared(); + inviwo::LogCentral::getPtr()->setVerbosity(inviwo::LogVerbosity::Error); + inviwo::LogCentral::getPtr()->registerLogger(logger); int ret = -1; { From 69a1b5454e056c21b740d047060fa47dd26fba60 Mon Sep 17 00:00:00 2001 From: Martin Falk Date: Wed, 1 Apr 2026 12:02:09 +0200 Subject: [PATCH 11/16] gha: packaging test --- .github/presets/linux.json | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/.github/presets/linux.json b/.github/presets/linux.json index dc1836c31..23e2a0905 100644 --- a/.github/presets/linux.json +++ b/.github/presets/linux.json @@ -38,7 +38,9 @@ "IVW_MODULE_MOLVISPYTHON": { "type": "BOOL", "value": "OFF" }, "IVW_MODULE_DATAFRAMECLUSTERING": { "type": "BOOL", "value": "OFF" }, "IVW_MODULE_C3D": { "type": "BOOL", "value": "OFF" }, - "IVW_MODULE_TTK": { "type": "BOOL", "value": "OFF" } + "IVW_MODULE_TTK": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_OPACTOPT": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_PYTHONTOOLS": { "type": "BOOL", "value": "OFF" } } }, { From f03ad48cf1d0d06ee47198bfbcd8e26ae2d84eb1 Mon Sep 17 00:00:00 2001 From: Martin Falk Date: Wed, 1 Apr 2026 12:25:59 +0200 Subject: [PATCH 12/16] AstroPhysics: minor changes --- .../include/infravis/astrophysics/astrophysicsmodule.h | 1 - infravis/astrophysics/python/ivwastrovis.py | 2 -- infravis/astrophysics/python/processors/AngleToVector2D.py | 2 +- infravis/astrophysics/python/processors/FitsVolumeSource.py | 2 +- 4 files changed, 2 insertions(+), 5 deletions(-) diff --git a/infravis/astrophysics/include/infravis/astrophysics/astrophysicsmodule.h b/infravis/astrophysics/include/infravis/astrophysics/astrophysicsmodule.h index 7d9813fed..83b864b1d 100644 --- a/infravis/astrophysics/include/infravis/astrophysics/astrophysicsmodule.h +++ b/infravis/astrophysics/include/infravis/astrophysics/astrophysicsmodule.h @@ -38,7 +38,6 @@ namespace inviwo { class IVW_MODULE_ASTROPHYSICS_API AstroPhysicsModule : public InviwoModule { public: explicit AstroPhysicsModule(InviwoApplication* app); - virtual ~AstroPhysicsModule() = default; pyutil::ModulePath scripts_; PythonProcessorFolderObserver pythonFolderObserver_; diff --git a/infravis/astrophysics/python/ivwastrovis.py b/infravis/astrophysics/python/ivwastrovis.py index 9cf7971c8..51d83c084 100644 --- a/infravis/astrophysics/python/ivwastrovis.py +++ b/infravis/astrophysics/python/ivwastrovis.py @@ -78,13 +78,11 @@ def createDepthMesh(params: FitsParams, positions_np = np.array(positions) min_max = [np.min(positions_np[:, 2]), np.max(positions_np[:, 2])] - # ivw.logInfo(f'min/max: {min_max}') depth_scaling = np.max(np.abs(min_max)) colors = [tf.sample(value / depth_scaling * 0.5 + 0.5) for value in positions_np[:, 2]] extent, offset, _ = astroviscommon.getLatLongBasis(params, lat_long) - ivw.logInfo(f'{extent=}\n{offset=}') positions_np = (positions_np * np.array([extent[0] / (dims[0] - 1), extent[1] / (dims[1] - 1), 1.0]) + np.array([offset[0], offset[1], 0.0])) diff --git a/infravis/astrophysics/python/processors/AngleToVector2D.py b/infravis/astrophysics/python/processors/AngleToVector2D.py index 2e22bed3d..33a14126b 100644 --- a/infravis/astrophysics/python/processors/AngleToVector2D.py +++ b/infravis/astrophysics/python/processors/AngleToVector2D.py @@ -47,7 +47,7 @@ def process(self): magnitude = volume_magnitude.data if use_degree: - angle *= np.pi / 180.0 + angle = angle * np.pi / 180.0 data = np.stack([np.cos(angle + np.pi / 2.0) * magnitude, np.sin(angle + np.pi / 2.0) * magnitude], axis=3) volumerep = ivw.data.VolumePy(np.copy(data.astype(dtype=np.float32), order='C')) diff --git a/infravis/astrophysics/python/processors/FitsVolumeSource.py b/infravis/astrophysics/python/processors/FitsVolumeSource.py index bc39de6c7..ccc7ea930 100644 --- a/infravis/astrophysics/python/processors/FitsVolumeSource.py +++ b/infravis/astrophysics/python/processors/FitsVolumeSource.py @@ -137,7 +137,7 @@ def process(self): filename: Path = ivwbase.io.downloadAndCacheIfUrl(self.filename.valueAsString()) if not filename.exists(): - return + raise FileNotFoundError(f"Could not locate {self.filename.valueAsString():s}") # astroviscommon.readFITSData.cache_clear() From 6394ed3bf135c637d7a761653d7f4fe4eafe3123 Mon Sep 17 00:00:00 2001 From: Martin Falk Date: Thu, 9 Apr 2026 17:21:10 +0200 Subject: [PATCH 13/16] AstroPhysics: minor fixes --- .../python/processors/FitsPolarizationSource.py | 4 ++-- infravis/astrophysics/python/processors/FitsVolumeSource.py | 6 ++++++ 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/infravis/astrophysics/python/processors/FitsPolarizationSource.py b/infravis/astrophysics/python/processors/FitsPolarizationSource.py index d2652035c..395ca7182 100644 --- a/infravis/astrophysics/python/processors/FitsPolarizationSource.py +++ b/infravis/astrophysics/python/processors/FitsPolarizationSource.py @@ -69,8 +69,8 @@ def createVectorField(angle: np.array, length: np.array, volume = ivw.data.Volume(volumerep) # TODO: extract extent from hdul[0].header - volume.basis = ivw.glm.mat3(1) - volume.offset = ivw.glm.vec3(-0.5, -0.5, -0.5) + volume.basis = ivw.glm.dmat3(1) + volume.offset = ivw.glm.dvec3(-0.5, -0.5, -0.5) max_len: float = np.nanmax(length) datarange = ivw.glm.dvec2(-max_len, max_len) diff --git a/infravis/astrophysics/python/processors/FitsVolumeSource.py b/infravis/astrophysics/python/processors/FitsVolumeSource.py index ccc7ea930..8e68706f2 100644 --- a/infravis/astrophysics/python/processors/FitsVolumeSource.py +++ b/infravis/astrophysics/python/processors/FitsVolumeSource.py @@ -182,6 +182,12 @@ def process(self): volume_intensity = FitsVolumeSource.createVolume(fits_data) volume_intensity.dataMap.valueAxis.name = fits_data.params.btype + if fits_data.params.bunit in ['deg', 'rad'] and fits_data.params.btype != 'Angle': + # overriding btype as it does not match the unit + volume_intensity.dataMap.valueAxis.name = 'Angle' + ivw.logWarn(f"{self.identifier}: " + f"Value axis set to 'Angle' since FITS btype '{fits_data.params.btype}' " + f"does not match bunit '{fits_data.params.bunit}'") # hdul[0].header['bunit'] is 'Jy/beam ', which is not a valid unit in Inviwo volume_intensity.dataMap.valueAxis.unit = ivw.data.Unit( fits_data.params.bunit.replace('/beam', '')) From 16a7a2e478ed0dd8c09c3dd3dfc2cd1c2e649a31 Mon Sep 17 00:00:00 2001 From: Martin Falk Date: Thu, 9 Apr 2026 18:24:56 +0200 Subject: [PATCH 14/16] AstroPhysics: added axis range scaling --- infravis/astrophysics/python/ivwastrovis.py | 22 +++++++++++++++++++ .../python/processors/FitsVolumeSource.py | 19 ++++++++++------ 2 files changed, 34 insertions(+), 7 deletions(-) diff --git a/infravis/astrophysics/python/ivwastrovis.py b/infravis/astrophysics/python/ivwastrovis.py index 51d83c084..afa00b996 100644 --- a/infravis/astrophysics/python/ivwastrovis.py +++ b/infravis/astrophysics/python/ivwastrovis.py @@ -6,6 +6,7 @@ from inviwopy import glm import numpy as np +import math from contextlib import suppress from collections.abc import Callable @@ -199,3 +200,24 @@ def createFitsCompositeProperty(identifier: str, crpix3, crval3, cdelt3, restFrequency, extent, offset, valueName, valueUnit, dataRange]) return prop + + +def scaleRange(r: glm.dvec2) -> tuple[glm.dvec2, int]: + """ + Determine an appropriate scaling factor for the input range based on thousand separators, + i.e. digits grouped by 3. + + Parameters + ---------- + r: glm.dvec2 + input range + + Returns + ------- + A tuple containing the scaled range and the exponent for base 10 + """ + s: int = 3 + exp: float = math.floor(math.log10(glm.compMax(glm.abs(r)))) + major_exp: float = math.floor(exp / s) * s if exp > 0 else math.round(exp / s) * s + factor: float = math.pow(10.0, -major_exp) + return (r * factor, int(major_exp)) diff --git a/infravis/astrophysics/python/processors/FitsVolumeSource.py b/infravis/astrophysics/python/processors/FitsVolumeSource.py index 8e68706f2..fe6beabbf 100644 --- a/infravis/astrophysics/python/processors/FitsVolumeSource.py +++ b/infravis/astrophysics/python/processors/FitsVolumeSource.py @@ -8,6 +8,7 @@ import ivwastrovis import numpy as np +import math from pathlib import Path from collections.abc import Callable from enum import Enum @@ -156,7 +157,7 @@ def process(self): extent_xy, offset_xy, unit = astroviscommon.getLatLongBasis( fits_data.params, LatLong(self.latLonCoords.value)) - axes = [("x", unit), ("y", unit), ("z", "")] + axes = [("x", ivw.data.Unit(unit)), ("y", ivw.data.Unit(unit)), ("z", ivw.data.Unit(""))] extent_z: float = 0.0 offset_z: float = 0.0 @@ -165,18 +166,22 @@ def process(self): case AxisType.SliceIndex: extent_z = float(dim[0] - 1) offset_z = 0 - axes[2] = ("Channel", "") + axes[2] = ("Channel", ivw.data.Unit("")) case AxisType.Dataset: freq_func = astroviscommon.frequencyFunc(fits_data.params) freq_range = (freq_func(0), freq_func(dim[0] - 1)) - extent_z = freq_range[1] - freq_range[0] - offset_z = freq_range[0] - axes[2] = (fits_data.params.ctype[2], fits_data.params.cunit[2]) + # rescale range, e.g. 345 GHz instead of 3.45e11 Hz + scaled_range, exponent = ivwastrovis.scaleRange(glm.dvec2(freq_range)) + extent_z = scaled_range.y - scaled_range.x + offset_z = scaled_range.x + axes[2] = (fits_data.params.ctype[2], + ivw.data.Unit(math.pow(10.0, exponent), + ivw.data.Unit(fits_data.params.cunit[2]))) case AxisType.Velocity: vel_range = (vel_func(0), vel_func(dim[0] - 1)) extent_z = vel_range[1] - vel_range[0] offset_z = vel_range[0] - axes[2] = ("Velocity", "km/s") + axes[2] = ("Velocity", ivw.data.Unit("km/s")) case _: raise ValueError(f"Unexpected axis type {self.axisType.value}") @@ -194,7 +199,7 @@ def process(self): for i, (label, unit) in enumerate(axes): volume_intensity.axes[i].name = label - volume_intensity.axes[i].unit = ivw.data.Unit(unit) + volume_intensity.axes[i].unit = unit m = ivw.glm.dmat4(1.0) m[0][0] = extent_xy[0] From 9186e48fa4add49a8a729c6e895bdf19cb7305ad Mon Sep 17 00:00:00 2001 From: Martin Falk Date: Fri, 10 Apr 2026 10:14:13 +0200 Subject: [PATCH 15/16] AstroPhysics: more robust data import --- .../astrophysics/python/astroviscommon.py | 26 ++++++++++++------- .../python/processors/FitsVolumeSource.py | 9 ++++--- 2 files changed, 21 insertions(+), 14 deletions(-) diff --git a/infravis/astrophysics/python/astroviscommon.py b/infravis/astrophysics/python/astroviscommon.py index 08d9b59a9..928be0232 100644 --- a/infravis/astrophysics/python/astroviscommon.py +++ b/infravis/astrophysics/python/astroviscommon.py @@ -59,19 +59,25 @@ def readFITSData(filepath: Path, def extractFITSParams(header: fits.header.Header) -> FitsParams: - rest_frequency: float = header['restfrq'] # equals 345_339_760_000 + def get_attrib(attrib: str, + default_value: str | int | float | None = None) -> str | int | float | None: + if attrib in header: + return header[attrib] + return default_value + + rest_frequency: float = get_attrib('restfrq', 1.0) # equals 345_339_760_000 # rest_frequency = 345_339_769_300 # [hz] return FitsParams( - objectName=str(header['object']), - naxis=tuple(header[f'naxis{i}'] for i in range(1, header['naxis'] + 1)), - ctype=tuple(header[f'ctype{i}'].lower() for i in range(1, 4)), - cunit=tuple(header[f'cunit{i}'] for i in range(1, 4)), - crpix=tuple(header[f'crpix{i}'] for i in range(1, 4)), - crval=tuple(header[f'crval{i}'] for i in range(1, 4)), - cdelt=tuple(header[f'cdelt{i}'] for i in range(1, 4)), - btype=header['btype'], - bunit=header['bunit'], + objectName=str(get_attrib('object')), + naxis=tuple(get_attrib(f'naxis{i}', 0) for i in range(1, get_attrib('naxis', 0) + 1)), + ctype=tuple(get_attrib(f'ctype{i}', '').lower() for i in range(1, 4)), + cunit=tuple(get_attrib(f'cunit{i}') for i in range(1, 4)), + crpix=tuple(get_attrib(f'crpix{i}', 0.0) for i in range(1, 4)), + crval=tuple(get_attrib(f'crval{i}', 0.0) for i in range(1, 4)), + cdelt=tuple(get_attrib(f'cdelt{i}', 0.0) for i in range(1, 4)), + btype=get_attrib('btype', ''), + bunit=get_attrib('bunit', ''), rest_frequency=rest_frequency) diff --git a/infravis/astrophysics/python/processors/FitsVolumeSource.py b/infravis/astrophysics/python/processors/FitsVolumeSource.py index fe6beabbf..d9947570f 100644 --- a/infravis/astrophysics/python/processors/FitsVolumeSource.py +++ b/infravis/astrophysics/python/processors/FitsVolumeSource.py @@ -150,10 +150,6 @@ def process(self): # print(astroviscommon.readFITSData.cache_info()) - vel_func: Callable[[float], int] = astroviscommon.frequencyToVelocityFunc(fits_data.params) - velocity: list = [vel_func(channel) for channel in range(dim[0])] - print(f'Velocity matching frequency: {velocity}') - extent_xy, offset_xy, unit = astroviscommon.getLatLongBasis( fits_data.params, LatLong(self.latLonCoords.value)) @@ -178,6 +174,11 @@ def process(self): ivw.data.Unit(math.pow(10.0, exponent), ivw.data.Unit(fits_data.params.cunit[2]))) case AxisType.Velocity: + vel_func: Callable[[float], int] = astroviscommon.frequencyToVelocityFunc( + fits_data.params) + velocity: list = [vel_func(channel) for channel in range(dim[0])] + print(f'Velocity matching frequency: {velocity}') + vel_range = (vel_func(0), vel_func(dim[0] - 1)) extent_z = vel_range[1] - vel_range[0] offset_z = vel_range[0] From 48a128eaaba11566a4b0c085a00068eb323a2c8f Mon Sep 17 00:00:00 2001 From: Martin Falk Date: Fri, 10 Apr 2026 15:25:24 +0200 Subject: [PATCH 16/16] AstroPhysics: optional port fix --- .../python/processors/AngleToVector2D.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/infravis/astrophysics/python/processors/AngleToVector2D.py b/infravis/astrophysics/python/processors/AngleToVector2D.py index 33a14126b..ba048926d 100644 --- a/infravis/astrophysics/python/processors/AngleToVector2D.py +++ b/infravis/astrophysics/python/processors/AngleToVector2D.py @@ -36,7 +36,7 @@ def initializeResources(self): def process(self): volume_angle = self.inports.angle.getData() - volume_magnitude = self.inports.magnitude.getData() + volume_magnitude = self.inports.magnitude.getData() if self.inports.magnitude.hasData() else None use_degree: bool = True # assume degree if anything other than radians is stated as unit @@ -44,7 +44,7 @@ def process(self): use_degree = False angle = volume_angle.data - magnitude = volume_magnitude.data + magnitude = volume_magnitude.data if volume_magnitude else 1.0 if use_degree: angle = angle * np.pi / 180.0 @@ -58,10 +58,13 @@ def process(self): volume.dataMap.dataRange = datarange volume.dataMap.valueRange = datarange - volume.dataMap.valueAxis = volume_magnitude.dataMap.valueAxis - volume.modelMatrix = volume_magnitude.modelMatrix - volume.worldMatrix = volume_magnitude.worldMatrix - volume.axes = volume_magnitude.axes - volume.copyMetaDataFrom(volume_magnitude) + volume.modelMatrix = volume_angle.modelMatrix + volume.worldMatrix = volume_angle.worldMatrix + volume.axes = volume_angle.axes + if volume_magnitude: + volume.dataMap.valueAxis = volume_magnitude.dataMap.valueAxis + volume.copyMetaDataFrom(volume_magnitude) + else: + volume.copyMetaDataFrom(volume_angle) self.outports.outport.setData(volume)