diff --git a/.github/presets/linux.json b/.github/presets/linux.json index 83abd6826..23e2a0905 100644 --- a/.github/presets/linux.json +++ b/.github/presets/linux.json @@ -13,10 +13,34 @@ "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" }, + "IVW_MODULE_OPACTOPT": { "type": "BOOL", "value": "OFF" }, + "IVW_MODULE_PYTHONTOOLS": { "type": "BOOL", "value": "OFF" } } }, { 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/)) diff --git a/infravis/astrophysics/CMakeLists.txt b/infravis/astrophysics/CMakeLists.txt new file mode 100644 index 000000000..1c764ab2e --- /dev/null +++ b/infravis/astrophysics/CMakeLists.txt @@ -0,0 +1,46 @@ +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/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 +) +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}) + +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 +#ivw_add_to_module_pack(${CMAKE_CURRENT_SOURCE_DIR}/glsl) diff --git a/infravis/astrophysics/depends.cmake b/infravis/astrophysics/depends.cmake new file mode 100644 index 000000000..5c984fd76 --- /dev/null +++ b/infravis/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/infravis/astrophysics/include/infravis/astrophysics/astrophysicsmodule.h b/infravis/astrophysics/include/infravis/astrophysics/astrophysicsmodule.h new file mode 100644 index 000000000..83b864b1d --- /dev/null +++ b/infravis/astrophysics/include/infravis/astrophysics/astrophysicsmodule.h @@ -0,0 +1,46 @@ +/********************************************************************************* + * + * 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); + + pyutil::ModulePath scripts_; + PythonProcessorFolderObserver pythonFolderObserver_; +}; + +} // namespace inviwo diff --git a/infravis/astrophysics/include/infravis/astrophysics/astrophysicsmoduledefine.h b/infravis/astrophysics/include/infravis/astrophysics/astrophysicsmoduledefine.h new file mode 100644 index 000000000..1c54b3532 --- /dev/null +++ b/infravis/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/infravis/astrophysics/include/infravis/astrophysics/processors/fitsvolumeraycaster.h b/infravis/astrophysics/include/infravis/astrophysics/processors/fitsvolumeraycaster.h new file mode 100644 index 000000000..3d0091903 --- /dev/null +++ b/infravis/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_; + EntryExitComponent entryExit_; + BackgroundComponent background_; + IsoTFComponent<1> isoTF_; + RaycastingComponent raycasting_; + CameraComponent camera_; + LightComponent light_; +}; + +} // namespace inviwo diff --git a/infravis/astrophysics/include/infravis/astrophysics/shadercomponents/fitsnonlineardepthcomponent.h b/infravis/astrophysics/include/infravis/astrophysics/shadercomponents/fitsnonlineardepthcomponent.h new file mode 100644 index 000000000..76f742763 --- /dev/null +++ b/infravis/astrophysics/include/infravis/astrophysics/shadercomponents/fitsnonlineardepthcomponent.h @@ -0,0 +1,79 @@ +/********************************************************************************* + * + * 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_; + StringProperty objectName_; + FloatProperty assumedDistance_; + FloatProperty velocityInf_; + FloatProperty velocityStar_; + FloatProperty crpix3_; + FloatProperty crval3_; + FloatProperty cdelt3_; + FloatProperty restFrequency_; + FloatVec3Property dataExtent_; + FloatVec3Property dataOffset_; + IntVec2Property validChannelRange_; + FloatVec3Property invalidColor_; + FloatProperty invalidAlpha_; +}; + +} // namespace inviwo diff --git a/infravis/astrophysics/python/astroviscommon.py b/infravis/astrophysics/python/astroviscommon.py new file mode 100644 index 000000000..928be0232 --- /dev/null +++ b/infravis/astrophysics/python/astroviscommon.py @@ -0,0 +1,233 @@ + +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 + 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] + crpix: tuple[float, float, float] + crval: tuple[float, float, float] + cdelt: tuple[float, float, float] + btype: str + bunit: str + rest_frequency: float # [Hz] + + +@dataclass +class FitsData: + params: FitsParams + data: np.ndarray + + +@lru_cache(maxsize=5) +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: + 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(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) + + +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]]: + 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, str]: + """ + 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, str] + lat/long extent and offset, and unit + + """ + 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 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 == 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, getBasisUnit(lat_long)) diff --git a/infravis/astrophysics/python/ivwastrovis.py b/infravis/astrophysics/python/ivwastrovis.py new file mode 100644 index 000000000..afa00b996 --- /dev/null +++ b/infravis/astrophysics/python/ivwastrovis.py @@ -0,0 +1,223 @@ +import astroviscommon +from astroviscommon import FitsParams, LatLong + +import inviwopy as ivw +import ivwbase +from inviwopy import glm + +import numpy as np +import math +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, (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) + + dims = np.array(params.naxis) + + 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 valid_channels: + 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]), (valid_channels[0], valid_channels[-1])) + + positions_np = np.array(positions) + min_max = [np.min(positions_np[:, 2]), np.max(positions_np[:, 2])] + 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) + 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[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), (valid_channels[0], valid_channels[-1])) + + +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 + + +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/AngleToVector2D.py b/infravis/astrophysics/python/processors/AngleToVector2D.py new file mode 100644 index 000000000..ba048926d --- /dev/null +++ b/infravis/astrophysics/python/processors/AngleToVector2D.py @@ -0,0 +1,70 @@ +# 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() if self.inports.magnitude.hasData() else None + + 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 volume_magnitude else 1.0 + + if use_degree: + 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')) + 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.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) diff --git a/infravis/astrophysics/python/processors/AxisSliceAnnotation.py b/infravis/astrophysics/python/processors/AxisSliceAnnotation.py new file mode 100644 index 000000000..58aa9935f --- /dev/null +++ b/infravis/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] - 1) + 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/infravis/astrophysics/python/processors/FitsDepthCorrection.py b/infravis/astrophysics/python/processors/FitsDepthCorrection.py new file mode 100644 index 000000000..adea554dd --- /dev/null +++ b/infravis/astrophysics/python/processors/FitsDepthCorrection.py @@ -0,0 +1,157 @@ +# 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) + + cb = ivw.properties.ConstraintBehavior + IntOption = ivw.properties.IntOption + semantics = ivw.properties.PropertySemantics + + self.addInport(ivw.PythonInport("fitsdata")) + self.addOutport(ivw.data.MeshOutport("depthmesh")) + + 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"), + 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)) + + 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.xaxisRange, self.yaxisRange, self.zaxisRange, + self.axisScaling, 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() + latlong_calc: LatLong = LatLong.RelativeArcsec + depth_mesh, min_max, valid_channel_range = ivwastrovis.createDepthMesh( + fits_data.params, self.tf.value, + 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'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) + 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(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) + + 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] + 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] = [- 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 + + self.outports.depthmesh.setData(depth_mesh) diff --git a/infravis/astrophysics/python/processors/FitsPolarizationSource.py b/infravis/astrophysics/python/processors/FitsPolarizationSource.py new file mode 100644 index 000000000..395ca7182 --- /dev/null +++ b/infravis/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.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) + 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/infravis/astrophysics/python/processors/FitsVolumeSource.py b/infravis/astrophysics/python/processors/FitsVolumeSource.py new file mode 100644 index 000000000..d9947570f --- /dev/null +++ b/infravis/astrophysics/python/processors/FitsVolumeSource.py @@ -0,0 +1,259 @@ +# Name: FitsVolumeSource + +import inviwopy as ivw +from inviwopy import glm +import ivwbase +import astroviscommon +from astroviscommon import LatLong +import ivwastrovis + +import numpy as np +import math +from pathlib import Path +from collections.abc import Callable +from enum import Enum + +import time + + +class AxisType(Enum): + SliceIndex = 0 + Dataset = 1 + Velocity = 2 + + +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.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")) + + 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 (arcsec)_ coordinates in arcseconds relative to the star/center +* _Relative (arcsec)_ coordinates in arcseconds relative to the star/center'''), + options=options, selectedIndex=2) + + 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)) + + 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") + + self.basis = ivwbase.properties.BasisProperty("basis", "Basis and offset") + self.information = ivwbase.properties.VolumeInformationProperty( + "information", "Data information") + + self.addProperties([self.filename, self.fitsParameters, + self.latLonCoords, self.axisType, + self.xaxisRange, self.yaxisRange, self.zaxisRange, + self.axisScaling, 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, FITS, Astrophysics"), + help=ivw.unindentMd2doc(FitsVolumeSource.__doc__) + ) + + def getProcessorInfo(self): + return FitsVolumeSource.processorInfo() + + def initializeResources(self): + pass + + @staticmethod + 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) + volumerep = ivw.data.VolumePy(np.copy(data_float, order='C')) + + volume = ivw.data.Volume(volumerep) + volume.swizzlemask = ivw.data.SwizzleMask.defaultData(1) + + datarange = ivw.glm.dvec2(np.nanmin(fits_data.data), np.nanmax(fits_data.data)) + volume.dataMap.dataRange = datarange + volume.dataMap.valueRange = datarange + return volume + + def process(self): + if not self.filename.value: + return + + filename: Path = ivwbase.io.downloadAndCacheIfUrl(self.filename.valueAsString()) + if not filename.exists(): + raise FileNotFoundError(f"Could not locate {self.filename.valueAsString():s}") + + # 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'loading FITS data: {elapsed/1_000_000.0} ms') + + # print(astroviscommon.readFITSData.cache_info()) + + extent_xy, offset_xy, unit = astroviscommon.getLatLongBasis( + fits_data.params, LatLong(self.latLonCoords.value)) + + 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 + # 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", ivw.data.Unit("")) + case AxisType.Dataset: + freq_func = astroviscommon.frequencyFunc(fits_data.params) + freq_range = (freq_func(0), freq_func(dim[0] - 1)) + # 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_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] + axes[2] = ("Velocity", ivw.data.Unit("km/s")) + case _: + raise ValueError(f"Unexpected axis type {self.axisType.value}") + + 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', '')) + + for i, (label, unit) in enumerate(axes): + volume_intensity.axes[i].name = label + volume_intensity.axes[i].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.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 + 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) + self.deserializing = True 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/infravis/astrophysics/src/astrophysicsmodule.cpp b/infravis/astrophysics/src/astrophysicsmodule.cpp new file mode 100644 index 000000000..9f405084b --- /dev/null +++ b/infravis/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/infravis/astrophysics/src/processors/fitsvolumeraycaster.cpp b/infravis/astrophysics/src/processors/fitsvolumeraycaster.cpp new file mode 100644 index 000000000..373b7b6f3 --- /dev/null +++ b/infravis/astrophysics/src/processors/fitsvolumeraycaster.cpp @@ -0,0 +1,70 @@ +/********************************************************************************* + * + * 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/infravis/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp b/infravis/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp new file mode 100644 index 000000000..e94b52100 --- /dev/null +++ b/infravis/astrophysics/src/shadercomponents/fitsnonlineardepthcomponent.cpp @@ -0,0 +1,301 @@ +/********************************************************************************* + * + * 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} + , 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) + .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)} + , 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); +} + +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(); + mat4 m = glm::scale(dataExtent_.get()); + m[3] = vec4{dataOffset_.get(), 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_); + shader.setUniform(buf.replace("{}.validChannelRange", prefix), validChannelRange_); + shader.setUniform("invalidColor", vec4{invalidColor_.get(), invalidAlpha_.get()}); +} + +std::vector> FitsNonlinearDepthComponent::getInports() { + return {{&volumePort_, std::string{"volumes"}}}; +} + +void FitsNonlinearDepthComponent::initializeResources(Shader&) { + // shader definitions... +} + +std::vector FitsNonlinearDepthComponent::getProperties() { + return {&fitsParameters_, &invalidColor_, &invalidAlpha_}; +} + +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] + + ivec2 validChannelRange; // range of valid channels based on line-of-sight velocity +}; + +const float invalidValue = -1.e20; +)"); + +constexpr std::string_view uniforms = util::trim(R"( +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"( +// 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 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 + //coord.xy *= params.assumedDistance; + return params.velocityInf * coord.z / length(coord) + params.velocityStar; +}} + +// 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, 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 || u < 0.0 || u > 1.0 + || v_normalized < 0.0 || v_normalized > 1.0) {{ + return vec4(invalidValue); + }} + + texCoord = vec3(texCoord.st, u); + 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"( +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, minMaxVelocity); +)"); + +constexpr std::string_view unsetColor = util::trim(R"( +color = vec4(0.0); +)"); + +constexpr std::string_view invalidValue = util::trim(R"( +if (volumeVoxel[channel] <= invalidValue) {{ + color = invalidColor; +}} +)"); + +} // 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(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}, + }; + + return segments; +} + +std::optional FitsNonlinearDepthComponent::channelsForVolume() const { + if (auto data = volumePort_.getData()) { + return data->getDataFormat()->getComponents(); + } + return std::nullopt; +} + +} // namespace inviwo diff --git a/infravis/astrophysics/tests/unittests/astrophysics-unittest-main.cpp b/infravis/astrophysics/tests/unittests/astrophysics-unittest-main.cpp new file mode 100644 index 000000000..8ddf58e4e --- /dev/null +++ b/infravis/astrophysics/tests/unittests/astrophysics-unittest-main.cpp @@ -0,0 +1,56 @@ +/********************************************************************************* + * + * 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) { + inviwo::LogCentral::init(); + auto logger = std::make_shared(); + inviwo::LogCentral::getPtr()->setVerbosity(inviwo::LogVerbosity::Error); + inviwo::LogCentral::getPtr()->registerLogger(logger); + + int ret = -1; + { + ::testing::InitGoogleTest(&argc, argv); + inviwo::ConfigurableGTestEventListener::setup(); + ret = RUN_ALL_TESTS(); + } + return ret; +}